Assembly of material forces

[Seyed Ali Mohseni]

Dear deal.II experts,

I have a question regarding the assembly of material forces also known as configurational forces. 

First of all, I have constructed my own solid mechanics code within deal.II which I already validated with our own C++ code. Hence, the deal.II code works until here.

Attached you find my results on material forces for one single element. It also works perfectly. 

Now I try to compute the material forces for the whole structure. Hence, I tried using the same procedure of assembly for our LHS or RHS matrices in deal.II, but it fails. Do I have to take a special approach in order to make it work? I assume my local material forces have to be assembled to the global nodes, but this assembly is it just a simple mapping to the nodes? Since there are quadtree constraints, it will not be that simple, I think.

The successful results are computed within our simple c++ code and I just show them for comparison reasons. 

Extra question: Would my work be of interest for your code gallery maybe in the near future?  

Kind regards,

S. A. Mohseni

[Wolfgang Bangerth]

> Now I try to compute the material forces for the whole structure. Hence, I
> tried using the same procedure of assembly for our LHS or RHS matrices in
> deal.II, but it fails.

"It fails" does not carry enough information to really say what may be the
cause. Please take a look at

https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#im-stuck

> Extra question: Would my work be of interest for your code gallery maybe in
> the near future?

If you have good questions, always!

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Seyed Ali Mohseni]

Dear Prof. Bangerth,

It fails, because I tried to use the same procedure for the standard assembly such as in step-18 or step-40 and I received some awkward contour plots which won't match with what I showed from the C++ code. Hence, fails is refered to my own failure in understanding the core structure of deal.II ;)

Can I assemble a vector such as my configurational force vector (8x1) after I solved the system? Hence, I am thinking about a procedure such as:

configurational_forces = cell_cf;

cell->get_dof_indices(local_dof_indices);

constraints.distribute_local_to_global(cell_cf, local_dof_indices, configurational_forces);

where configurational_forces represents a Vector<double> variable I defined globally within my class.

I also found some MeshWorker approach such as shown in step-12. Is it also possible to use such a method to reach my aim?

Conclusion: It would be great, if you can at least guide me to the correct direction of how to assemble a vector after I solved my system. Assembly of the cell vector should also take all existing boundary conditions, constraints etc. into account.

Thank you :)

[Daniel Arndt]

Seyed,

What exactly are you trying to do? Which equations are you considering? How would what you want to do look in formulas?
The information you provide are much too high-level to expect a useful answer.

[...]

Conclusion: It would be great, if you can at least guide me to the correct direction of how to assemble a vector after I solved my system. Assembly of the cell vector should also take all existing boundary conditions, constraints etc. into account.

What exactly do you mean by that?

Best,
Daniel

[Seyed Ali Mohseni]

My problem is not the theory. I am just not a deal.II expert and therefore it is still quite difficult for me to implement my approach. Going in to theoretical detail would take maybe much longer. Maybe it is better instead, if I reformulate my question like this:

There is a vector called cell_cf for instance. It contains all integrated material forces (summed up material forces from integration points) for each cell. Each cell has 4 vertices and 8 DoFs in 2D. Hence, it has to be a 8x1 vector.

For instance:

CONFIGURATIONAL FORCES WITHIN ELEMENT

   0.010126  

   0.007502  

   0.011509  

   -0.012249  

   -0.008940  

   -0.009762  

   -0.012695  

   0.014508 

Now imagine you have this vector which represents forces for each vertex, e.g. vertex 1 has 0.010126 in X-direction and 0.007502 in Y-direction, vertex 2 has  0.011509 in X-direction and -0.012249  in Y-direction and so on...

How is one able now to take these vectors for each cell and assemble them correspondent to the DoF and summarize the values for coinciding nodes. Remember the values are computed within Gauss points and are just summarized over integration points to one single vector for each cell. It means there should be nothing stored on the nodes yet, but we have to take the information to the global structure and node/vertices level.

I already tried to take the approach we implemented within our C++ code, but it is quite difficult to store data in deal.II since there is no processinfo or such a thing which tracks all data created and stored over the computation. I cannot save my configurational force vector and access it anywhere except I store it by using the setup_quadrature_point_history()  and update_quadrature_point_history() functions from step-18. But this means I have to store my elastic energy then update and compute my material forces which are based on the energy and update again. This is not very nice to call a function several times for every little thing you want to save.

Isn't there a better way to store data efficiently during runtime?

Sorry for causing headaches ^^'

[Daniel Arndt]

Seyed,

There is a vector called cell_cf for instance. It contains all integrated material forces (summed up material forces from integration points) for each cell. Each cell has 4 vertices and 8 DoFs in 2D. Hence, it has to be a 8x1 vector.

[...]

Now imagine you have this vector which represents forces for each vertex, e.g. vertex 1 has 0.010126 in X-direction and 0.007502 in Y-direction, vertex 2 has  0.011509 in X-direction and -0.012249  in Y-direction and so on...

So you have a FESystem<2>(FE_Q<2>(1),2) element, i.e. linear continuous elements with two components, and the vector represents corresponding local DoF values?
 

How is one able now to take these vectors for each cell and assemble them correspondent to the DoF and summarize the values for coinciding nodes. Remember the values are computed within Gauss points and are just summarized over integration points to one single vector for each cell. It means there should be nothing stored on the nodes yet, but we have to take the information to the global structure and node/vertices level.

For Q1 elements you have a nodal basis. In particular, assigning the value to the DoF with the corresponding support point and component is all you have to do. The Gauss quadrature points never include the vertices [1] so saying that the values represent forces for each vertex and at the same time that the values are computed within Gauss points sounds contradictory.
You can find out about the (global) DoF corresponding to the (vertex_no)th vertex via
cell->vertex_dof_index(vertex_no, i)
where i loops through the number of components, i.e i=0 for the x-component and i=1 for the y-component. The position of the (vertex_no)th can be obtained by:
cell->vertex(vertex_no)
Note that this approach is only suitable for Q1 elements. For higher polynomial degree, you have to calculate the values for all the support points of the finite element.
 

Best,
Daniel

[1] https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature

[Jean-Paul Pelteret]

Dear Seyed,

 

My problem is not the theory. 

Providing the theory is for our benefit (and therefore yours). No one here is necessarily an expert in what you're trying to accomplish. So, for all of your explanation, an equation or two might go a long way to help us understand exactly what you're trying to achieve and, therefore, help us suggest how you would accomplish it. If you can't express your problem clearly then please don't be surprised if you don't receive any meaningful assistance. Help us help you.

It is quite difficult to store data in deal.II since there is no processinfo or such a thing which tracks all data created and stored over the computation. ...

Deal.II does not attempt to be the beginning and end of provided functionality. Yes, there are times when deal.II does not provide some functionality that you require. There are a zoo of other libraries that you could also use to extend your code and fill in the gaps that deal.II does not provide. Of course, we always welcome extensions to existing features (e.g the relatively new CellDataStorage class and TransferableQuadraturePointData class) or the addition of new ones.

Now I try to compute the material forces for the whole structure. Hence, I tried using the same procedure of assembly for our LHS or RHS matrices in deal.II, but it fails.

Is this a deal.II error or a user error? Have you repopulate the local_dof_indices vector for each new cell that you're doing this "assembly" on? It looks like all contributions are going to a fixed set of global DoFs.

Regards,

Jean-Paul

[Seyed Ali Mohseni]

@Daniel:

So you have a FESystem<2>(FE_Q<2>(1),2) element, i.e. linear continuous elements with two components, and the vector represents corresponding local DoF values?

 Exactly. I defined first 

FESystem<dim> fe;

and then 

template<int dim>

SolidMechanics<dim>::SolidMechanics(): ..., fe(FE_Q<dim>(1), dim), ...

as constructor initializer. 

Here dim=2 would result in what you just mentioned.

In particular, assigning the value to the DoF with the corresponding support point and component is all you have to do.

Is there a command that allows me to assign values to the DoF. Can you please explain this a bit more.

The Gauss quadrature points never include the vertices [1] so saying that the values represent forces for each vertex and at the same time that the values are computed within Gauss points sounds contradictory.

Yes, you are right. I described it a bit mixed up. The vector is just shaped and designed for our global DoFs structure. You initialize it by:

 

Vector<double> cell_cf(dofs_per_cell);


Of course integration points and vertices are different. Thanks for the reference, it may come useful.

Note that this approach is only suitable for Q1 elements. For higher polynomial degree, you have to calculate the values for all the support points of the finite element.


Would you suggest another approach? Since it may be useful to implement this in a general way, useable for higher degrees of freedom. 

@Jean-Paul:

Providing the theory is for our benefit (and therefore yours). No one here is necessarily an expert in what you're trying to accomplish. So, for all of your explanation, an equation or two might go a long way to help us understand exactly what you're trying to achieve and, therefore, help us suggest how you would accomplish it. If you can't express your problem clearly then please don't be surprised if you don't receive any meaningful assistance. Help us help you.

Of course, but I just wanted to save you from additional headaches :) 

The computation of configurational forces is accomplished by the following formula:

Deal.II does not attempt to be the beginning and end of provided functionality. Yes, there are times when deal.II does not provide some functionality that you require. There are a zoo of other libraries that you could also use to extend your code and fill in the gaps that deal.II does not provide. Of course, we always welcome extensions to existing features (e.g the relatively new CellDataStorage class and TransferableQuadraturePointData class) or the addition of new ones.

How can I use these new features? Is there a step example or tutorial? Is this TransferableQuadraturePointData the same like in step-18?

 Is this a deal.II error or a user error? Have you repopulate the local_dof_indices vector for each new cell that you're doing this "assembly" on? It looks like all contributions are going to a fixed set of global DoFs.

I think more of a user error. Deal.II is nice and very elegantly designed. I just don't think it will automatically assemble it correctly in this case. First I have to understand it, before telling the program. ;)

I just did this:

std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();

for (; cell != endc; ++cell)

{

    if ( cell->is_locally_owned() )

    {
       fe_values.reinit(cell);

     ... configurational forces computation ...

    }
   
    configurational_forces = cell_cf;
   
    cell->get_dof_indices(local_dof_indices);

    constraints.distribute_local_to_global(cell_cf, local_dof_indices, configurational_forces);
}

I am really not sure, if the hanging node constraints have anything to do with my desired assembly at the moment. I assume they play a role when I refine the mesh.

So what exactly do you mean by repopulating each cell?

Many thanks, this group here is full of nice people. Hopefully, I can help you also as a reward in future. :)

[Wolfgang Bangerth]

Seyed,

But it doesn't work when we don't know what you want to do :-)

> The computation of configurational forces is accomplished by the
> following formula:
>
>

> <https://lh3.googleusercontent.com/-B38RUw135-c/WIjDUOvqsfI/AAAAAAAAAEI/4rD24zbgMeoI9O-2MFpXt77UTbGAHjmLQCLcB/s1600/config_forces.png>

>
>
> Deal.II does not attempt to be the beginning and end of provided
> functionality. Yes, there are times when deal.II does not provide
> some functionality that you require. There are a zoo of other
> libraries that you could also use to extend your code and fill in
> the gaps that deal.II does not provide. Of course, we always welcome
> extensions to existing features (e.g the relatively
> new CellDataStorage class

> <https://www.dealii.org/developer/doxygen/deal.II/classCellDataStorage.html> and TransferableQuadraturePointData
> class
> <https://www.dealii.org/developer/doxygen/deal.II/classTransferableQuadraturePointData.html>)

> or the addition of new ones.
>
>
> How can I use these new features? Is there a step example or tutorial?
> Is this TransferableQuadraturePointData the same like in step-18?
>
> Is this a deal.II error or a user error? Have you repopulate
> the local_dof_indices vector for each new cell that you're doing
> this "assembly" on? It looks like all contributions are going to a
> fixed set of global DoFs.
>
>
> I think more of a user error. Deal.II is nice and very elegantly
> designed. I just don't think it will automatically assemble it correctly
> in this case. First I have to understand it, before telling the program. ;)

Having read through this thread, I think that you are stuck because you
are familiar with one code and the corresponding notation, and you
expect that deal.II uses the same notation and does the same thing. So
you try to do in deal.II what you do there, but you can't find the same
terminology, the same functions, etc.

As an example:

> I just did this:
>
> |
> std::vector<types::global_dof_index>local_dof_indices(dofs_per_cell);
>
> typenameDoFHandler<dim>::active_cell_iterator cell
> =dof_handler.begin_active(),endc =dof_handler.end();
> for(;cell !=endc;++cell)
> {
> if(cell->is_locally_owned())
> {
> fe_values.reinit(cell);
>
> ...configurational forces computation ...
>
> }
>
> configurational_forces =cell_cf;

We don't know what "configurational forces" are. Or what, in your
formula, the various terms are.

I think that the only way to make progress is if you write down the
mathematics you want to implement, and then compare it with what the
various tutorial programs do. For example, for integrating the local
right hand side you may want to see steps 3 and 4, as well as 8. For how
exactly this interacts with hanging nodes, you may want to see step-6.

The point is that you need to expect to see similar mathematics, but
different terms than you are used to.

[Jean-Paul Pelteret]

Dear Seyed,

On top of what Wolfgang has already said, I have a few more direct comments about what you posted a few days ago.

 

The computation of configurational forces is accomplished by the following formula:

If one reinterprets the Eshelby (energy momentum) stress tensor as any other "regular" stress tensor, then the computation of $\mathbf{g}$ is exactly the same as the assembly of the RHS vector in nonlinear elasticity (without body forces and tractions). You can compare your implementation to that of step-42 and step-44.

 

 

How can I use these new features? Is there a step example or tutorial? Is this TransferableQuadraturePointData the same like in step-18?

CellStorageData is now used in step-44 (development version); for TransferableQuadraturePointData and ContinuousQuadratureDataTransfer you'll have to look at the unit tests; see this post for details on how to do that.

 

I just did this:

std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();

for (; cell != endc; ++cell)

{

    if ( cell->is_locally_owned() )

    {
       fe_values.reinit(cell);

     ... configurational forces computation ...

    }
   
    configurational_forces = cell_cf;
   
    cell->get_dof_indices(local_dof_indices);

    constraints.distribute_local_to_global(cell_cf, local_dof_indices, configurational_forces);
}

There are a few potential points that you can check here, and a few to be concerned about:

- Is the Eshelby stress tensor non-trivial on all (locally owned) cells?

- Is  cell_cf zero on cells that aren't locally owned? Why do you still assemble on non-locally-owned cells?

- Is cell_cf non-trivial on cells that are locally owned?

- Why do you set "configurational_forces = cell_cf;" ? Isn't the one a global vector and the other one with local contributions? This is suggested by the call to distribute_local_to_global. Does this even run in debug mode?

So what exactly do you mean by repopulating each cell?

I queried whether you called cell->get_dof_indices(local_dof_indices); for each cell before distribution. 

I agree with Wolfgang that its probably most helpful to you if you become more familiar with the "deal.II way" by studying the tutorials (even some not directly related to your topic of interest) more thoroughly.

Best regards,

Jean-Paul

[Seyed Ali Mohseni]

First of all, thank you for all your suggestions. They helped me a lot, especially from Prof. Bangerth and Jean-Paul. Many thanks ;)

I solved the assembly by listening to Prof. Bangerth's advice that I should try to implement things in a deal.II way. And after Jean-Paul told me to take a look at step-42 I realized how things are working in deal.II.

Unfortunately, I have still some issues with the values. They are not the same. I believe this is dependent on the B-operator, namely the jacobians. They differ by factor 2 probably since our unit cell in deal.II is within [0,1]^dim and the C++ code uses [-1,1]^dim.

But shouldn't the assembled configurational forces still be the same independent from the mapping?

Here is what I got so far :)

Kind regards,

Seyed Ali

[Wolfgang Bangerth]

On 01/28/2017 04:56 AM, 'Seyed Ali Mohseni' via deal.II User Group wrote:
>
> I somehow figured the assembly out by listening to Prof. Bangerth's advice
> that I should try to solve it the deal.II way. And after Jean-Paul told me to

> take a look at step-42 I realized how things are working in deal.II.
> Unfortunately, I have still some issues with the values. They are not the
> same. I believe this is dependent on the B-operator, namely the jacobians.
> They differ by factor 2 probably since our unit cell in deal.II is within
> [0,1]^dim and the C++ code uses [-1,1]^dim.
> But shouldn't the assembled configurational forces still be the same
> independent from the mapping?
>
> Here is what I got so far :)

I really have no idea what the pictures show, nor which ones correspond to the
deal.II solution and your own code, respectively.

But I do see that one seems to be computed on a finer mesh -- maybe a
triangular mesh, whereas the other one is computed on a quadrilateral mesh.

[Seyed Ali Mohseni]

@Prof. Bangerth: The upper one where pseudocolor is written, is from VisIt and deal.II. Below as a reference I showed my results from the C++ code. Both simulations are done using Q4 quadrilateral elements. There are no triangular elements here and the meshes are identical. They are intentionally created for comparison.

Here I show a more refined version, maybe you see the solution more clearly:

Additionally, I colored the configurational force arrows from the deal.II solution to black, so it becomes more visible.

I think I found the last small detail which seemed to be wrong, but am still not 100 % sure. There is something I don't understand completely. Surprisingly, the JxW value in deal.II stays the same value like in my C++ code, although the jacobian matrix differs by a factor of 2. Exactly this factor of 2 was the difference in my previous results. Now as you can see I get identical results in both approaches. 

I also fixed some other problem, namely the storage of my data for each element separately. Earlier I stored everything in Gauss points only, but the global vector for the PointHistory was overwritten since I had now several elements. For a single element it worked like a charm which is why I didn't recognize it directly. Now I am curious, if this is the best way to store data in a vector for large number of elements. I think you have probably some deal.II way of storing data on cells or at least efficiently in a vector using pointer arithmetics. Is the CellStorageData Jean-Paul mentioned what I seek? I really don't want to lose the scalability of my approach since it should be applicable to large problems also.

[Seyed Ali Mohseni]

The last piece of the puzzle is solved. The small problem I mentioned was due to the problem discussed here: https://groups.google.com/d/msg/dealii/5YXHbjPs3ls/bsGMAPwXBQAJ

Thank you everyone.

The material force computation is now working 100% correctly for any case.