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.
https://en.wikipedia.org/wiki/Rotation_matrix#Axis_and_angle
<https://en.wikipedia.org/wiki/Rotation_matrix#Axis_and_angle>
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?
>
>
> const double tan_angle = std::sqrt(0.5 * curl * 0.5 * curl);
> or
>
>
>
> 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
> <mailto: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
>> <mailto: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 <mailto:jppelte...@gmail.com>
>> Sent: Saturday, June 26, 2021 2:57 PM
>> To: dealii@googlegroups.com <mailto: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
>> <https://antoniobilotta-structuralengineer.github.io/MyWebSite/SMbook_en/infinitesimal_sec_strain_chap_en.html>”
>> 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
>> <mailto: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 (), since ?
>>
>> Can you check this as well?
>>
>>
>> Best,
>> Michael
>>
>>
>>
>>
>>
>>
>> On Sat, Jun 26, 2021 at 1:48 AM Jean-Paul Pelteret <jppelte...@gmail.com
>> <mailto: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?
>>
>> https://en.wikipedia.org/wiki/Angular_velocity#Angular_velocity_as_a_vector_field
>>
>> <https://en.wikipedia.org/wiki/Angular_velocity#Angular_velocity_as_a_vector_field>
>>
>> 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
>> <mailto: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
>> <https://github.com/dealii/dealii/wiki/Contributing>
>>
>> Best
>> W.
>>
>>
>>
>> On Sat, Jun 26, 2021 at 1:48 AM Jean-Paul Pelteret <jppelte...@gmail.com
>> <mailto: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?
>>
>> https://en.wikipedia.org/wiki/Angular_velocity#Angular_velocity_as_a_vector_field
>>
>> <https://en.wikipedia.org/wiki/Angular_velocity#Angular_velocity_as_a_vector_field>
>>
>> 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
>> <mailto: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
>> <https://github.com/dealii/dealii/wiki/Contributing>
>>
>> Best
>> W.
>>
>>
>> --
>> ------------------------------------------------------------------------
>> Wolfgang Bangerth email: bange...@colostate.edu
>> <mailto:bange...@colostate.edu>
>> www: http://www.math.colostate.edu/~bangerth/
>> <http://www.math.colostate.edu/~bangerth/>
>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> <http://www.dealii.org/>
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> <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
>> <mailto: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
>>
>> <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/
>> <http://www.dealii.org/>
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> <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
>> <mailto: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
>>
>> <https://groups.google.com/d/msgid/dealii/1D258AD6-FC4C-43B7-9056-A4B21114D694%40gmail.com?utm_medium=email&utm_source=footer>.
>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> <http://www.dealii.org/>
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> <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
>> <mailto: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
>>
>> <https://groups.google.com/d/msgid/dealii/CAEyr2zhNeyCO_tyXpw%2Bzgh4DAwzCmAGHG5W9y8RnLMswq_O10A%40mail.gmail.com?utm_medium=email&utm_source=footer>.
>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> <http://www.dealii.org/>
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> <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
>> <mailto: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
>>
>> <https://groups.google.com/d/msgid/dealii/1ADF3B07-54FB-4DDA-A848-1D7F9168D4A6%40gmail.com?utm_medium=email&utm_source=footer>.
>>
>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> <http://www.dealii.org/>
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> <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
>> <mailto: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
>>
>> <https://groups.google.com/d/msgid/dealii/99FD763E-ABDF-46F5-A62C-D4427C562103%40hxcore.ol?utm_medium=email&utm_source=footer>.
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> <http://www.dealii.org/>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <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
> <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
> <mailto: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
>
> <https://groups.google.com/d/msgid/dealii/1AB07F53-AD43-4065-BDA7-C03B1D8F0E57%40gmail.com?utm_medium=email&utm_source=footer>.
>
> --
> The deal.II project is located at http://www.dealii.org/
> <http://www.dealii.org/>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <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
> <mailto: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
>
> <https://groups.google.com/d/msgid/dealii/CAEyr2zjnynMo7Ar2eJ-MB%2Bo8PQVkdB270abKHZ7Bq-fbWWMRFg%40mail.gmail.com?utm_medium=email&utm_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/BFBD4892-0C15-4692-8510-E0009522D4D0%40gmail.com.
