← Back to team overview

yade-dev team mailing list archive

Re: computation of MomentBending

 

Please go on. 

I am attaching a full example. The toAxisAngle routine implemented by
hand is only to check the calculation results. It is not how it
_should_ be implemented. In fact I didn't look at wm3 or at any other
source code , I just looked at the quaternion formulas and coded them down.

The output of attached program is following:

check with boost:
sub products:
(1,-1.56132e-11,7.47311e-11,-1.7066e-10) * (1,9.50188e-19,-2.45305e-19,0) - Bruno's
(1,7.62493e-11,-3.53675e-12,-8.53304e-11) * (1,9.50187e-19,0,0) - result
total:
(1, -1.56132e-11, 7.47311e-11, -1.7066e-10) - Bruno's
(1,7.62494e-11,-3.53677e-12,-8.53305e-11) - boost result

check with wm3:
sub products:
(1,-1.56132e-11,7.47311e-11,-1.7066e-10) * (1,9.50188e-19,-2.45305e-19,0) - Bruno's
1 7.62493e-11 -3.53675e-12 -8.53304e-11 * 1 9.50187e-19 -2.45305e-19 0 - result
total:
1, -1.56132e-11, 7.47311e-11, -1.7066e-10 - Bruno's
1 7.62494e-11 -3.53677e-12 -8.53305e-11 - wm3 result

Bruno's: (1 0 0), 0   - with wm3
by hand: (1 0 0), 0
wm3    : (1 0 0), 0
            

summarising: it is possible to obtain nan's, if toAxisAngle is not
taking care of corner cases, and if quaternion diverges too much
from identity. In your example numbers, the nan was coming from
arcus_cos(1.00001), which of course must be nan, since acos doesn't
accept arguments higher than 1.0           

best regards
Janek Kozicki


Bruno Chareyre said:     (by the date of Tue, 26 Oct 2010 17:09:12 +0200)

> Below is an example (from revno 2514). If somebody see an obvious reason 
> for getting a NaN, please let me know; if not, I'll send this to Eigen 
> (the first line should be enough).
> 
> Bruno
> 
> 
> NaN angle found in angleAxisr(q), for quaternion 1 -1.56132e-11 
> 7.47311e-11 -1.7066e-10, after quaternion product
> rbp1.ori * (initialOrientation1.conjugate())) * (initialOrientation2 * 
> (rbp2.ori.conjugate()) is:
> 0.745315 -0.0308916 -0.665997 -1.14489e-10 * 0.745315 0.0308916 0.665997 
> -0 * 0.746671 -0.166278 -0.644077 0 * 0.746671 0.166278 0.644077 
> -1.47527e-18
> with sub-products :
> 1 -1.56132e-11 7.47311e-11 -1.7066e-10 * 1 9.50188e-19 -2.45305e-19 0
> 
> 
> On 26/10/10 15:17, Václav Šmilauer wrote:
> >> If I'm the author of this comment in the sources (I am, probably),
> >> "close to identity" means close to null rotation. In other words the
> >> corresponding rotation matrix is very close to identity matrix.
> >> So, let M be a rotation matrix defining a rotation of 10^-8 degrees.
> >> Theoretically, M*M^(-1)=Id.
> >> Now, do the same numerically with quaternions : let q be the quaternion
> >> representing the same rotation as matrix M, and compute q*q.conjugate().
> >> It seems, sometimes in Eigen, the result will be NaN.
> >> As a workaround, I set angle=0 when angle=NaN.
> >>      
> > Bruno, could you post that issue to eigen's forums
> > http://eigen.tuxfamily.org/index.php?title=Main_Page#Get_support? They
> > are very responsive as far as my experience goes and will give us proper
> > mathematical explanation about the cause.
> >
> > v
> >
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~yade-dev
> > Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
> > Unsubscribe : https://launchpad.net/~yade-dev
> > More help   : https://help.launchpad.net/ListHelp
> >
> >    
> 
> 
> -- 
> _______________
> Bruno Chareyre
> Associate Professor
> ENSE³ - Grenoble INP
> Lab. 3SR
> BP 53 - 38041, Grenoble cedex 9 - France
> Tél : +33 4 56 52 86 21
> Fax : +33 4 76 82 70 43
> ________________
> 
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev
> More help   : https://help.launchpad.net/ListHelp
> 


-- 
Janek Kozicki                               http://janek.kozicki.pl/  |

Attachment: 32-quaternions-example.tar.bz2
Description: application/bzip


Follow ups

References