Re: [deal.II] How to find the cell index by using global dof index

[Lance Zhang]

Hello Daniel,

thanks for your message.

I know how to find the cell index.



Could I ask one more question?

The question is about the grid.

In code I used the function below to create the grid.

 template 
  Point grid_y_transform (const Point _in)
  {
const double  = pt_in[0];
const double  = pt_in[1];
Point pt_out = pt_in;
return pt_out;
  }

  template 
  void Solid::make_grid()
  {
// Divide the beam, but only along the x- and y-coordinate directions
std::vector< unsigned int > repetitions(dim,
parameters.elements_per_edge);
// Only allow one element through the thickness
// (modelling a plane strain condition)
if (dim == 3){
 repetitions[dim-1] = 8;

//repetitions[0] = 0;
//repetitions[1] = 2;
//repetitions[2] = 2;
}

repetitions[1] = 8;
repetitions[2] = 8;
const Point top_right =Point(0.0, 0.0, -10);
const Point bottom_left =Point(40,20, 10);
GridGenerator::subdivided_hyper_rectangle(triangulation,
  repetitions,
  bottom_left,
  top_right);

// Since we wish to apply a Neumann BC to the right-hand surface, we
// must find the cell faces in this part of the domain and mark them
with
// a distinct boundary ID number. The faces we are looking for are on
the
// +x surface and will get boundary ID 11.
// Dirichlet boundaries exist on the left-hand face of the beam (this
fixed
// boundary will get ID 1) and on the +Z and -Z faces (which correspond
to
// ID 2 and we will use to impose the plane strain condition)
const double tol_boundary = 1e-6;
typename Triangulation::active_cell_iterator cell =
  triangulation.begin_active(), endc = triangulation.end();


for (; cell != endc; ++cell)
  for (unsigned int face = 0;
   face < GeometryInfo::faces_per_cell; ++face)
if (cell->face(face)->at_boundary() == true)
  {
if (std::abs(cell->face(face)->center()[0] - 0.0) <
tol_boundary)
  cell->face(face)->set_boundary_id(1); // -X faces
else if (std::abs(cell->face(face)->center()[0] - 40.0) <
tol_boundary)
  cell->face(face)->set_boundary_id(11); // +X faces
else if (dim == 3 &&
std::abs(std::abs(cell->face(face)->center()[2]) - 10) < tol_boundary)
  cell->face(face)->set_boundary_id(2); // +Z and -Z faces
  }


// Transform the hyper-rectangle into the beam shape
GridTools::transform(_y_transform, triangulation);

GridTools::scale(parameters.scale, triangulation);

vol_reference = GridTools::volume(triangulation);
vol_current = vol_reference;
std::cout << "Grid:\n\t Reference volume: " << vol_reference <<
std::endl;
  }


My question is whether I can use NURBS to create the same grids and how I
could create the grids with the NURBS.


Thanks in advance!

Best regards
Lance

Daniel Arndt  于 2023年12月4日周一 16:08写道：

> Lance,
>
> Have a look at GridTools::vertex_to_cell_map(
> https://www.dealii.org/current/doxygen/deal.II/namespaceGridTools.html#a9b7e2ca8ecd26a472e5225ba91a58acb
> ).
>
> Best,
> Daniel
>
> On Sun, Dec 3, 2023 at 11:38 PM Lance Zhang  wrote:
>
>> Hello team,
>>
>> may I know how to find the cell index with the information of its global
>> Dofs index?
>>
>> I have the global dofs index ,I would like to find which cell this global
>> dof belongs to.
>>
>> Thanks in advance!
>> Best regards
>> Lance
>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> ---
>> You received this message because you are subscribed to the Google Groups
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an
>> email to dealii+unsubscr...@googlegroups.com.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/4453e03d-2511-43c4-a05d-71e695a9276cn%40googlegroups.com
>> <https://groups.google.com/d/msgid/dealii/4453e03d-2511-43c4-a05d-71e695a9276cn%40googlegroups.com?utm_medium=email_source=footer>
>> .
>>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/U7Il_6jPOCQ/unsubscribe.
> To unsubscribe from this gro

[deal.II] How to find the cell index by using global dof index

[Lance Zhang]

Hello team,

may I know how to find the cell index with the information of its global 
Dofs index?

I have the global dofs index ,I would like to find which cell this global 
dof belongs to.

Thanks in advance!
Best regards
Lance

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4453e03d-2511-43c4-a05d-71e695a9276cn%40googlegroups.com.

Re: [deal.II] How to find the neighbor of one cell

[Lance Zhang]

Hello Daniel,

thanks for your message.

I wrote a statement to check if the cell  exists.


 for (const auto  : solid_3d.dof_handler_ref.active_cell_iterators())
  {
//***//
if(cell->at_boundary()){
  continue;
  }//to check whether the cell exists
//**filter if the cell is at boundary ,if yes, the cell 
iterator  will  move forward*//
//***otherwise if the cell exists,the cell index will be 
calculated*//
  const unsigned int cell_index = cell->active_cell_index();

  if (visited_cells[cell_index])
 continue; // Skip if the cell has been visited by its all 
neighbors

  
// Loop over all degrees of freedom on the current cell
for (unsigned int i = 0; i < solid_3d.dofs_per_cell; ++i)

  {
 
// Loop over all neighboring cells
for (unsigned int f = 0; f < cell->n_faces(); ++f)
{
auto neighbor = cell->neighbor(f);//find the neighbor ,but not 
sure if it exists

if(neighbor->at_boundary())
continue;//check if the neighbor of cell exists
int neighbor_index = neighbor->active_cell_index();//get the 
global index of the cell if it exists
Best regards
Lance
On Wednesday, November 29, 2023 at 6:04:45 PM UTC+2 d.arnd...@gmail.com 
wrote:

> Lance,
>
> You would check for the cell being at a boundary first before asking for 
> its neighbor:
>
> // Loop over all neighboring cells
> for (unsigned int f = 0; f < cell->n_faces(); ++f)
> {
> if (cell->at_boundary(f)) 
>   continue;
> auto neighbor = cell->neighbor(f);
> int neighbor_index = neighbor->active_cell_index();//get the 
> global index of the cell if it exists
> [...]
>
>
> Best,
> Daniel
>
> On Wed, Nov 29, 2023 at 9:43 AM Lance Zhang  wrote:
>
>> Hello Wolfgang,
>>
>> thanks for your message.
>>
>> May I know which method or keywords can be used to stand for the 
>> non-existing cell and could I use cell->at_boundary() to check whether 
>> the cell exists or not?
>>
>> I would like to know how we can identify if the object 
>> of 
>> dealii::IteratorRange> 3, false> > >::IteratorOverIterators does not exist. I have tried ==NULL 
>> and nullptr, however, it doe snot work.
>>
>> There is the code snippet below, please take a look:
>>
>>  for (const auto  : solid_3d.dof_handler_ref.active_cell_iterators())
>>   {
>>
>> if(cell->at_boundary()){
>>   continue;
>>   }//to check whether the cell exists
>>
>>   const unsigned int cell_index = cell->active_cell_index();
>>
>>   if (visited_cells[cell_index])
>>  continue; // Skip if the cell has been visited by its all 
>> neighbors
>>   
>> // Loop over all degrees of freedom on the current cell
>> for (unsigned int i = 0; i < solid_3d.dofs_per_cell; ++i)
>>   {
>>  
>> // Loop over all neighboring cells
>> for (unsigned int f = 0; f < cell->n_faces(); ++f)
>> {
>> auto neighbor = cell->neighbor(f);//find the neighbor ,but 
>> not sure if it exists
>>
>> if(neighbor->at_boundary())
>> continue;//check if the neighbor of cell exists
>>     int neighbor_index = neighbor->active_cell_index();//get the 
>> global index of the cell if it exists
>> ...
>>
>> The whole code can be found via the link :
>> https://github.com/Lancelof2019/cook_membrane_snippet/blob/main/code_snippet.cpp
>>
>> Thanks in advance!
>>
>> Best regards
>> Lance
>>
>> Wolfgang Bangerth  于2023年11月28日周二 01:54写道：
>>
>>>
>>> On 11/27/23 14:52, Lance Zhang wrote:
>>> > 
>>> > May I know how to find the neighbor of a cell?
>>> > 
>>> > Here is one part of my code:
>>> > 
>>> --
>>> > for (const auto  : 
>>> solid_3d.get_dof_handler().active_cell_iterators()
>>> > {
>>> >  const unsigned int cell_index = cell->active_cell_index();
>>> > 
>>> >  // Loop over all degrees of freedom on the current cell
&g

Re: [deal.II] How to find the neighbor of one cell

[Lance Zhang]

Hello Wolfgang,

thanks for your message.

May I know which method or keywords can be used to stand for the
non-existing cell and could I use cell->at_boundary() to check whether the
cell exists or not?

I would like to know how we can identify if the object
of dealii::IteratorRange > >::IteratorOverIterators does not exist. I have tried ==NULL
and nullptr, however, it doe snot work.

There is the code snippet below, please take a look:

 for (const auto  : solid_3d.dof_handler_ref.active_cell_iterators())
  {

if(cell->at_boundary()){
  continue;
  }//to check whether the cell exists

  const unsigned int cell_index = cell->active_cell_index();

  if (visited_cells[cell_index])
 continue; // Skip if the cell has been visited by its all
neighbors

// Loop over all degrees of freedom on the current cell
for (unsigned int i = 0; i < solid_3d.dofs_per_cell; ++i)
  {

// Loop over all neighboring cells
for (unsigned int f = 0; f < cell->n_faces(); ++f)
{
auto neighbor = cell->neighbor(f);//find the neighbor ,but not
sure if it exists

if(neighbor->at_boundary())
continue;//check if the neighbor of cell exists
int neighbor_index = neighbor->active_cell_index();//get the
global index of the cell if it exists
...

The whole code can be found via the link :
https://github.com/Lancelof2019/cook_membrane_snippet/blob/main/code_snippet.cpp

Thanks in advance!

Best regards
Lance

Wolfgang Bangerth  于2023年11月28日周二 01:54写道：

>
> On 11/27/23 14:52, Lance Zhang wrote:
> >
> > May I know how to find the neighbor of a cell?
> >
> > Here is one part of my code:
> >
> --
> > for (const auto  :
> solid_3d.get_dof_handler().active_cell_iterators()
> > {
> >  const unsigned int cell_index = cell->active_cell_index();
> >
> >  // Loop over all degrees of freedom on the current cell
> >  for (unsigned int i = 0; i < dofs_per_cell.dofs_per_cell; ++i)
> >  {
> >  // Loop over all neighboring cells
> >  for (const auto  : Find_neighbors?)){
>
> You can write this as
>for (unsigned int f=0; fn_faces(); ++f)
>{
>  DoFHandler::cell_iterator neighbor = cell->neighbor(f);
>  ...
>}
>
> Best
>   W.
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/XYNlOZYskis/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/ef1fdaad-0bde-4478-8054-74df1a572c5f%40colostate.edu
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6P%3Dy7uJAdD5q8DwpUnK-3iQBgTnVh5ioi-Pr62w54D0ikw%40mail.gmail.com.

[deal.II] How to find the neighbor of one cell

[Lance Zhang]

Hello team,

I have one quick question.

May I know how to find the neighbor of a cell?

Here is one part of my code:
--
for (const auto  : solid_3d.get_dof_handler().active_cell_iterators()
{
const unsigned int cell_index = cell->active_cell_index();

// Loop over all degrees of freedom on the current cell
for (unsigned int i = 0; i < dofs_per_cell.dofs_per_cell; ++i)
{
// Loop over all neighboring cells
for (const auto  : Find_neighbors?)){
 ...

-
Thanks in advance!
Best regards
Lance   

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4ae3d66e-faaa-401d-8027-12b6f7602bf0n%40googlegroups.com.

Re: [deal.II] How to convert stiffness matrix in global degrees of freedom into stiffness matrix in global cells

[Lance Zhang]

Hello Wolfgang,

thanks very much for your information.

I got the point.

Best regards
Lance

Wolfgang Bangerth  于 2023年11月13日周一 21:01写道：

>
> On 11/13/23 03:05, Lance Zhang wrote:
> >
> > I wonder whether the stiffness matrix could be converted to a matrix
> > with rows = cell_numbers and columns = cell_numbers.
>
> Lance:
> I think I mentioned this before: You can't. The entries of the stiffness
> matrix are sums over contributions from different cells. So, in 2d with
> Q1 elements for example, every matrix entry A_{ij} has the form
>
>A_{ij} = a+b+c+d
>
> where a...d are the contributions of the four cells adjacent to vertex i
> (or j). Once you add these four numbers together, you can no longer
> determine what the four numbers were. If you need the contributions of
> individual cells, you need to store them before you add them together.
>
> Best
>   W.
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/1sVmdUIFY3k/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/b50df74b-e825-4a66-b6bb-4fbd72aa43f9%40colostate.edu
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6PmNwW2ZimY097wge8DP-jGRV_GBCYEJk1eC_W6Q035cKA%40mail.gmail.com.

Re: [deal.II] How to convert stiffness matrix in global degrees of freedom into stiffness matrix in global cells

[Lance Zhang]

Hello Wolfgang,

thanks for your reply.

I wonder whether the stiffness matrix could be converted to a matrix with
rows = cell_numbers and columns = cell_numbers.

I tried to write the gradient with rows = degree_of_freedom into with rows
= cell_numbers.

###

for (unsigned int i = 0; i < solid_3d.solution_n.size(); ++i) {
solver_vector(i) = solid_3d.solution_n[i];
obj_det(i)= solid_3d.solution_n[i];
norm2_value+=std::pow(solid_3d.solution_n[i],2);
}

  norm2_value=std::sqrt(norm2_value);



   Eigen::VectorXd
solver_gradient=-1*(obj_det/norm2_value).transpose()*inv_dense_matrix*solver_vector*solver_vector.transpose()*100;
   // the length of rows and columns in solver_gradient is respectively
equal to the number of degree of freedom



 /*-new update on degree of finite into
cells */
  //change the length of rows and column into the number of cells//
  Eigen::VectorXd
cell_gradients_solid(solid_3d.triangulation.n_active_cells());
  FE_Q fe_tmp(parameters.poly_degree);
  const QGauss qf_cell(parameters.quad_order);
  const UpdateFlags uf_cell(update_gradients | update_JxW_values);
  FEValues fe_values_ref_tmp(fe_tmp, qf_cell, uf_cell);
  unsigned int cell_index_solid = 0;
  double dof_value = 0;
  double gradient_tmp = 0;
  std::cout<<"***"<is_locally_owned()) {
fe_values_ref_tmp.reinit(cell);

std::vector
local_dof_indices(solid_3d.dofs_per_cell) ;

   cell->get_dof_indices(local_dof_indices);

   for (unsigned int i = 0; i < local_dof_indices.size(); ++i) {

   dealii::types::global_dof_index global_index =
local_dof_indices[i];


   dof_value = solver_gradient(global_index);
   gradient_tmp += dof_value;
   dof_value=0;

 }

 cell_gradients_solid(cell_index_solid) = gradient_tmp;
 gradient_tmp=0;
 ++cell_index_solid;
 }
 }

 //change the length of rows and column into the number of cells
 std::cout< 于2023年10月31日周二 01:32写道：

>
> On 10/30/23 15:54, Lance Zhang wrote:
> >
> > may I know how the stiffness matrix  in global degrees of freedom could
> > be converted into  stiffness matrix in global cells?
>
> Lance:
> How do you define "stiffness matrix in global cells" mathematically?
>
> Best
>   W.
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/1sVmdUIFY3k/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/b56b296b-14bf-46e0-bc39-da2f3b1de3b8%40colostate.edu
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6P%3DRHcfEe63gDG5m23MZXPL8g%2BstczDEeUiZzAqwzFP8Jg%40mail.gmail.com.

[deal.II] How to convert stiffness matrix in global degrees of freedom into stiffness matrix in global cells

[Lance Zhang]

Hello team,

may I know how the stiffness matrix  in global degrees of freedom could be 
converted into  stiffness matrix in global cells?

I would like to use cook_memebrane.cc to get a stiffness matrix in a way of 
global cells.

Currently,I only got the stiffness matrix (solid_3d.tangent_matrix) in 
global degree of freedom.

How could I use cell matrix or stiffness matrix in global degree of freedom 
to reconstruct cells stiffness matrix?

Thanks in advance!

Best regards
Lance

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/80384c58-26f7-426c-b491-47605b38393fn%40googlegroups.com.

Re: [deal.II] Question about how to calculate derivative of stiffness matrix with density

[Lance Zhang]

Hello Wolfgang,

I would like to use this method to get the optimization of a original
object beam to get a shape of optimized object, for example , the shape of
bridge.

Best regards
Lance

Wolfgang Bangerth  于 2023年9月19日周二 23:27写道：

> On 9/19/23 09:21, Lance Zhang wrote:
> >
> > I would like to use adjoint equation to get the gradiant dJ/dΘ.A stands
> for
> > stiffness matrix and Θ stands for density,J stands for destination
> function.
> > Finally I will get the value of sensitivity.
>
> Right. But what's the operation you want to do with dJ/dΘ? Do you actually
> have to have it element by element, or do you want to do a product with
> this, etc?
>
> Best
>   W.
>
>
> --
> 
> Wolfgang Bangerth  email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/d-09gOQBL4A/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/ea732ffc-c612-8e18-fe98-c1088bdc0a8c%40colostate.edu
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6Pnd-jid6t8NONieuEmh59UXkXQTiGsW7v1ggXHM74UNMw%40mail.gmail.com.

[deal.II] Question about how to calculate derivative of stiffness matrix with density

[Lance Zhang]

Dear group,

thanks for your time to review this email.

I have queation about the ways to get derivative of stiffness matrix with
density.

Currently,I have the stiffness matrix A an density vector v.

If I know the constructing process of cell matrix in code,how should I
calculate the derivative of stiffness matrix with density?

A is a 4131*4131 matrix,v is 1024*1 vector. After derivative calculation
the result is a 4131*4131*1024 matrix and a Jacobin matrix,if I am wrong
,please provide your suggestions.

Is there any similar examples written in tutorial in dealii?

Could you please provide any hint or suggestions?

This program is based on cook_membrance.cc
https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc

Through the program the detail of cell matrix with density can be seen.

Best regards
Lance

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6P%3DPoVPixMBUw8Zo%2BJLSmBkyB%2BYNvArQhF7UGeHb4C6NOA%40mail.gmail.com.

[deal.II] Re: How to use cell stiffness matrix to build global stiffness matrix

[Lance Zhang]

Hello Abbas,

Thanks for your reply.

I will take a look at the step13.

Howerver, I would like to know if the cell matrix transform seems like the 
process below:

const Eigen::MatrixXd& cell_matrix,
const Eigen::VectorXd& cell_youngs_modulus,
const std::vector& local_dof_indices,
Eigen::MatrixXd& tangent_matrix,
Eigen::VectorXd& global_youngs_modulus)

   const int N = cell_matrix.rows();
const int M = tangent_matrix.rows();

// Loop over the local degrees of freedom
for (int i = 0; i < N; ++i) {
const int global_i = local_dof_indices[i]; // Global index of the 
i-th local dof

// Add the contribution of the i-th local dof to the global Young's 
modulus vector
global_youngs_modulus(global_i) += cell_youngs_modulus(i);

for (int j = 0; j < N; ++j) {
const int global_j = local_dof_indices[j]; // Global index of 
the j-th local dof

// Add the contribution of the (i,j)-th entry of the cell 
matrix multiplied by the i-th and j-th Young's modulus to the global 
tangent matrix
tangent_matrix(global_i, global_j) += cell_matrix(i, j) * 
cell_youngs_modulus(i) * cell_youngs_modulus(j);
}
}

I am not quite sure if the assembly process for the global stiffness matrix 
is correct.

Best regards
Lance

On Tuesday, August 1, 2023 at 11:47:12 AM UTC+2 abbas.b...@gmail.com wrote:

> Lance, 
>
> I am not sure if I understood you correctly. 
> Maybe looking at the local_assemble_matrix() function in step:13 might 
> help. 
>
> Abbas
>
> On Tuesday, August 1, 2023 at 11:00:10 AM UTC+2 dim...@gmail.com wrote:
>
>> Hello dear group,
>>
>> I have one question which is about how to build global stiffness matrix 
>> with cell stiffness matrix.
>>
>> In our project,the density is a vector (actually we need a mu vector but 
>> in the code of deallii mu is a scalar)with different values,I would like to 
>> use the element of density vector to build cell matrix as written in the 
>> link below:
>>
>>  
>> https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
>> (starting from line 1782)
>>
>> And finally the cell stiffness matrix is used to build global stiffness 
>> matrix,but I don' know how the global stiffness matrix is bulit with cell 
>> stiffness matrix.
>>
>>
>> https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
>> (starting from line 1599)
>> //>>
>>  constraints.distribute_local_to_global(
>> data.cell_matrix, 
>> data.cell_rhs,
>> data.local_dof_indices,
>> tangent_matrix, 
>> system_rhs);
>> }
>>
>>
>> //
>>
>>
>> I have a idea that I use mu[0] to obtain the first cell_matrix and use 
>> mu[1] to get the second cell stiffness matrix ,step by step then I got 
>> final mu[m] to obtain the last cell stiffness matrix and in the end I use 
>> some function like  distribute_local_to_global to construct the global 
>> stiffness matrix,and system rhs and even the solutions.
>>
>>I checked the code in 
>> AffineConstraints::distribute_local_to_global 
>> .
>>  
>> (
>> https://www.dealii.org/current/doxygen/deal.II/affine__constraints_8h_source.html
>> )
>> But the hint information guides me to read the cm.templates.h file.
>> ewcfp I was a little confused ,I did not find the file from dealii lib. 
>> could anyone provide any information or hint? Thanks in advance!
>> Best regards
>> Lance
>> Could 
>>
>> I
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0e042fbd-fab4-408e-a67e-a43f612c7490n%40googlegroups.com.

[deal.II] How to use cell stiffness matrix to build global stiffness matrix

[Lance Zhang]

Hello dear group,

I have one question which is about how to build global stiffness matrix 
with cell stiffness matrix.

In our project,the density is a vector (actually we need a mu vector but in 
the code of deallii mu is a scalar)with different values,I would like to 
use the element of density vector to build cell matrix as written in the 
link below:

 
https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
(starting from line 1782)

And finally the cell stiffness matrix is used to build global stiffness 
matrix,but I don' know how the global stiffness matrix is bulit with cell 
stiffness matrix.

https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
(starting from line 1599)
//>>
 constraints.distribute_local_to_global(
data.cell_matrix, 
data.cell_rhs,
data.local_dof_indices,
tangent_matrix, 
system_rhs);
}


//


I have a idea that I use mu[0] to obtain the first cell_matrix and use 
mu[1] to get the second cell stiffness matrix ,step by step then I got 
final mu[m] to obtain the last cell stiffness matrix and in the end I use 
some function like  distribute_local_to_global to construct the global 
stiffness matrix,and system rhs and even the solutions.

   I checked the code in 
AffineConstraints::distribute_local_to_global 
.
 
(
https://www.dealii.org/current/doxygen/deal.II/affine__constraints_8h_source.html
)
But the hint information guides me to read the cm.templates.h file.
ewcfp I was a little confused ,I did not find the file from dealii lib. 
could anyone provide any information or hint? Thanks in advance!
Best regards
Lance
Could 

I

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e3f16c37-7670-40ee-b0d2-821fdcb7b307n%40googlegroups.com.

Re: [deal.II] Question about how to get the Stiffness matrix partial derivative with respect to density

[Lance Zhang]

Hello Wolfgang,

thanks for your reply.

Currently,I would like to get the sensitivity by using gradient of
objective function with respect to density.

The Objective function is a vector of movement which is calculated from
Ax=b.

J=|x|norm1 is the objective function.
dJ(x,Θ)/dΘ is the gradient that we will obtain in the end.

We use adjoint method written in previous emails to calculate dJ(x,Θ)/dΘ.

Best regards
Lance




Wolfgang Bangerth  于 2023年7月24日周一 21:49写道：

> On 7/24/23 11:24, Lance Zhang wrote:
> >
> > One moire question,how could I get  ∂A(θ)/∂θx,because I did not find any
> > information about density vector like  θ=[θ1,θ2,...,θm],may I know if I
> have
> > to set density vector value in this finite element?
>
> Lance -- I have no idea. I don't know the program, but looking at it, I
> see no
> mention of the word "density" anywhere. Apparently, the person who wrote
> the
> program implemented an equation in which density does not appear. If you
> want
> that the density is part of the formulation, you should first write down
> what
> the specific problem is you want to solve. We can't know what it is you
> want
> to do. You should then compare what you want to do with what the existing
> program does, and change it so that it matches to what you want it to do.
>
> Best
>   W.
>
> --
> 
> Wolfgang Bangerth  email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/5IN4HUV1UhY/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/8819fb69-2fc9-a80c-d37d-a45411d0a89b%40colostate.edu
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6PmT1dR%2B_ni60J-oYC792qDpWsimgX%3D5Gacs%3DJK6j0aftQ%40mail.gmail.com.

Re: [deal.II] Question about how to get the Stiffness matrix partial derivative with respect to density

[Lance Zhang]

Hello Wolfgang,

The value of  A^-1(θ))(-∂A(θ)/∂θ)x is what I would like to obtain.

Thanks in advance!
Best regards
Lance

On Monday, July 24, 2023 at 7:24:53 PM UTC+2 Lance Zhang wrote:

> Hello Wolfgang,
>
> thanks for your reply.
>
> I will follow your point to see if I could find the solution.
>
> One moire question,how could I get  ∂A(θ)/∂θx,because I did not find any 
> information about density vector like  θ=[θ1,θ2,...,θm],may I know if I 
> have to set density vector value in this finite element?
>
> Code Link:
> https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
>
> mu is invoked from prm file at the position of line 252
> Best regards
> Lance
> On Monday, July 24, 2023 at 6:49:12 PM UTC+2 Wolfgang Bangerth wrote:
>
>> > The given conditions are stiffness matrix,destination function and x 
>> from 
>> > Ax=b,x is movement ,b is force,A is global stiffness matrix. 
>> > J(x)=||x||norm-1 
>>
>> This is a questionable choice because dJ/dx is not defined at x=0. 
>>
>> > There is one equation: 
>> > 
>> > dJ(x,θ)/dθ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x 
>> > *∂x/∂θ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x*[A^-1(θ))(-∂A(θ)/∂θ)x],because b has no 
>> relation 
>> > with density. 
>> > 
>> > I don't know how to get [A^-1(θ))(-∂A(θ)/∂θ)x],A is matrix with nxn,if 
>> θ is 
>> > mx1 vector,∂A(θ)/∂θ is nxnxm matrix. 
>>
>> ∂A(θ)/∂θ)x is an nxm matrix. You want to apply A^{-1} to it. The way you 
>> do 
>> that is to call the result 
>> Y = A^{-1} [∂A(θ)/∂θ)x] 
>> and then each column of Y, call it y_k (k=1...m) is computed as 
>> y_k = A^{-1} {[∂A(θ)/∂θ)x]_k} 
>> which requires solving a linear system 
>> A y_k = {[∂A(θ)/∂θ)x]_k} 
>>
>> You never need A^{-1} explicitly. 
>>
>> Best 
>> W. 
>>
>> -- 
>>  
>> Wolfgang Bangerth email: bang...@colostate.edu 
>> www: http://www.math.colostate.edu/~bangerth/ 
>>
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/385ed5b7-2f0e-4343-b3c9-be3c8b5a73d9n%40googlegroups.com.

Re: [deal.II] Question about how to get the Stiffness matrix partial derivative with respect to density

[Lance Zhang]

Hello Wolfgang,

thanks for your reply.

I will follow your point to see if I could find the solution.

One moire question,how could I get  ∂A(θ)/∂θx,because I did not find any 
information about density vector like  θ=[θ1,θ2,...,θm],may I know if I 
have to set density vector value in this finite element?

Code 
Link:https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc

mu is invoked from prm file at the position of line 252
Best regards
Lance
On Monday, July 24, 2023 at 6:49:12 PM UTC+2 Wolfgang Bangerth wrote:

> > The given conditions are stiffness matrix,destination function and x 
> from 
> > Ax=b,x is movement ,b is force,A is global stiffness matrix.
> > J(x)=||x||norm-1
>
> This is a questionable choice because dJ/dx is not defined at x=0.
>
> > There is one equation:
> > 
> > dJ(x,θ)/dθ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x 
> > *∂x/∂θ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x*[A^-1(θ))(-∂A(θ)/∂θ)x],because b has no 
> relation 
> > with density.
> > 
> > I don't know how to get [A^-1(θ))(-∂A(θ)/∂θ)x],A is matrix with nxn,if θ 
> is 
> > mx1 vector,∂A(θ)/∂θ is nxnxm matrix.
>
> ∂A(θ)/∂θ)x is an nxm matrix. You want to apply A^{-1} to it. The way you 
> do 
> that is to call the result
> Y = A^{-1} [∂A(θ)/∂θ)x]
> and then each column of Y, call it y_k (k=1...m) is computed as
> y_k = A^{-1} {[∂A(θ)/∂θ)x]_k}
> which requires solving a linear system
> A y_k = {[∂A(θ)/∂θ)x]_k}
>
> You never need A^{-1} explicitly.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/0a3ade68-cbed-4551-8cba-67e49970217en%40googlegroups.com.

[deal.II] Question about how to get the Stiffness matrix partial derivative with respect to density

[Lance Zhang]

Hello ,

Thanks for your time to see this question.

The given conditions are stiffness matrix,destination function and x from 
Ax=b,x is movement ,b is force,A is global stiffness matrix.
J(x)=||x||norm-1
There is one equation:

dJ(x,θ)/dθ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x 
*∂x/∂θ=∂J(x,θ)/∂θ+∂J(x,θ)/∂x*[A^-1(θ))(-∂A(θ)/∂θ)x],because b has no 
relation with density.

I don't know how to get [A^-1(θ))(-∂A(θ)/∂θ)x],A is matrix with nxn,if θ is 
mx1 vector,∂A(θ)/∂θ is nxnxm matrix. 

But I don't know how to calculate (A^-1(θ))(-∂A(θ)/∂θ),I can simply get the 
value of A^-1.
I use cook_membrane.cc this file for dJ(x,θ)/dθ.

I found the density was just included from the mu invoked from 
parameters.prm.

>


public:
Material_Compressible_Neo_Hook_One_Field(const double mu,
 const double nu)
  :
  kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))),
  c_1(mu / 2.0)
{
  Assert(kappa > 0, ExcInternalError());
}

const double kappa;
const double c_1;

// Value of the volumetric free energy
NumberType
get_Psi_vol(const NumberType _F) const
{
  return (kappa / 4.0) * (det_F*det_F - 1.0 - 2.0*std::log(det_F));
}

NumberType
get_dPsi_vol_dJ(const NumberType _F) const
{
  return (kappa / 2.0) * (det_F - 1.0 / det_F);
}

get_d2Psi_vol_dJ2(const NumberType _F) const
{
  return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F)));
}NumberType
get_Psi_iso(const SymmetricTensor<2,dim,NumberType> _bar) const
{
  return c_1 * (trace(b_bar) - dim);
}
SymmetricTensor<2,dim,NumberType>
get_tau_bar(const SymmetricTensor<2,dim,NumberType> _bar) const
{
  return 2.0 * c_1 * b_bar;
}


 assemble_system_one_cell(const typename 
DoFHandler::active_cell_iterator ,
 ScratchData_ASM ,
 PerTaskData_ASM )
{
  // Due to the C++ specialization rules, we need one more
  // level of indirection in order to define the assembly
  // routine for all different number. The next function call
  // is specialized for each NumberType, but to prevent having
  // to specialize the whole class along with it we have inlined
  // the definition of the other functions that are common to
  // all implementations.
  assemble_system_tangent_residual_one_cell(cell, scratch, data);
  assemble_neumann_contribution_one_cell(cell, scratch, data);
}



<



Unfortunately ,I did not find density vector in this file  to calculate 
-∂A(θ)/∂θ.

Could anyone provide any hint or suggestions?

Thanks in advance!
Best regards
Lance

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f67b2b60-c697-4235-bc57-d0db7890dfccn%40googlegroups.com.

Re: [deal.II] change the value of item modul

[Lance Zhang]

Hello Wolfgang,

The error shows below:

[ 66%] Built target cook_membrane
[100%] Run cook_membrane with Release configuration
///
// Test the GCMMA algorithm
///
terminate called after throwing an instance of 
'dealii::ParameterHandler::ExcNoSubsection'
  what():  

An error occurred in line <1848> of file 
<./source/base/parameter_handler.cc> in function
void dealii::ParameterHandler::scan_line(std::string, const string&, 
unsigned int, bool)
The violated condition was: 
skip_undefined || entries->get_child_optional( 
get_current_full_path(subsection))
Additional information: 
Line <3> of file  On 7/17/23 15:57, Lance Zhang wrote:
> > 
> > The information above is what I wrote,but the errors come out .
>
> Lance -- nobody can help you if you don't say what exactly the error is. 
> Carefully reading the error message often is 50% of the process of 
> figuring 
> out what the problem is.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4dbc0008-6ffb-460d-b0bc-a9aab6e15b35n%40googlegroups.com.

[deal.II] change the value of item modul

[Lance Zhang]

Hello,

I would like to change the value ofset Shear modulus   to 1e6 using a 
parameter in function. so I can control the value of modulus and this 
parameter will be used and changed in each iteration of loop.

There are two ideas
1st)I want to change the value of modulus and the new value will be written 
back to prm file.And in next step the file will read as input to execute 
the next steps.

2nd)I would like to create a object prm and changed the value of 
modulus,but not written to back to  prm file,use prm object in next steps;

in function:
   prm.parse_input("parameters.prm");
  
   prm.enter_subsection("Material properties");

   prm.set("Shear modulus", x[0]);

  prm.leave_subsection();

  
prm.print_parameters("parameters.prm",dealii::ParameterHandler::Text);

 {
  deallog.depth_console(0);

  Parameters::AllParameters parameters("parameters.prm");
  parameters.mu=x[0];
//x_density_percell=parameters.mu;(I found there are errors in 
umpack)
  if (parameters.automatic_differentiation_order == 0)
 {
  std::cout << "Assembly method: Residual and linearisation are 
computed manually." << std::endl;

  // Allow multi-threading
 // Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
   //   
dealii::numbers::invalid_unsigned_int);

  typedef double NumberType;
  Solid solid_3d(parameters);
 solid_3d.run();
   }
 }

The information above is what I wrote,but the errors come out . I am not 
quite familiar with prm part,I am still learning dealii for more knowledge 
and skills.

May I know if anyone could provide any hint or suggestions?Which idea is 
better  for implementation?


Thanks in advance!
Best regards
Lance

#>>
# Listing of Parameters
# -
subsection Assembly method
  # The automatic differentiation order to be used in the assembly of the 
linear system.
  # Order = 0: Both the residual and linearisation are computed manually.
  # Order = 1: The residual is computed manually but the linearisation is 
performed using AD.
  # Order = 2: Both the residual and linearisation are computed using AD. 
  set Automatic differentiation order = 0
end

subsection Finite element system
  # Displacement system polynomial order
  set Polynomial degree = 1

  # Gauss quadrature order
  set Quadrature order  = 2
end


subsection Geometry
  # Number of elements per long edge of the beam
  set Elements per edge  =20 

  # Global grid scaling factor
  set Grid scale = 1e-3
end


subsection Linear solver
  # Linear solver iterations (multiples of the system matrix size)
  set Max iteration multiplier  = 10 

  # Linear solver residual (scaled by residual norm)
  set Residual  = 1e-6

  # Preconditioner type
  set Preconditioner type= ssor

  # Preconditioner relaxation value
  set Preconditioner relaxation  = 0.65

  # Type of solver used to solve the linear system
  set Solver type   = Direct
end


subsection Material properties
  # Poisson's ratio
  set Poisson's ratio = 0.3

  # Shear modulus
  # set Shear modulus   = 0.4225e6
   set Shear modulus   = 1e6
end


subsection Nonlinear solver
  # Number of Newton-Raphson iterations allowed
  set Max iterations Newton-Raphson = 10

  # Displacement error tolerance
  set Tolerance displacement= 1.0e-6

  # Force residual tolerance
  set Tolerance force   = 1.0e-9
end


subsection Time
  # End time
  set End time   = 1

  # Time step size
  set Time step size = 0.1
end

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/b52ea25f-776e-49e0-b303-2c8b66ff8e2bn%40googlegroups.com.

[deal.II] adding openmp in CMakeLists.txt

[Lance Zhang]

Hello everyone,

I have one problem when I tried to add openmp lib to CMakeLists.txt to use 
parallel functionality.

May I know how I could add openmp lib?

And one more question,I wrote two code files which are called gcmma and 
mma,may I know how could I add these two files in this CMakeLists as well?

I tried to find the lib and add them in this txt file,but the error still 
exists.

#
SET(CLEAN_UP_FILES
  # a custom list of globs, e.g. *.log *.vtk
  *.vtk
)

# Usually, you will not need to modify anything beyond this point...

CMAKE_MINIMUM_REQUIRED(VERSION 3.3.0)

include_directories(/usr/include/trilinos)
include_directories(/usr/lib/x86_64-linux-gnu/openmpi/include)
include_directories(/usr/include/suitesparse)
include_directories(/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-real/include)

add_library(mma_gcmma OBJECT gcmma/GCMMASolver.cpp)
add_library(mma_mma OBJECT mma/MMASolver.cpp)

# find Eigen
find_package(Eigen3 REQUIRED)

# include Eigen
include_directories(${EIGEN3_INCLUDE_DIR})



#find_package(OpenMP REQUIRED)
#target_link_libraries(cook_membrane PUBLIC OpenMP::OpenMP_CXX)



FIND_PACKAGE(deal.II 9.3.0 QUIET
  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
  )
IF(NOT ${deal.II_FOUND})
  MESSAGE(FATAL_ERROR "\n"
"*** Could not locate deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to 
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this 
path."
)
ENDIF()
#target_link_libraries(cook_membrane mma_gcmma mma_mma ${DEAL_II_LIBRARIES})

DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()

SET(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.ucd *.d2 *.vtu 
*.pvtu)
MESSAGE(STATUS "deal.II_DIR: ${deal.II_DIR}")
MESSAGE(STATUS "DEAL_II_DIR: ${DEAL_II_DIR}")

#<<<
Could anyone provide any hint or suggestions?

Thanks
Best regards
Lance


 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5f9f021d-6dbb-47e6-84d6-bd9542146da1n%40googlegroups.com.

[deal.II] Density of rectangle and objective function in code file cook_membrane.cc

[Lance Zhang]

Hello,

If I would like to get the density of the rectangle,how should I calculate 
the value from rectangle in code cook_membrane.cc?

How should I design the function J(x,θ) to get the objective function J(x,θ
),?This function indicates the minimum movement.θ is density of this 
rectangle object. 
>>
Conditions:Ax=b is given,A,b are known, solver x can be calculated.
I would like to get the vector of dJ/dθ=∂J/∂θ+(∂J/∂x.A^-1(θ) ){∂b(θ)/∂θ)-∂A(
θ)/∂θ}.
The expression is calculated with adjoint method.
>>  
However,I didn't find the value of density in this code file ,or did I miss 
something?

I wish someone could provide hint or suggestions.

Thanks in advance！
Best regards
Lance

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/41b88072-30b1-484f-8859-c46f498ceb10n%40googlegroups.com.

[deal.II] Re: finite difference method to calculate the gradient of density in cells by reference file of cook_membrane.cc

[Lance Zhang]

Hello,
 
later , I tried to use the number of active cells,
density.reinit(triangulation.n_active_cells());
gradient.reinit(triangulation.n_active_cells());

But I still got the same issue.

In last code pasted in the first email I think I have to replace 
dof_handler_ref.n_dofs() with triangulation.n_active_cells(),because 
dof_handler_ref.n_dofs()  is about the number of dofs instead of activae 
cells.


Could anyone provide any suggestions?

Best regards
Lance
On Sunday, July 2, 2023 at 2:36:50 PM UTC+2 Lance Zhang wrote:

> Hello ,
>
> I would like to use finite difference method to calculate the gradient of 
> density in cells.
>
> But I got issue when the program was processed at the part of 
> compute_gradient();
>
> Here is the code below  about density and gradient.
>
>
> //>>>>>>>>>>>>>>>Density>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
> template 
>
> void Solid::compute_density()
> {
>   typename DoFHandler::active_cell_iterator cell = 
> dof_handler_ref.begin_active(),
>  endc = 
> dof_handler_ref.end();
>
>   for (; cell != endc; ++cell)
>   {
> // Compute the grid density for each cell,
> density[cell->active_cell_index()] = 1.0 / cell->measure();
>   }
> }
>
>
>
> //>>>>>>>>>>>>>>>>>>>Gradient>>>>>>>>>>>>>>>>>>>>>>>>
> template 
>
>
> void Solid::compute_gradient()
> {
>   typename DoFHandler::active_cell_iterator cell = 
> dof_handler_ref.begin_active(),
>  endc = 
> dof_handler_ref.end();
>
>   for (; cell != endc; ++cell)
>   {
> for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; 
> ++face)
> {
>   if (cell->face(face)->at_boundary())
>   {
> // Compute the gradient at each boundary face by differencing the 
> neighboring cell densities
> const typename DoFHandler::cell_iterator neighbor = 
> cell->neighbor(face);
> gradient[cell->face(face)->index()] = 
> density[neighbor->active_cell_index()] - density[cell->active_cell_index()];
>   }
> }
>   }
> }
> //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>
> The error comes out at the part of compute_gradient().
>
> 
> Error information:
> Thread 1 "cook_membrane" received signal SIGSEGV, Segmentation fault.
> 0x74bfcdaa in dealii::CellAccessor<3, 3>::active_cell_index() 
> const () from /lib/x86_64-linux-gnu/libdeal.ii.so.9.3.2
>
> 
> I tried to reduce the numbers of grids and dofs ,but I sill got the same 
> issue.
>
> in the part of void Solid::run(),
> I add  the code below  
>
>  template 
>   void Solid::run(){
>  .//other codes
>  .//other codes
>  .//other codes
> density.reinit(dof_handler_ref.n_dofs());
> gradient.reinit(dof_handler_ref.n_dofs());
> //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
> print_vertical_tip_displacement();
>
> std::cout<<"Density claculation"< compute_density();//sucessfully
> std::cout<<"Gradient calculation"< compute_gradient();//error comes out
> }
>
> Could anyone provide any hint or suggestions?
>
> Thanks in advance!
> Best regards
> Lance
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/989423f2-b29d-4199-82c5-6a121b120d01n%40googlegroups.com.

[deal.II] finite difference method to calculate the gradient of density in cells by reference file of cook_membrane.cc

[Lance Zhang]

Hello ,

I would like to use finite difference method to calculate the gradient of 
density in cells.

But I got issue when the program was processed at the part of 
compute_gradient();

Here is the code below  about density and gradient.


//>>>Density>>>
template 

void Solid::compute_density()
{
  typename DoFHandler::active_cell_iterator cell = 
dof_handler_ref.begin_active(),
 endc = 
dof_handler_ref.end();

  for (; cell != endc; ++cell)
  {
// Compute the grid density for each cell,
density[cell->active_cell_index()] = 1.0 / cell->measure();
  }
}



//>>>Gradient
template 


void Solid::compute_gradient()
{
  typename DoFHandler::active_cell_iterator cell = 
dof_handler_ref.begin_active(),
 endc = 
dof_handler_ref.end();

  for (; cell != endc; ++cell)
  {
for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; 
++face)
{
  if (cell->face(face)->at_boundary())
  {
// Compute the gradient at each boundary face by differencing the 
neighboring cell densities
const typename DoFHandler::cell_iterator neighbor = 
cell->neighbor(face);
gradient[cell->face(face)->index()] = 
density[neighbor->active_cell_index()] - density[cell->active_cell_index()];
  }
}
  }
}
//<<<

The error comes out at the part of compute_gradient().

Error information:
Thread 1 "cook_membrane" received signal SIGSEGV, Segmentation fault.
0x74bfcdaa in dealii::CellAccessor<3, 3>::active_cell_index() const 
() from /lib/x86_64-linux-gnu/libdeal.ii.so.9.3.2

I tried to reduce the numbers of grids and dofs ,but I sill got the same 
issue.

in the part of void Solid::run(),
I add  the code below  

 template 
  void Solid::run(){
 .//other codes
 .//other codes
 .//other codes
density.reinit(dof_handler_ref.n_dofs());
gradient.reinit(dof_handler_ref.n_dofs());
//<<
print_vertical_tip_displacement();

std::cout<<"Density claculation"

Re: [deal.II] Control points are missing in changed code by reference of cook_membrance.cc

[Lance Zhang]

Hello  Simranjeet,

thanks for your reply.

I got the solution to the issue.

Best regards
Lance

Simranjeet Singh  于2023年7月2日周日 09:37写道：

> Hi Lance,
> I am not sure if I understand your question correctly, but if you're
> asking "Why is there only one division along the z-axis", I believe that is
> what you're essentially doing with your code-
> "1099 // Divide the beam, but only along the x- and y-coordinate
> directions"
> ...
> "1101 // Only allow one element through the thickness"
>
> Also, in this case, the point cloud looks fine, there are 8 points
> defining one cell -> six surfaces per cell.
> You have *32*32*1* cells in total.
> Regards,
> Simranjeet
>
>
> On Saturday, July 1, 2023 at 10:22:25 PM UTC+5:30 dim...@gmail.com wrote:
>
>> Hello Wolfgang,
>>
>> thanks for your reply.
>>
>> One quick question,is the issue about cells in zox(only 32 cells) surface
>> related to the paraview? (in xoy surface there is 32*32 surfaces)
>>
>> Best regards
>> Lance
>>
>> Wolfgang Bangerth  于 2023年7月1日周六 18:43写道：
>>
>>> On 7/1/23 10:00, Lance Zhang wrote:
>>> >
>>> > I used Paraview to review this rectangle and use  the item of
>>> > Representation->Points to display these points.
>>>
>>> Lance -- I do not know what Paraview does to select and visualize these
>>> points. As you already found out, all of the cells are there when you
>>> visualize the mesh. You'll have to figure out how Paraview defines the
>>> points
>>> it shows when you do Representation->Points -- but that's a Paraview
>>> question.
>>>
>>> Best
>>>   W>
>>>
>>> --
>>> 
>>> Wolfgang Bangerth  email: bang...@colostate.edu
>>> www:
>>> http://www.math.colostate.edu/~bangerth/
>>>
>>>
>>> --
>>> The deal.II project is located at http://www.dealii.org/
>>> For mailing list/forum options, see
>>> https://groups.google.com/d/forum/dealii?hl=en
>>> ---
>>> You received this message because you are subscribed to a topic in the
>>> Google Groups "deal.II User Group" group.
>>> To unsubscribe from this topic, visit
>>> https://groups.google.com/d/topic/dealii/2YU5TDG9WfQ/unsubscribe.
>>> To unsubscribe from this group and all its topics, send an email to
>>> dealii+un...@googlegroups.com.
>>>
>> To view this discussion on the web visit
>>> https://groups.google.com/d/msgid/dealii/e2afd181-6f92-ed6b-a1be-9c6b60b77078%40colostate.edu
>>> .
>>>
>> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/2YU5TDG9WfQ/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/c0aa103d-be31-4bea-8fad-3f604bad1e8an%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/c0aa103d-be31-4bea-8fad-3f604bad1e8an%40googlegroups.com?utm_medium=email_source=footer>
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6P%3D2xxVmXzrrZmy%2Bymqafe_gRUC6UeFPLRzHk7eGMvKyPw%40mail.gmail.com.

Re: [deal.II] Control points are missing in changed code by reference of cook_membrance.cc

[Lance Zhang]

Hello Wolfgang,

thanks for your reply.

One quick question,is the issue about cells in zox(only 32 cells) surface
related to the paraview? (in xoy surface there is 32*32 surfaces)

Best regards
Lance

Wolfgang Bangerth  于 2023年7月1日周六 18:43写道：

> On 7/1/23 10:00, Lance Zhang wrote:
> >
> > I used Paraview to review this rectangle and use  the item of
> > Representation->Points to display these points.
>
> Lance -- I do not know what Paraview does to select and visualize these
> points. As you already found out, all of the cells are there when you
> visualize the mesh. You'll have to figure out how Paraview defines the
> points
> it shows when you do Representation->Points -- but that's a Paraview
> question.
>
> Best
>   W>
>
> --
> 
> Wolfgang Bangerth  email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/2YU5TDG9WfQ/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/e2afd181-6f92-ed6b-a1be-9c6b60b77078%40colostate.edu
> .
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CADwW6Pnu7RWjCc5q7oq86-TJzte5KNK%2BMbWGUkaD6fL6Mh1CGQ%40mail.gmail.com.

[deal.II] Re: Control points are missing in changed code by reference of cook_membrance.cc

[Lance Zhang]

each surface has 32*32 cells.

On Saturday, July 1, 2023 at 12:10:06 PM UTC+2 Lance Zhang wrote:

> Hello,
>
> I updated some parts in cook_membrance.cc to fulfill requests.
>
> But after updating I have one issue.Normally the rectangle has 6 surfaces 
> and control points on each surface as while.In my updated code,I can only 
> see control point in two surfaces.
>
> May I know which parts in code I should change?
>
> [image: image_2023-07-01_120629896.png]
> [image: Screenshot from 2023-07-01 11-57-07.png]
>
> Here is the code from cc file.
> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
> 1096   template 
> 1097   void Solid::make_grid()
> 1098   {
> 1099 // Divide the beam, but only along the x- and y-coordinate 
> directions
> 1100 std::vector< unsigned int > repetitions(dim, 
> parameters.elements_per_edge);
> 1101 // Only allow one element through the thickness
> 1102 // (modelling a plane strain condition)
> 1103 if (dim == 3){
> 1104  repetitions[dim-1] = 1;
> 1105 
> 1106 //repetitions[0] = 0;
> 1107 //repetitions[1] = 0;
> 1108 //repetitions[2] = 1;
> 1109 }
> 1110 //const Point bottom_left = (dim == 3 ? Point(0.0, 0.0, 
> -0.5) : Point(0.0, 0.0));
>    //  const Point top_right = (dim == 3 ? Point(48.0, 44.0, 
> 0.5) : Point(48.0, 44.0));
> 1112 const Point top_right =Point(0.0, 0.0, -10);
> 1113 const Point bottom_left =Point(40,20, 10);
> 1114 GridGenerator::subdivided_hyper_rectangle(triangulation,
> 1115   repetitions,
> 1116   bottom_left,
> 1117   top_right);
> 1118 
> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>
> Can anyone provide any hint or suggestions?
> Thanks in advance!
> Lance
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/31e05843-01ee-4db0-a04a-da0215485c5en%40googlegroups.com.

Re: [deal.II] inverse matrix of matrix with dealii::BlockSparseMatrix

[Lance Zhang]

Hello  Wolfgang,

thanks for your message.

May I know if there is any method which can be used for  adjoint 
sensitivity method?

Best regards
Lance
On Thursday, June 22, 2023 at 7:29:18 PM UTC+2 Wolfgang Bangerth wrote:

> On 6/22/23 06:35, Lance Zhang wrote:
> > 
> > May I know whether there is any other method to change the matrix with 
> data 
> > type BlockSparseMatrix into a matrix with 
> BlockSparseMatrix?
>
> No. The inverse of a sparse matrix is, in general, not sparse. For all 
> problems of a realistic size, even if you can store a matrix A, you will 
> not 
> be able to store A^{-1}. As a consequence, we never use the inverse of a 
> sparse matrix in deal.II, and there are no functions to compute inverses 
> of 
> sparse matrices.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/7de63c70-15c5-4e2f-9b79-ce8770970889n%40googlegroups.com.

[deal.II] Re: matrix multiply matrix

[Lance Zhang]

Hello ,

I changed the code into :

  BlockSparseMatrix output_vector;
 SparseDirectUMFPACK A_direct;

 
A_direct.initialize(solid_3d.tangent_matrix.block(solid_3d.u_dof,solid_3d.u_dof));
 A_direct.vmult(output_vector.block(solid_3d.u_dof), 
solid_3d.system_rhs.block(solid_3d.u_dof));

But I got error information:

  | ^~~
/usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc: In 
function ‘int main(int, char**)’:
/usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc:2510:44: 
error: no matching function for call to 
‘dealii::BlockSparseMatrix::block(Cook_Membrane::Solid<3, 
double>::)’
 2510 |  A_direct.vmult(output_vector.block(solid_3d.u_dof), 
solid_3d.system_rhs.block(solid_3d.u_dof));
  | ~~~^~~~
In file included from /usr/include/deal.II/lac/block_sparse_matrix.h:22,
 from /usr/include/deal.II/lac/sparse_direct.h:25,
 from 
/usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc:38:
/usr/include/deal.II/lac/block_matrix_base.h:1567:1: note: candidate: 
‘dealii::BlockMatrixBase::BlockType& 
dealii::BlockMatrixBase::block(unsigned int, unsigned int) 
[with MatrixType = dealii::SparseMatrix; 
dealii::BlockMatrixBase::BlockType = 
dealii::SparseMatrix


I don't understand why the error comes out,the similar code comes out in 
line of 2180 in cook_membrane.cc.

BlockVector _update
 SparseDirectUMFPACK A_direct;
  A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
  A_direct.vmult(newton_update.block(u_dof), 
system_rhs.block(u_dof));

solid3d.u_dof=0 which I got by writing in main function.

Could anyone provide any hint or suggestions?

Thanks in advance!
Best regards
Lance
On Thursday, June 22, 2023 at 6:01:44 PM UTC+2 Lance Zhang wrote:

> the code is related to cook_membrane.cc.
>
> On Thursday, June 22, 2023 at 6:00:55 PM UTC+2 Lance Zhang wrote:
>
>> Hello ,
>>
>> When I use vmult to multiply matrix with vector,I found error information:
>> Error information:
>>  BOOST_HEADER_DEPRECATED("")
>>   | ^~~
>> /usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc: 
>> In function ‘int main(int, char**)’:
>> /usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc:2510:24:
>>  
>> error: no matching function for call to 
>> ‘dealii::SparseDirectUMFPACK::vmult(dealii::BlockSparseMatrix&, 
>> dealii::BlockVector&)’
>>  2510 |  A_direct.vmult(output_vector,solid_3d.system_rhs);
>>   |  ~~^~~
>>
>> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
>> code:
>>   
>>  BlockSparseMatrix output_vector;
>>  SparseDirectUMFPACK A_direct;
>>  A_direct.initialize(solid_3d.tangent_matrix);
>>  A_direct.vmult(output_vector,solid_3d.system_rhs);
>>
>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>
>> BlockSparseMatrix is the  data type of solid_3d.tangent_matrix
>> The data type of system_rhs is BlockSparseMatrix .
>>
>> Could anyone provide any hint or sugesstion?
>>
>> Thanks in advance!
>> Best regards
>> Lance
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4403f4c0-9110-4c5f-a19b-505d0e3e0494n%40googlegroups.com.

[deal.II] Re: matrix multiply matrix

[Lance Zhang]

the code is related to cook_membrane.cc.

On Thursday, June 22, 2023 at 6:00:55 PM UTC+2 Lance Zhang wrote:

> Hello ,
>
> When I use vmult to multiply matrix with vector,I found error information:
> Error information:
>  BOOST_HEADER_DEPRECATED("")
>   | ^~~
> /usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc: In 
> function ‘int main(int, char**)’:
> /usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc:2510:24:
>  
> error: no matching function for call to 
> ‘dealii::SparseDirectUMFPACK::vmult(dealii::BlockSparseMatrix&, 
> dealii::BlockVector&)’
>  2510 |  A_direct.vmult(output_vector,solid_3d.system_rhs);
>   |  ~~^~~
>
> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
> code:
>   
>  BlockSparseMatrix output_vector;
>  SparseDirectUMFPACK A_direct;
>  A_direct.initialize(solid_3d.tangent_matrix);
>  A_direct.vmult(output_vector,solid_3d.system_rhs);
>
> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>
> BlockSparseMatrix is the  data type of solid_3d.tangent_matrix
> The data type of system_rhs is BlockSparseMatrix .
>
> Could anyone provide any hint or sugesstion?
>
> Thanks in advance!
> Best regards
> Lance
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/69565fa3-c884-4716-870f-3ff03359dab1n%40googlegroups.com.

[deal.II] matrix multiply matrix

[Lance Zhang]

Hello ,

When I use vmult to multiply matrix with vector,I found error information:
Error information:
 BOOST_HEADER_DEPRECATED("")
  | ^~~
/usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc: In 
function ‘int main(int, char**)’:
/usr/share/doc/libdeal.ii-doc/examples/cook_membrane3/cook_membrane.cc:2510:24: 
error: no matching function for call to 
‘dealii::SparseDirectUMFPACK::vmult(dealii::BlockSparseMatrix&, 
dealii::BlockVector&)’
 2510 |  A_direct.vmult(output_vector,solid_3d.system_rhs);
  |  ~~^~~


code:
  
 BlockSparseMatrix output_vector;
 SparseDirectUMFPACK A_direct;
 A_direct.initialize(solid_3d.tangent_matrix);
 A_direct.vmult(output_vector,solid_3d.system_rhs);

<<<

BlockSparseMatrix is the  data type of solid_3d.tangent_matrix
The data type of system_rhs is BlockSparseMatrix .

Could anyone provide any hint or sugesstion?

Thanks in advance!
Best regards
Lance

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/858ef3d1-651a-4a5d-9258-91effd184ac3n%40googlegroups.com.

[deal.II] Re: inverse matrix of matrix with dealii::BlockSparseMatrix in file cook_membrane.cc

[Lance Zhang]

Hello,

The code is related to the file cook_membrane.cc.


I added the content below in main function
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
double **get_matrix=new double*[get_rows];
for(int i=0;i 
>::operator() (j=24, i=0, this=0x7fffdde8) at 
/usr/include/deal.II/lac/block_matrix_base.h:2112
2112  col_index.second);
(gdb) exit
A debugging session is active.

Inferior 1 [process 113432] will be killed.


But from the process information it seems the compiling is successful  and 
some result can be displayed.

In code,I wrote the content 2d array to obtain the value from 
solid_3d.tangent_matrix for inverse matrix:
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

  std::cout<< "L1 norm:"<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Could anyone provide any suggestion or hint?

Thanks in advance!
Lance
On Thursday, June 22, 2023 at 2:35:52 PM UTC+2 Lance Zhang wrote:

> Hello,
>
> Currently, I met one problem when I calculate the inverse matrix with 
> datatype of class dealii::BlockSparseMatrix.
>
> I searched the website of dealii for the method,but it seems there is no 
> method to get the inverse matrix in BlockSparseMatrix.
>
> May I know whether there is any other method to change the matrix with 
> data type BlockSparseMatrix into a matrix with 
> BlockSparseMatrix?
>
> Best regards
> Lance
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/4cae4027-81b5-439c-a2f3-2c63bfd23ee2n%40googlegroups.com.

[deal.II] inverse matrix of matrix with dealii::BlockSparseMatrix

[Lance Zhang]

Hello,

Currently, I met one problem when I calculate the inverse matrix with 
datatype of class dealii::BlockSparseMatrix.

I searched the website of dealii for the method,but it seems there is no 
method to get the inverse matrix in BlockSparseMatrix.

May I know whether there is any other method to change the matrix with data 
type BlockSparseMatrix into a matrix with BlockSparseMatrix?

Best regards
Lance

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/966846af-5a12-44ee-b651-8c94efaed98fn%40googlegroups.com.

Re: [deal.II] Questions about cook_membrane.cc

[Lance Zhang]

Hello  Wolfgang,

thanks for your message.

Currently,I tried to implement adding local stiffness matrix into global  
stiffness matrix.

1.May I know which parameter can control the index of cells and add 
the local stiffness matrix of cell into global matrix?

2.Whether is each cell matrix a kind of FullMatrix   and 
local_dof_indices a local index parameter?
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  template 
  struct Assembler_Base
  {
virtual ~Assembler_Base() {}

// Here we deal with the tangent matrix assembly structures. The
// PerTaskData object stores local contributions.
struct PerTaskData_ASM
{
  const Solid  *solid;
  FullMatrix   cell_matrix;
  Vector   cell_rhs;
  std::vector local_dof_indices;

  PerTaskData_ASM(const Solid *solid)
:
solid (solid),
cell_matrix(solid->dofs_per_cell, solid->dofs_per_cell),
cell_rhs(solid->dofs_per_cell),
local_dof_indices(solid->dofs_per_cell)
  {}
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
I will appreciate it if you could provide any hint or suggestions.
I wish I could receive your reply.

Best regards
Lance
On Monday, June 12, 2023 at 12:43:40 AM UTC+2 Wolfgang Bangerth wrote:

> On 6/11/23 05:00, Lance Zhang wrote:
> > 
> > If I want to describe a deformation when the rectangular cuboid object 
> > receives a force F , how could I control the mesh moving during the 
> force 
> > deformation ? Should I use cell or node in the surface of rectangular 
> cuboid 
> > to represent the displacement?I s there any function in dealII for mesh 
> control?
>
> Lance:
> You will want to read step-18.
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/ba64ffc1-edca-4164-8af0-530b006d99fdn%40googlegroups.com.

Re: [deal.II] Questions about cook_membrane.cc

[Lance Zhang]

Hello  Daniel,

thanks for your message.

With the link I have created  rectangular cuboid and using 
subdivided_hyper_rectangle() I got points rectangular cuboid in as well.

But I still have one question.

If I want to describe a deformation when the rectangular cuboid  object 
receives a force F , how could I control the mesh moving during the force 
deformation ? Should I use cell or node in the surface of  rectangular 
cuboid to represent the displacement?I s there any function in dealII for 
mesh control?

Thanks in advance!
Best regards
Lance


On Friday, June 9, 2023 at 3:57:08 PM UTC+2 d.arnd...@gmail.com wrote:

> Lance,
>
> Have a look at the documentation for 
> GridGenerator::subdivided_hyper_recyangle (
> https://www.dealii.org/current/doxygen/deal.II/namespaceGridGenerator.html#ac76417d7404b75cf53c732f456e6e971
> ).
> You need to specify the coordinates of the corner with the lowest 
> coordinates in all space dimensions, the corner with the highest 
> coordinates in all space dimensions, as well as the subdivisions for all 
> space dimensions.
>
> Best,
> Daniel
>
> On Fri, Jun 9, 2023 at 5:27 AM Lance Zhang  wrote:
>
>> Hello,
>>
>> I am a novice in dealII and I am interested in dealII.
>>
>> Currently, I am trying to use the code of cook_membrane.cc to create  a 
>> cantilever.
>> Please take a look at the link for your reference.
>>
>> https://www.dealii.org/current/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html
>>
>> 1)May I know how I could create a Rectangular cuboid like the object a) 
>> below? and then change the Rectangular  into a mesh?
>>
>>
>> [image: image_2023-06-09_111218320.png]
>>
>> 2)I tried to configure the parameter of const Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html> 
>> top_right = (dim == 3 ? Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html>(48.0, 
>> 44.0, 0.5) : Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html>(48.0,
>>  
>> 44.0)) ,but I got failed information.
>>
>>
>>
>> template 
>> void Solid::make_grid()
>> {
>>
>> Divide the beam, but only along the x- and y-coordinate directions
>> std::vector< unsigned int > repetitions(dim, 
>> parameters.elements_per_edge);
>>
>> Only allow one element through the thickness (modelling a plane strain 
>> condition)
>> if (dim == 3)
>> repetitions[dim-1] = 1;
>> const Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html> 
>> bottom_left = (dim == 3 ? Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html>(0.0, 
>> 0.0, -0.5) : Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html>(0.0, 
>> 0.0));
>> const Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html> 
>> top_right = (dim == 3 ? Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html>(48.0, 
>> 44.0, 0.5) : Point 
>> <https://www.dealii.org/developer/doxygen/deal.II/classPoint.html>(48.0,
>>  
>> 44.0));
>> GridGenerator::subdivided_hyper_rectangle 
>> <https://www.dealii.org/developer/doxygen/deal.II/namespaceGridGenerator.html#ac76417d7404b75cf53c732f456e6e971>
>> (triangulation 
>> <https://www.dealii.org/developer/doxygen/deal.II/p4est__wrappers_8cc.html#ace00f2f80d9780ef9aa1007e1c22c6a4>
>> ,
>> repetitions,
>> bottom_left,
>> top_right);
>> [image: image_2023-06-09_112028575.png]n In the code lines:
>> 1067  template 
>> 1068   void Solid::make_grid()
>> 1069   {
>> 1070 // Divide the beam, but only along the x- and y-coordinate 
>> directions
>> 1071 std::vector< unsigned int > repetitions(dim, 
>> parameters.elements_per_edge);
>> 1072 // Only allow one element through the thickness
>> 1073 // (modelling a plane strain condition)
>> 1074 if (dim == 3)
>> 1075   repetitions[dim-1] = 1;
>> 1076 
>> 1077 const Point bottom_left = (dim == 3 ? Point(0.0, 0.0, 
>> -0.9) : Point(0.0, 0.0));
>> 1078 const Point top_right = (dim == 3 ? Point(48.0, 44.0, 
>> 20) : Point(48.0, 44.0));
>>
>>
>> Could anyone provide any hint or suggestions?
>>
>> Thank you in advance!
>>
>> Best regards
>> Lance
>>
>> -- 
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see 
>> https://groups.google.com/d/forum/