Update algorithm of shear displacement

Hello Jibril,

thanks for checking the documentation and code so much in detail. After
checking the code, the quick answer is that you are right and documentation
and actual implementation do differ.

Another discussion is which one is correct and how much it influence the
result or if the difference is negligible (Bruno?).

Also interesting note, the parts of code and documentation you mentioned
are unchanged since 2010 :-)

cheers
Jan

2016-05-04 11:13 GMT+02:00 Jibril Coulibaly <
<email address hidden>>:

> New question #293235 on Yade:
> https://answers.launchpad.net/yade/+question/293235
>
> Hi everyone,
>
> My questions concern the algorithm used for the shear dispalcement update
> and its implementation is the class ScGeom. I would like to fully
> understand the update procedure but it seems like there are discrepancies
> between the documentation and the code.
>
> In the documentation (
> https://yade-dem.org/doc/formulation.html#fig-shear-2d ), the update is
> said to be "done by perpendicular projection to the plane first (which
> might decrease |\uT|) and adding what corresponds to spatial rotation of
> the interaction instead". The formulae for the corresponding increments
> don't indicate a projection and look more like infinitesimal rotations.
> These 2 rotations, acting of the previous shear dispalcement vector, are
> eventually linearly added, not composed, to form the current shear
> dispalcement vector.
>
> Getting into the corresponding code from the class ScGeom (
> https://github.com/yade/trunk/blob/master/pkg/dem/ScGeom.cpp ), the
> ScGeom::rotate function clearly performs the composition of the two
> approximate rotations, not the addition, about what is defined as the
> orthonormal axis first (line 17) and the twist axis second (line 18).
> The determination of these axes is performed in the ScGeom::precompute
> function. While the definition of the orthonormal axis (line 27)
> corresponds to what the documentation says, the definition of the twist
> angle and axis (lines 28,29) differs from the formula in the documentation
> in that is uses the variable "normal" for the current normal vector n° and
> not "currentNormal".
>
> Regarding the documentation, shouldn't the second rotation of the shear
> dispalcement (around the twist axis) be applied to the shear displacement
> already updated by the first rotation (around the orthonormal axis), as it
> is in the function ScGeom::rotate ? That is :
>
> (\Delta \uT)_2&=-(\prevuT + (\Delta \uT)_1) \times\left(\frac{\Delta t}{2}
> \currn \cdot (\pprev{\vec{\omega}}_1+\pprev{\vec{\omega}}_2)\right) \currn
>
> instead of
>
> (\Delta \uT)_2&=-\prevuT\times\left(\frac{\Delta t}{2} \currn \cdot
> (\pprev{\vec{\omega}}_1+\pprev{\vec{\omega}}_2)\right) \currn
>
> And keeping the summation of the 3 increments as it is.
>
> Regarding the ScGeom::precompute function, shouldn't the twist axis be
> around the current normal vector n° as it is in the documentation ? That is
> replacing the variable "normal" by "currentNormal" in lines (28,29) ? This
> way the rotations would be properly composed :
> - rotation around the orthonormal axis first brings the previous normal
> vector together with the current normal vector.
> - rotation around the twist axis (that is colinear to the current normal)
> then performs the remaining rotation.
>
> I would really appreciate your help and hope I have been understandable
> and clear.
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Jibril,

thanks a lot for the checking and analyzing the code and documentation.
As Jan said, we need to update outdated documentation. If you have
a wish to do it, feel free to join us. We are always opened for new
collaborators.

Regarding the correctness of ScGeom, we need to discuss it in more
details.

Best regards

Anton

2016-05-04 11:13 GMT+02:00 Jibril Coulibaly
<email address hidden>:
> New question #293235 on Yade:
> https://answers.launchpad.net/yade/+question/293235
>
> Hi everyone,
>
> My questions concern the algorithm used for the shear dispalcement update and its implementation is the class ScGeom. I would like to fully understand the update procedure but it seems like there are discrepancies between the documentation and the code.
>
> In the documentation ( https://yade-dem.org/doc/formulation.html#fig-shear-2d ), the update is said to be "done by perpendicular projection to the plane first (which might decrease |\uT|) and adding what corresponds to spatial rotation of the interaction instead". The formulae for the corresponding increments don't indicate a projection and look more like infinitesimal rotations. These 2 rotations, acting of the previous shear dispalcement vector, are eventually linearly added, not composed, to form the current shear dispalcement vector.
>
> Getting into the corresponding code from the class ScGeom ( https://github.com/yade/trunk/blob/master/pkg/dem/ScGeom.cpp ), the ScGeom::rotate function clearly performs the composition of the two approximate rotations, not the addition, about what is defined as the orthonormal axis first (line 17) and the twist axis second (line 18).
> The determination of these axes is performed in the ScGeom::precompute function. While the definition of the orthonormal axis (line 27) corresponds to what the documentation says, the definition of the twist angle and axis (lines 28,29) differs from the formula in the documentation in that is uses the variable "normal" for the current normal vector n° and not "currentNormal".
>
> Regarding the documentation, shouldn't the second rotation of the shear dispalcement (around the twist axis) be applied to the shear displacement already updated by the first rotation (around the orthonormal axis), as it is in the function ScGeom::rotate ? That is :
>
> (\Delta \uT)_2&=-(\prevuT + (\Delta \uT)_1) \times\left(\frac{\Delta t}{2} \currn \cdot (\pprev{\vec{\omega}}_1+\pprev{\vec{\omega}}_2)\right) \currn
>
> instead of
>
> (\Delta \uT)_2&=-\prevuT\times\left(\frac{\Delta t}{2} \currn \cdot (\pprev{\vec{\omega}}_1+\pprev{\vec{\omega}}_2)\right) \currn
>
> And keeping the summation of the 3 increments as it is.
>
> Regarding the ScGeom::precompute function, shouldn't the twist axis be around the current normal vector n° as it is in the documentation ? That is replacing the variable "normal" by "currentNormal" in lines (28,29) ? This way the rotations would be properly composed :
> - rotation around the orthonormal axis first brings the previous normal vector together with the current normal vector.
> - rotation around the twist axis (that is colinear to the current normal) then performs the remaining rotation.
>
> I would really appreciate your help and hope I have been understandable and clear.
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp