RE: [deal.II] Questions on Step-18

Wed, 07 Jul 2021 20:43:05 -0700

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: dealii@googlegroups.com
Subject: Re: [deal.II] Questions on Step-18

 

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

 



On 3. Jul 2021, at 02:56, Michael Lee <lianxi...@gmail.com> wrote:

 

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

 

2.      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

 

On Mon, Jun 28, 2021 at 2:21 AM Andrew McBride <mcbride.and...@gmail.com> wrote:

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



On 27 Jun 2021, at 02:19, Michael Li <lianxi...@gmail.com> wrote:

 

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: dealii@googlegroups.com
Subject: Re: [deal.II] Questions on Step-18

 

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

 

On 26. Jun 2021, at 16:14, Michael Lee <lianxi...@gmail.com> wrote:

 

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

 

 

 

On Sat, Jun 26, 2021 at 1:48 AM Jean-Paul Pelteret <jppelte...@gmail.com> wrote:

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

 

 

On 26. Jun 2021, at 05:58, Wolfgang Bangerth <bange...@colostate.edu> wrote:

 

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.

 

 

On Sat, Jun 26, 2021 at 1:48 AM Jean-Paul Pelteret <jppelte...@gmail.com> wrote:

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

 

 

On 26. Jun 2021, at 05:58, Wolfgang Bangerth <bange...@colostate.edu> wrote:

 

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:                 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 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/32840123-c6f5-abd9-18bb-03f06b5f918e%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/1D258AD6-FC4C-43B7-9056-A4B21114D694%40gmail.com.

 

-- 
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/CAEyr2zhNeyCO_tyXpw%2Bzgh4DAwzCmAGHG5W9y8RnLMswq_O10A%40mail.gmail.com.

 

-- 
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/1ADF3B07-54FB-4DDA-A848-1D7F9168D4A6%40gmail.com.

 

 

-- 
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/99FD763E-ABDF-46F5-A62C-D4427C562103%40hxcore.ol.

 

 

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

 

-- 
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/CAEyr2zjnynMo7Ar2eJ-MB%2Bo8PQVkdB270abKHZ7Bq-fbWWMRFg%40mail.gmail.com.

 

--
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/BFBD4892-0C15-4692-8510-E0009522D4D0%40gmail.com.

 

--
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/799BC111-6DF5-4DDF-9198-B48A3657CA1F%40hxcore.ol.

Reply via email to