yade-dev team mailing list archive

[Anton Gladky <gladky.anton@xxxxxxxxx>]

Hi, Janek.

I dont see any weird particle/clumps behavior now.

If you, Janek, have an idea how to improve/stabilize YADE results, you are
welcome to create separate branch, like it has been done here
https://code.launchpad.net/~yade-dev/yade/wm3-eigen for testing purposes.
You can test there your ideas, and when we are sure, it is ok (like it was
with eigen-migration process), we can simply merge branches.

Please, do not commit to upstream controversial thoughts, it is very sadly
to see early in the morning, that something is broken.

Thank you for understanding and have a luck in a future YADE development.
____________________________

Anton Gladky


2010/6/4 Janek Kozicki <janek_listy@xxxxx>

> Hi,
>
> you shouldn't be running wild through the code adding
> quaternion.normalize() every second line in the code.
>
> I've went through this stage already, and let me explain what is the
> problem and how to proceed.
>
>
> Quaternions:
>
> Quaternions, are just like complex numbers, except that instead of
> single imaginary "i", they have three imaginary numbers: "i", "j", "k".
>
> We are using them to perform rotations in 3D space. Namely, we are
> using them as a representation of an SO(3) rotation group:
> http://en.wikipedia.org/wiki/SO(3)<http://en.wikipedia.org/wiki/SO%283%29>.
> Fortunately multiplying the
> quaternions is just like composition of two rotations.
> (BTW: orthogonal rotation matrices 3x3 are also a representation of SO(3)
> ).
>
> Quaternions, are quite like complex numbers, and can be of ANY
> magnitude (length) that you could imagine. You can calculate sin(q),
> ln(q), e^q, solve polynomial equations with them, and other stuff.
> Which we are not interested in.
>
>
> The problem:
>
> The problem is that any given quaternion belongs to the SE(3)
> rotation group only, and only, if its length equals ONE.
> If any quaternion suddenly doesn't have length=1, it does not belong
> any more to SE(3).
>
> Analytically multiplying quaternions with length=1 will always
> produce a quaternion with length=1.
>
> Numerically multiplying two quaternions requires 12 additions and 16
> multiplications. There is a numerical error of order 1e-17 which can
> build up, making the resultant quaternion<double> to be not unitary
> any more.
>
>
> The way to solve this:
>
> Analytically + theoretically if a quaternion's length is not 1, then
> this is for us a critical error since it is NOT a rotation any more.
>
> Numerically + practically there is a numerical error introduced by
> the computer. It is of magnitude std::numeric_limit<TYPE>::epsilon()
>
> http://www.cplusplus.com/reference/std/limits/numeric_limits/
> http://msdn.microsoft.com/en-us/library/6x7575x3%28VS.71%29.aspx
>
> (where used TYPE is usually float, double or long double, for double
> IIRC it was 1e-17).
>
> We have to cope with this error, and due to this error it is
> impossible to have a quaternion's length always equal to 1.
>
>
> How this epsilon() propagates through calculations (12 additions, 16
> multiplications) depends exactly on the way in which it was
> implemented. The worst case would be: a single epsilon() becomes
> 12*16*epsilon() larger after single multiplication. I am sure that
> this is never the case. But you get the picture. We can call that
> 12*16 a MAX_ERROR and use 200 instead of 12*16=192.
>
>
> The problem again:
>
> So, we want to use UNITARY quaternions, but the library (wm3 and
> eigen) provide quaternions ONLY. Not the unitary quaternions.
>
> Our problem is inherent to the mathematics of SO(3) and computer's
> numerical errors. The end user must never worry about it.
>
>
> The solution:
>
> We need a library that provides UNITARY quaternions.
>
> If we don't have such library, then we need to normalize quaternions
> every second line in the code. And suddenly YADE's code base is
> bigger by 5000 lines.
>
> The only correct way to cope with this, is to extend existing
> quaternion library with a simple function that checks if quaternion
> is UNITARY.
>
> if quaternion's length differs from 1, then:
>
> - if it differs less that MAX_ERROR*epsilon(). Then the error is not
>  critical. We have to normalize that quaternion (divide it by its
>  length) and silently move on.
>
> - if the difference from 1 is bigger than MAX_ERROR*epsilon()
>  then YADE has to report a critical error, throw exception and instantly
>  terminate all calculations. Since whatever calculations there were,
>  they ARE wrong.
>
> This simple function has to be called at the end of each non-const
> member function in quaternion class.
>
>
>
> Let me repeat the solution:
>
> 1) Do NOT run wild through the code normalizing quaternions everywhere.
>
> 2) Do check if quaternion is unitary at the end of each non-const
> member function of quaternion class.
>
> 3) quaternions are in eigen lib and we can not modify it? Then we
> have to either:
>   - send them a patch, witch enables the library to provide unitary
>     quaternions, or
>   - derive from their class and extend it.
>
>
> I hope that it is clear now.
>
> This error is NOT of type: OMG all our calculations were wrong!
>
>
> The error build-up exists, usually is negligible. It shows itself only if:
>
> - same quaternion is used a lot in a single loop. lading to lots
>  of multiplications in single iteration. This is not the case in
>  DEM. Such multiplication in DEM with spheres is done just once per
>  iteration per quaternion.
>
> - a quaternion is mistakenly initialized to be non-unitary. Which CAN
>  be the case quite often, if you don't read what you are writing.
>
>
> I stumbled on this error with lattice model, because rotating a
> single rod required about 10 to 20 quaternion multiplications per
> iteration per quaternion. And I thought that it should be fixed in miniwm3,
> but:
> - I was too lazy to modify miniwm3
> - I was scared of computation overhead that this is going to introduce
> - I saw that in fact this problem manifests itself rarely.
>
>
> In fact that function instead of doing unitarity check on every call,
> could count the number of times that it was called, and perform
> normalization only every 10th call. That solution doesn't seem to be
> "perfect", but I tend to like it.
>
> Maybe it is better to do all those checks ALWAYS (each time) in DEBUG mode.
>
> And do them every 10th time in non debug mode. You cannot skip them
> totally in non-debug mode, because in non-debug mode the numerical
> buildup error occurs too!
>
>
> well, I hope that now you understand what's going on here, and aren't
> scared anymore
> --
> Janek Kozicki                               http://janek.kozicki.pl/
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev<https://launchpad.net/%7Eyade-dev>
> Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev<https://launchpad.net/%7Eyade-dev>
> More help   : https://help.launchpad.net/ListHelp
>