Questions on Step-18

154 views
Skip to first unread message

Michael Lee

unread,
Jun 22, 2021, 9:24:54 AM6/22/21
to deal.II User Group
Hello,
I have two questions when studying Step-18. 

1) Should there be a factor 1/2 when calculating the rotation matrix angle (code in tutorial: angle = std::atan(curl) ) ? This is a mathematical problem. Can anyone give some material for the formula to clear out my doubt?

2) In the introduction, it says "we will consider a model in which we produce a small deformation, deform the physical coordinates of the body by this deformation, and then consider the next loading step again as a linear problem. This isn't consistent, since the assumption of linearity implies that deformations are infinitesimal and so moving around the vertices of our mesh by a finite amount before solving the next linear problem is an inconsistent approach." My question is how to make it consistent. Can I just add the nonlinear part of the strain tensor like the following (of course in both get_strain functions)? I just want to consider the geometrical nonlinearity.

  template <int dim>
  inline SymmetricTensor<2, dim>
  get_strain(const std::vector<Tensor<1, dim>> &grad)
  {
    Assert(grad.size() == dim, ExcInternalError());

    SymmetricTensor<2, dim> strain;
    for (unsigned int i = 0; i < dim; ++i)
      strain[i][i] = grad[i][i];

    for (unsigned int i = 0; i < dim; ++i)
      for (unsigned int j = i + 1; j < dim; ++j)
        strain[i][j] = (grad[i][j] + grad[j][i] +       // linear part
                              grad[i][0] * grad[j][0] +    //nonlinear part
                              grad[i][1] * grad[j][1] +   
                              grad[i][2] * grad[j][2]) / 2;

    return strain;
  }

Thanks,
Michael

Jean-Paul Pelteret

unread,
Jun 24, 2021, 4:03:13 PM6/24/21
to dea...@googlegroups.com
Hi Michael,

I cannot comment on the first question, but might be able to assist a bit with the second. But may I first ask, what precisely are you trying to achieve with this extension? 

As interesting as it is, in the past I had found step-18 to deviate too significantly from the “classical” approach to elasticity to be a natural candidate extend towards finite strain elasticity, for example (this is kind-of implied by the caveat that you partially quoted). I’ve been looking at some of my textbooks (e.g. reviewing the topic of the updated Lagrangian formulation in Holzapfel’s “Nonlinear solid mechanics” and Wrigger’s “Nonlinear finite element methods”) to try to answer (2), but cannot trivially reconcile the two approaches. I think that there’s a bit too much going on to be able to correctly deduce by eye what the requisite changes are. I think that you’d need to reformulate the balance laws and consider their implications for the implemented weak forms — in particular, I think that you’d be missing a contribution that (for finite deformation) looks like the geometric stiffness, but there could be further differences that extend from the overall approach taken to the problem. 

If you’re interested in an examples that are more closely aligned with what you might see in the literature, then you can take a look at step-44 or the “Quasi-Static Finite-Strain Compressible Elasticity” code-gallery example, which is effectively step-44 reduced to the one-field total Lagrangian formulation that you’d find in many standard textbooks. It would be easy enough to rework this to use the updated Lagrangian approach, if that it what you desire.

I hope that this helps you a little.

Best,
Jean-Paul


--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/02f465df-3cab-4d84-bd48-92ca4a7981efn%40googlegroups.com.

Andrew McBride

unread,
Jun 24, 2021, 4:23:06 PM6/24/21
to deal.II User Group
Hi both,

Jean-Paul has addressed the second point nicely. On the first point, I think there is a 1/2 missing. The curl of the velocity gradient is the vorticity which is twice the angular velocity - hence I think you need a 1/2. Happy to be corrected on this.

Best,
Andrew

Michael Li

unread,
Jun 24, 2021, 8:32:58 PM6/24/21
to dea...@googlegroups.com

Andrew, thanks for confirming that. The missing 1/2 does not affect the demonstration of functionalities of deal.II but it may change the results.

 

Jean-Paul, thanks for commenting on my second question. I want to study the pure geometrical nonlinear elasticity (large deformation with linear material). It should be the simplest nonlinear model in elasticity. Step 18 looks like a good start as an updated Lagrangian formulation; but it does not include the nonlinear part of the finite strain. I spent some on Step 44 because “ Quasi-Static Finite-Strain Compressible Elasticity” also mentions it is based on Step 44. To be honest, I have a hard time to understand the three-filed formulation in Step 44. So maybe I should get into “Quasi-Static … Elasticity” first which uses classical one-field formulation. However, I’m still interested in implementing the pure geometrical large finite deformation model. I believe that will make things simpler and help me understand more complex nonlinear models. Meanwhile, I try to compare the results with a reference which only consider the geometrical nonlinearity.

 

-Michael

Jean-Paul Pelteret

unread,
Jun 25, 2021, 1:43:20 AM6/25/21
to dea...@googlegroups.com
HI Michael.

To be honest, I have a hard time to understand the three-filed formulation in Step 44. 

That’s a fair comment, and a known deficit in our current tutorial list. I hope that by the end of the year we will have a tutorial that uses one-field elasticity, and would be a much better starting point for you and others (WIP: https://github.com/dealii/dealii/pull/10394).

However, I’m still interested in implementing the pure geometrical large finite deformation model. I believe that will make things simpler and help me understand more complex nonlinear models. Meanwhile, I try to compare the results with a reference which only consider the geometrical nonlinearity.

It sounds to me like you’re looking to implement finite strain elasticity with a St Venant-Kirchhoff constitutive law. That is, as far as I understand, the simplest constitutive law for finite strain elasticity.
In this model, the elastic tangent looks like that for linear elasticity, but the strain tensor of interest is no longer the linear/small strain tensor, but rather the (nonlinear) Green-Lagrange strain tensor. (Nevertheless, even though the material law is linear, you still have the geometric tangent contribution that appears after linearisation of the residual, which comes from the linearisation of the now nonlinearly deformation-dependent test function “de” in  \int de : \sigma dv .)

If this is indeed what you want, then it is possible to swap out the Neo-Hookean constitutive law that is implemented in the code-gallery example with the St-Venant Kirchhoff constitutive law. You’d get “for free” the geometric stiffness contribution to the tangent, and would just have to worry about the material tangent and definition of the stress. I think that the only complication is the default parameterisation for the constitutive law — it was parameterised in terms of the left Cauchy-Green tensor as the whole problem is set up with the Kirchhoff stress (that’s the Cauchy stress, but per unit reference volume). Nevertheless, with the help of these “push-forward” transformations you would be able to use whatever parameterisation you’d like for the constitutive law (e.g. in terms of the Green-Lagrange strain tensor) and transform the resulting stress and material tangent to the correct (i.e. the spatial) configuration as is required for that implementation of the weak forms. The important thing to remember is that, despite all of the transformations involved, you’re still solving the same BVP with the same constitutive law irrespective of which choice of stress-strain (energetic) conjugate pairs you use, and irrespective of whether its expressed in the total or updated Lagrangian form.

But anyway, this is just a suggestion. Irrespective of whether or not you choose to go that route, I wish you luck in implementing what you’re needing to.

Best,
Jean-Paul

Michael Li

unread,
Jun 25, 2021, 11:36:05 AM6/25/21
to dea...@googlegroups.com

Hi Jean-Paul,

Thank you so much for the deliberate suggestions! It seems the code-gallery example is a good starting point for me to implement St-Venant Kirchhoff hyperelastic model with the help of “push-forward” transformations. I need to dig into this example to understand how the geometric stiffness tangent is embedded for free.

 

It’s good to know a tutorial that uses one-filed elasticity will be available. I believe that will help a lot for beginners. AD is an interesting tool for nonlinear FEM. In  an early discussion, Alex shared his similar codes using AD at Integrated material and spatial traction forces on boundary not equal (google.com). But since I’m a deal.II beginner and not familiar with AD, I have not studied his code further.

 

So I am going to study your code-gallery example first as it is well documented. I may ask for your further help and will let you know if I get the code roll.

 

Thank you again for your wonderful help!

 

Best,

Michael

Wolfgang Bangerth

unread,
Jun 25, 2021, 11:58:11 PM6/25/21
to dea...@googlegroups.com
On 6/24/21 6:32 PM, Michael Li wrote:
> Andrew, thanks for confirming that. The missing 1/2 does not affect the
> demonstration of functionalities of deal.II but it may change the results.

If this is wrong in the deal.II code, would you be interested in writing a
patch to correct this? It should be an easy fix, and a good first patch to
contribute to the library!

If you're curious about ho to write a patch, take a look at
https://github.com/dealii/dealii/wiki/Contributing

Best
W.


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

Jean-Paul Pelteret

unread,
Jun 26, 2021, 3:48:53 AM6/26/21
to dea...@googlegroups.com
Following Andrew’s explanation, I suppose that this is relation for which we’re lacking the factor of 1/2, right?


If so, then maybe we should link to this equation in the tutorial documentation too.

If this is wrong in the deal.II code, would you be interested in writing a patch to correct this? It should be an easy fix, and a good first patch to contribute to the library!

I concur — it would be great if you’d be willing to write a patch that fixes this! It would be interesting to see the effect on the results of the tutorial.

Best,
Jean-Paul


--
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+un...@googlegroups.com.

Michael Lee

unread,
Jun 26, 2021, 10:15:01 AM6/26/21
to dea...@googlegroups.com
I would be happy to do the comparison and make the fix. But before doing that, I want to make sure I understand the formula correctly.

When we calculate the rotation matrix, it seems that the du (incremental displacement) is used.
          // Then initialize the <code>FEValues</code> object on the present
          // cell, and extract the gradients of the displacement at the
          // quadrature points for later computation of the strains
          fe_values.reinit(cell);
          fe_values.get_function_gradients(incremental_displacement,
                                           displacement_increment_grads);

Should we use the velocity (\vec{v} = \frac{d \vec{u}}{dt}), since \vec{\omega} = \frac{1}{2}\nabla \times \vec{v} = \frac{1}{2}\nabla \times \frac{d\vec{u}}{dt} ?

Can you check this as well?


Best,
Michael




To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1D258AD6-FC4C-43B7-9056-A4B21114D694%40gmail.com.

Jean-Paul Pelteret

unread,
Jun 26, 2021, 4:57:38 PM6/26/21
to dea...@googlegroups.com
Hi Michael,

What’s written in the implementation suggests that computing the curl of the displacement (increments) is the right thing to do in this instance:

Nevertheless, if the material was prestressed in a certain direction, then this direction will be rotated along with the material. To this end, we have to define a rotation matrix R(Δun) that describes, in each point the rotation due to the displacement increments. It is not hard to see that the actual dependence of R on Δun can only be through the curl of the displacement, rather than the displacement itself or its full gradient (as mentioned above, the constant components of the increment describe translations, its divergence the dilational modes, and the curl the rotational modes). 

I think that this is because the displacement field for solids is analogous to the velocity field for fluids, for which that equation in the Wiki link was expressly written (Compare incompressible linear elasticity to Stokes flow — the equations have the same form). Furthermore, having dug through my continuum mechanics notes, I now see that this specific rotation matrix is called the “infinitesimal rotation tensor” and is defined as the antisymmetric part of \grad u. Since step-18 is dealing with incremental updates, the incremental form of the rotation tensor is computed.

Sorry for the confusion. So, to my understanding it’s just the factor of 1/2 that needs to be added.

Best,
Jean-Paul

Michael Li

unread,
Jun 26, 2021, 9:19:30 PM6/26/21
to dea...@googlegroups.com

Hi Jean-Paul,

Thank you for the nice explanation and information! It cleared out my concern and helped me understand the related concepts of vorticity and rotation. Vorticity (1/2 curl of velocity) stands for the rotation rate (angular velocity) but the infinitesimal rotation tensor gives the rotation angle (curl of displacement). Do correct me if I’m wrong.

 

Now I’m excited to learn how to make my first patch.

 

- Michael

 

From: Jean-Paul Pelteret
Sent: Saturday, June 26, 2021 2:57 PM
To: dea...@googlegroups.com
Subject: Re: [deal.II] Questions on Step-18

 

Hi Michael,

Andrew McBride

unread,
Jun 28, 2021, 4:21:36 AM6/28/21
to deal.II User Group
My view here is that the problem is quasi-static and hence time serves to order events. Hence a displacement increment can be viewed as equivalent to the velocity. 

For example, let’s fix the number of time-steps to 10. If the total time is 1s or 10s we should get the same results (assuming the timing of the application of loading is scaled based on the assumption that time simply order events).

Best,
Andrew

Michael Lee

unread,
Jul 2, 2021, 8:56:51 PM7/2/21
to dea...@googlegroups.com

Hi Jean-Paul and All,

 

I would like to discuss the formulas in step-18 further.

 

  • In the code, the angles are computed as follows.

For 2D,

    const double curl = (grad_u[1][0] - grad_u[0][1]);

    const double angle = std::atan(curl);

 

For 3D,

    const double tan_angle = std::sqrt(curl * curl);

    const double angle     = std::atan(tan_angle);

 

I don’t really understand why we need to do the atan operation to get the angle. It seems to me that half of the curl of displacement itself is the rotation angle. 


Moreover, for 2D, do we need to do the same as 3D   tan_angle = std::sqrt(curl * curl)? I mean we first calculate the magnitude of the curl then calculate the angle. I’m concerned about the sign of the curl in 2D angle = std::atan(curl) (we need the magnitude of the curl).

 

I wonder what references I can get from these formulas. I hope I can be fully convinced of the stuff.

 

  • To be specific, I want to make sure where I should put the factor of ½ to make the patch. I.e., which is correct in the following?

 

  1. const double tan_angle = std::sqrt(0.5 * curl * 0.5 * curl);

or


  1. const double tan_angle = 0.5 * std::sqrt(curl * curl);


  • In my opinion, the code should be as follows. 


For 2D,

    const double curl = (grad_u[1][0] - grad_u[0][1]);


    const double angle = 1.0 / 2.0 * std::sqrt(curl * curl); // Half of the magnitude of curl is the rotation angle.

 

For 3D,

    const double angle = 1.0 / 2.0 * std::sqrt(curl * curl);

    const double angle     = std::atan(tan_angle);

 


Any comments will be of great help.


Best,

Michael


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/B-jNRYdzqzs/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/1AB07F53-AD43-4065-BDA7-C03B1D8F0E57%40gmail.com.

Jean-Paul Pelteret

unread,
Jul 6, 2021, 1:50:46 AM7/6/21
to dea...@googlegroups.com
HI Michael,

I’m not the original author of step-18, but I think that I’ve found some sources that can explain both the construction of the non-normalised rotation axis w and the angle θ calculation.

1. Axial vector ω (the non-normalised axis of rotation)

For this, I’m summarising the salient parts of the references that I list later; specifically for following equations: 

Chadwick1998a: p 29 eq. 71, the paragraph on p79 between equations 46 and 47
Holzapfel2007a: p 17 eq. 1.118, p 18 eq. 1.124, p 49 eq. 1.276

The axial vector w is related to some skew symmetric tensor W in the following manner:
W.u = ω x u  for any vector u.

It can then be shown that
ω = -1/2[ ε_{ijk} W_{ij} e_{k}  ]
   = 1/2 curl(u)

So how this factor comes into the calculation would be most easily understood by introducing an intermediate quantity,

const Point<3> curl = …;
const Tensor<1,3> axial_vector = 0.5*curl;
const double tan_angle = std::sqrt(axial_vector * axial_vector);
const double angle = std::atan(tan_angle);

With some generalisation, this should hold in 2-d and 3-d (and addresses your one observation, that I comment on again in the next section).


@Book{Chadwick1998a,
  title     = {Continuum Mechanics: Concise Theory and Problems},
  publisher = {Dover Publications Inc.},
  year      = {1998},
  author    = {Chadwick, P.},
  address   = {Mineola, New York, USA},
  edition   = {2$^{\text{nd}}$},
  isbn      = {9780486401805},
  owner     = {Jean-Paul Pelteret},
  timestamp = {2016.01.20},
}

@Book{Holzapfel2007a,
  title     = {Nonlinear Solid Mechanics: A Continuum Approach for Engineering},
  publisher = {John Wiley \& Sons Ltd.},
  year      = {2007},
  author    = {Holzapfel, G. A.},
  address   = {West Sussex, England},
  isbn      = {0-471-82304-X},
  file      = {Holzapfel2007a.pdf:References/Books/Holzapfel2007a.pdf:PDF},
  owner     = {Jean-Paul Pelteret},
  timestamp = {2014.09.26},
}


2. Angle from ω

For the angle itself, I defer to this wikipedia entry. 
There you can see that the rotation angle θ is the arctangent of the norm of the non-normalised axis magnitude, |ω|.

As for why there is the inverse trigonometric operation involved, this can be understood most easily by thinking about the units involved. The norm |ω| is dimensionless, and we’re needing an angle in radians — an inverse trigonometric function does this conversion. As to why its the arctan and not something else… maybe some of the references on that page have a proof for the formula that would then provide that insight.

Moreover, for 2D, do we need to do the same as 3D   tan_angle = std::sqrt(curl * curl)? I mean we first calculate the magnitude of the curl then calculate the angle. I’m concerned about the sign of the curl in 2D angle = std::atan(curl) (we need the magnitude of the curl).

I think that your observation in 2-d is correct. According to the wikipedia entry, one should likely always use the magnitude of the curl. I also think that the where the rotation matrix is returned via the call to   
Rotations::rotation_matrix_<2,3>d(axis, -angle);   
, we should be passing in the angle and not its negation. The axial vector ω supposedly encodes both the axis direction and “sign" of the rotation angle, and should be treated as a pair. That reference says that one should construct the rotation matrix from either (ω, θ) or (-ω, -θ), which would ultimately lead to the same result.

Thanks for being inquisitive and for asking all of these questions. It looks like its lead to a better understanding (which we now use to improve the documentation of the tutorial), and might have uncovered a couple of bugs!

Best,
Jean-Paul


Michael Li

unread,
Jul 7, 2021, 11:43:02 PM7/7/21
to dea...@googlegroups.com

Hi Jean-Paul,

 

Thank you for your patience and detailed explanation! It did help me understand the formulations. The only thing after reading the Wikipedia (https://en.wikipedia.org/wiki/Rotation_matrix#Axis_and_angle ) I still have doubts is that the angle should be calculated by ½ mag(curl) and we do not need to do the atan operation. And I think the magnitude of curl has the dimension of radians which actually is dimensionless.

 

Quoted from Wikipedia:

“Near multiples of 180°, care is needed to avoid numerical problems: in extracting the angle, a two-argument arctangent with atan2(sin θ, cos θ) equal to θ avoids the insensitivity of arccos; and in computing the axis magnitude in order to force unit magnitude, a brute-force approach can lose accuracy through underflow .

”

 

My understanding on the atan above is that it is only used to deal with the case when the angle is near multiples of 180°. So the code I suggest will look like as follows.

       const Point<3> curl = …;

       const double mag_curl = std::sqrt(curl * curl);

       const double  angle = 0.5 * mag_curl;

Thank you!

Michael

 

 

From: Jean-Paul Pelteret
Sent: Monday, July 5, 2021 11:50 PM
To: dea...@googlegroups.com
Subject: Re: [deal.II] Questions on Step-18

 

HI Michael,

1.      const double tan_angle = std::sqrt(0.5 * curl * 0.5 * curl);

or

 

2.      const double tan_angle = 0.5 * std::sqrt(curl * curl);

Reply all
Reply to author
Forward
0 new messages