[SOLVED] bouncing on joints = EXPLOSION

[Julio Jerez]

Ok before I go on, I recreated the best I could the demo by looking at the set up in Blackbird demo.
It is all commited, alone with some other changes that do no have anythong to do, but si time to get them in.

the only part if not identical is the scale he uses, in both shape and bodies, bu I do no thing this is part for the problem
I can see teh bug whe I collide the contraction with obstacle.

The one thong we can rule out right away is the is no teh joint solver, because I have teh disable for teh work with a solver that does Implicit integrations. so tha we can have more accuracy handling gyro and Coriolis forces naturally.

The bug must be in teh collision Joint, that joint do no get much attention, and I would be surprise broke there.

One of teh problem I have wit newton is that whe I wrote the first time the major emphacis was in collision, Joints were an after thought, but then Joints start becoming more important and the data structure is the same for both contacts and bilateral, I just repurposed the variables.

The fine as long as you work and remember what you did, but after 14 years you forget and a change in one thing start causing side effects in another side, and I suspect this is what happen.

anyway I am debuging it.

[Julio Jerez]

OK I believe I found the problem.
There were actually two bugs that happen when a few moth ago I removed the exact direct solver for joint in favor of the gauss sidle only. I forget that high mass ratio is a big issue with gausssidle cenvergece.

I will put that back and add an options for these cases where, there are problematic mass ratios issues.
I though we were not going to need that since we have exact solver for joint and internal contacts, but this is not that case.

It can of course do automatic, by checking Is the island has at least one high mass ratio joint, and set the exact solver option, or I can make it and also forced by the user.
anyway I will try to bring it back from some older version first and see if it solve the bug.

Ther second bug was that the accumulation error I change from
err = fabs (accel)
to
err = accel * accel
and that is a big mistake, it should be

err += accel * accel

this is in line 555 file ..\sdk\dgPhysics\dgWorldDynamicsSimpleSolver.cpp
maxAccel += a * a;

David Gravel mentioned to me that the bug seemed random with the number of passes, he said when he made the number of passes high the bug seem random.
fixing the secong bug simple made it happen every time, because the error accumulation is more accurate and the solver makes more passes, this is what told me that this is a condition number bug that can never be solved with a Gauss or Jacobi matrix solver.

so I will bring back the direct solver again as an option.

I has to do with me removing the joint exact solve for perforemance reason.

[Sweenie]

I've been so focused on experimenting with the "bouncing" problem that I've forgot to test joints in general. Something weird is going on with them right now.
Noticed that in the standard joints scene in Unity the joints connected to world(null parent) just drops to the floor.
Then I tried the standard joints in the sandbox(reenabled the default joints) and the same thing seems to happen there, but not for all joints.

[blackbird_dream]

I found a new release on github tagged

debuging Blackbirh forklift demo.

I don't know if it is a fix but It doesn't correct the bug by me.

[Sweenie]

aha, I guess the joint solver is currently being worked on right now. Let's wait until Julio gives the go ahead to test it.

[Julio Jerez]

Yes what is committed does not work at all.
I added some debug code so that I can find the real reason of this bug.

The condition number issue is part of the problem, but that does not explain the behaviour, in fact condition number issues act to make joint softer.

David Gravel mentioned to me that he obderved that the problem seems to come from the cached forces.
He seems to be right, he is shy and almost never post, he prefer sending me private message, but his insights are very valuable and he is right almost every time.
Some people now know the engine detail better than I do.

Anyway the code that is checked in is for debugging the bug. I am ruling out cond number since this can be tested but doing unlimited passes until convergence error is very small, what is committed does that but the bug is still there. So continue debugging.

[Julio Jerez]

I've been so focused on experimenting with the "bouncing" problem that I've forgot to test joints in general. Something weird is going on with them right now.

yes this is related to the changes I am making to support gyro forces at eth solver level.
This is quiet challenging, and last Saturday, I have very encouraging result, but I still do not know if I will be able to approximate implicit derivation on the right side of the mass matrix,
worse case I can do a hybrid, where I assume some value changes are small enough the
approximation y' (t + dt) ~= y'(t) + y''(t) * dt
can be use in some places and y' (t + dt) ~= y'(t)
in some other cases.
this is why you see some joints malfunction, I commute them out until get to write the expression.

The good news is that last Saturday I run the first test when the solver was doing a two axis joints with and the wheel spun independelly in two axis with perfect conservation of momentum and without any penalty compensation.

any way I am now debugging the problem and I set the demo so that the intimidate results make sense to me. David was right, there is a bug in the cache forces.

these are the joint acceleration before the collision

0.000000 0.191282 0.168300 0.000001 0.000000 0.171406 0.000001 0.000000 0.165653 0.000000 0.000000 0.159045 0.000000 0.000000 0.152583 0.000009 0.000000 0.146362 0.000001 0.000000 0.140393 0.000004 0.000000 0.134673

0.000001 0.128661 0.080147 0.000024 0.225584 0.083925 0.000004 0.241372 0.081317 0.000001 0.234883 0.078095 0.000009 0.225744 0.074919 0.000001 0.216599 0.071866 0.000015 0.207775 0.068945 0.000001 0.199305 0.066135

0.000004 0.000000 0.027856 0.000004 0.084257 0.021304 0.000009 0.062206 0.019850 0.000002 0.057452 0.018965 0.000015 0.054825 0.018183 0.000012 0.052552 0.017434 0.000006 0.050405 0.016746 0.000015 0.048350 0.016061

0.000006 0.040799 0.000105 0.000024 0.000419 0.000002 0.000015 0.000002 0.000006 0.000015 0.000015 0.000029 0.000022 0.000018 0.000021 0.000005 0.000017 0.000022 0.000002 0.000017 0.000009 0.000013 0.000016 0.000005

these are the joint acceleration after the collision

0.000000 0.110358 0.000023 0.000000 0.000001 0.078683 0.000000 0.117226 0.011798 0.000001 0.040474 0.007478 0.000004 0.034665 0.006794 0.000019 0.032824 0.006471 0.000000 0.031437 0.006181 0.000003 0.030149 0.005939

0.000002 0.000000 0.030058 0.000003 0.140628 0.018406 0.000003 0.089813 0.016589 0.000014 0.080975 0.015785 0.000016 0.077009 0.015123 0.000000 0.073782 0.014505 0.000001 0.070763 0.013915 0.000013 0.067876 0.013348

0.000020 0.077343 0.012711 0.000011 0.063169 0.012291 0.000001 0.060006 0.011778 0.000001 0.057481 0.011300 0.000006 0.055127 0.010842 0.000004 0.052878 0.010397 0.000002 0.050722 0.009973 0.000001 0.048654 0.009566

0.000000 0.044361 0.008993 0.000003 0.044234 0.008783 0.000005 0.042877 0.008444 0.000004 0.041187 0.008103 0.000006 0.039515 0.007771 0.000000 0.037905 0.007454 0.000001 0.036360 0.007151 0.000003 0.034878 0.006863

as you can see all the acceleration are bellow the threshold, these mean that the solve we able to find a solution before and after collision, however the answer is wrong simply because the input is wrong some where.

these is the free diagram for the force a velocity of the small box, notice the velocity is reduce form
-21 to -6.0 m/s

f(74478.914063 7.232492 6.261780) v(-21.264746 -0.001171 -0.000834)
f(126650.648438 15.642369 11.401161) v(-18.161457 -0.000869 -0.000573)
f(152333.890625 -0.280977 -0.428631) v(-12.884346 -0.000217 -0.000098)
f(151849.890625 -8.846862 -6.789292) v(-6.537100 -0.000229 -0.000116)

whoever after the collision, the velocity should be a very small value and remind more of less constants during the sub steps, but it does no instead encloses monotonically and the is wrong.
and I believe the reason is that the force on the right size after collision remind too high and acct as if it was a push external force. This explain what Dave was saying that the effect seem random, because depend of the velocity on the last solution, the positive velocity will the be different.
if the velocity at the end of the last step before collision was zero, the it will look right, if for example velocity was -10, the it will bounce with a stronger force.

f(129267.429688 -31.507263 -21.148853) v(-0.209144 -0.000595 -0.000397)
f(95454.710938 -1.373553 0.618637) v(5.176999 -0.001908 -0.001278)
f(66807.406250 1.046663 1.224865) v(9.154279 -0.001965 -0.001253)
f(42558.292969 2.099880 1.421654) v(11.937921 -0.001922 -0.001202)

I still do not know why this happen, but I am sure is a bug that was introduced the pass few month.
I am investigation now.

[Julio Jerez]

I run some more test, and I determine this is a condition number bug.
I keep talking of condition number and some people do no really know what that is

essentially this is a metric that can be use to determine how sensitive a linear system of equation is to small perturbation of the input.
it is expressed as Cond(A) = abs (det(A) * abs (det(invA));

It is hard to understand what is mean for the definition and the literature by now with you tube there are intuitive and simple explanations than make sense
https://www.youtube.com/watch?v=Y-Ih1j-LA4I

why this is a problem for iterative solver? because iterate solve are base on multiples passes of the over the system each one calculation an incremental step, but because each pass rely on the value of a previous guess, you can see that if the Mass matrix A has a large condition number

the different between the privies and current step is in fact a small permutation, so when you get is a out of control magnification of small perturbations on the input.

It may be that a sufficient high condition number lead to a random error, by a medium condition number lead to a solution that has so much error that is not to be trusted.

The good knew is that there is a solution and we have before, but since we added the desire solve for Joints and no for contacts, the problem start happen again.

basically the solution to the problem is to apply a matrix preconditioner, that si intead fo solving

A * X = B

we define the identity A = C * M

where M is a matrix such that is easy to solve by direct method, a the has and C has a lower condition number that A, so now we get

C * M * X = B

let vector Y = M * X

now we solve for Y the system

C * Y = B

where now C is a more well behave matrix than A, after we get Y we can easy calculate X using a direct method.
This make the solution more accurate, if the condition number is manageable, and I already run test running an unlimited number of passes and the system do converge, so I expect the implementation will work to make it better. I will do that the Saturday,
Just remember this is not a bug, it is a fundamental problem of numerical analysis.

The reason was not there before is because the old system has that build in, but in my quest of faster solve I remove some of those checks.




The bad knew si that this is no a simple bug, I have to go over the solve

[Julio Jerez]

Ok guys I believe I have a working version now.

It was not as simple as I though it was going to be, but I think this is good enough.
I tested with mass ratios as high as 1000 to 1, I do not recomend it, but it does work.

Sweenie and Blackbird, when you read this, please try syncing an check it out before we move to Unity.
Sweenie you said you saw joints functioning, could you check that again?

what is committed has the Gyro torque disabled, let us get everything back working before with enable that funtionality.

I still have to add a function to the contact joint because it could make the solver much slower if a contraction of joints collides with too many contacts and non simulation class code may not necessarily requires that level of accuracy. The option is on by default, it needs a way to disable in a contact callback base of the type of contact.

Ok I think this is working now.

edit:
if you synced before this edit, please do it again I committed code with the hinge disable, so joint build on hinge will malfunction.
The reason that I was using eh Hinge for calibrating the Gyro functionality, and I forget to make the change under if def. I will try to be more careful about that.

[Julio Jerez]

Ok Sweeney, there is another change that will probably required compiled the Plugin

I added the function to disable contact collision using the joint solver.
but this requires change the on OnAABBOverlap interface from

Code: Select all

typedef int  (*NewtonOnAABBOverlap) (const NewtonMaterial* const material, const NewtonBody* const body0, const NewtonBody* const body1, int threadIndex);

to
typedef int (*NewtonOnAABBOverlap) (const NewtonJoint* const contact, dFloat timestep, int threadIndex);

I was goin to do that anyway, in previous version Contact Joint were no calculated unless OnAABBOverlap return true, In theory this use less memory, but after many test I realized that
this is more trouble that it is worth having.

In general we can think of the problem the same we program Rom memory.
a memory contain all possible permutations of the input, by the program ere do no have to use all combinations.
an the is how Newton 3.14 works now, Joint are created each time the broad phase determine the are close, then there is a pass that removed dead joints.
that way we get the joint even of the bodies do not collide and this simplify many complicated algorithms where a link between the bodies is needed but it does no exist yet.
therefore the joint is made, but the joint is deleted right away because not contact was generated.
Now the joint is there as long as there are close enough form the broad phase.

anyway I will link this to the Unity plugin tonight if you do no get to it first.

[Julio Jerez]

and one last check, I added the interface to deal with contacts contact to joint with high condition number.
in file ...\applications\demosSandbox\sdkDemos\demos\UsingNewtonMeshTool.cpp

function is the callback

Code: Select all

static int OnAABBOverlap(const NewtonJoint* const contact, dFloat timestep, int threadIndex)
{
//   NewtonContactJointResetSelfJointsCollision(contact);
   return 1;
}


as I said before all contacts are active by default, so the quality is always the highest possible.
The application will have to implement the On AABB in order to rest the contact base of the conditions.
The reason of this is that solve contacts as joint is quiet expensive O(n^3) when N is the number or rows. contact have many row for example that simple test will have 12 rows, therefore it si in eh best interest for the client app to be smart as to what collide what not.

here is an example:
on that test one box is very light weight and the other is very heavy.
we know that collision form eth heavy box with other object can be deal with by the iterative solver with out major problems, but when the small box collide with know that the heavy box push it, but the contact is no aware of that do the iterative solver which work by solving a joint at a time will responded with a small reaction and it will required several hundred of iteration to converge to a solution (over 1000 in this test).

here is where the engineering decision is made, The designer tag the collision shape with unique labels and in the call back check for those labels and disable any contact for which the mass is not the small body, it could also do it be check that mass, but I prefect using labels.

I know, the logic is negative, but that how all functionality is Newton is set by design.

Later when I get back to the vehicle this is the last piece that I need for the solver, to make extreme realistic vehicles.
so thanks Blackbird for prioritizing this issue.

[blackbird_dream]

doesn't compile here:

Erreur (active) E0386 aucune instance de fonction surchargée "dNewton::OnBodiesAABBOverlap" ne correspond au type requis dNewton c:\Users\me\VSProjects\NewtonDynamics\sdk\dNewton\dNewton.cpp 101

Erreur C2664 'void NewtonMaterialSetCollisionCallback(const NewtonWorld *const ,int,int,NewtonOnAABBOverlap,NewtonContactsProcess)' : impossible de convertir l'argument 4 de 'overloaded-function' en 'NewtonOnAABBOverlap' dNewton C:\Users\me\VSProjects\NewtonDynamics\sdk\dNewton\dNewton.cpp 101

version loaded:
https://github.com/MADEAPPS/newton-dyna ... e0b9cf55f0


neither does version:
https://github.com/MADEAPPS/newton-dyna ... a633a5104c

Avertissement LNK4221 Ce fichier objet ne dÚfinit pas de symboles publics jusqu'ici non dÚfinis. Par consÚquent, il ne sera utilisÚ par aucune opÚration de lien utilisant cette bibliothÞque dgCore C:\Users\me\VSProjects\NewtonDynamics\sdk\projects\visualStudio_2015_static_mt\dgGeneralVector.obj 1

same error with dAnimationPlayer.obj, dAnimationNode.obj, dAnimationClip.obj,dAnimationBlendNode.obj, dCustomControllerManager.obj

Erreur C2664 'void NewtonMaterialSetCollisionCallback(const NewtonWorld *const ,int,int,NewtonOnAABBOverlap,NewtonContactsProcess)' : impossible de convertir l'argument 4 de 'overloaded-function' en 'NewtonOnAABBOverlap' dNewton C:\Users\me\VSProjects\NewtonDynamics\sdk\dNewton\dNewton.cpp 101

[Sweenie]

Yep, the dNewton project hasn't been modified to reflect the new onBodiesAABBOverlap interface.
Julio will have to fix that first.

I made a local change for now to be able to move on and test the Unityplugin.

in dNewton.h change this

Code: Select all

CNEWTON_API static int OnBodiesAABBOverlap (const NewtonMaterial* const material, const NewtonBody* const body0, const NewtonBody* const body1, int threadIndex);

to this

Code: Select all

CNEWTON_API static int OnBodiesAABBOverlap(const NewtonJoint* const contact, dFloat timestep, int threadIndex);


and in dNewton.cpp change this

Code: Select all

int dNewton::OnBodiesAABBOverlap (const NewtonMaterial* const material, const NewtonBody* const body0, const NewtonBody* const body1, int threadIndex)
{
   dAssert (NewtonBodyGetWorld (body0) == NewtonBodyGetWorld (body1));
   dNewton* const world = (dNewton*) NewtonWorldGetUserData(NewtonBodyGetWorld (body0));
   dNewtonBody* const dBody0 = (dNewtonBody*) NewtonBodyGetUserData (body0);
   dNewtonBody* const dBody1 = (dNewtonBody*) NewtonBodyGetUserData (body1);

   return world->OnBodiesAABBOverlap(dBody0, dBody1, threadIndex);
}

to this

Code: Select all

int dNewton::OnBodiesAABBOverlap (const NewtonJoint* const contact, dFloat timestep, int threadIndex)
{
   return true;
   //dAssert (NewtonBodyGetWorld (body0) == NewtonBodyGetWorld (body1));
   //dNewton* const world = (dNewton*) NewtonWorldGetUserData(NewtonBodyGetWorld (body0));
   //dNewtonBody* const dBody0 = (dNewtonBody*) NewtonBodyGetUserData (body0);
   //dNewtonBody* const dBody1 = (dNewtonBody*) NewtonBodyGetUserData (body1);
   //return world->OnBodiesAABBOverlap(dBody0, dBody1, threadIndex);
}

Note, this isn't a correct fix, just a workaround to get the code compiling again.

The UnityPlugin was also modified to reflect this change, however I had to temporarily disable the option to disable collisions between bodies since the callback used the material to store the pointer to the dNewtonWorld. I'm not sure how to get the dNewtonWorld instance in the callback anymore.

Now the results

The joints are working again and as far as I can tell the excessive bouncing is gone.
Demo 6(connected hinge actuators) is still having problem with the joints not being as strong as before but let's leave that for later.

[blackbird_dream]

ok thks
I made the changes. I managed to compile.
The bug is still here, excessive bounce on slider.

[Sweenie]

Hmm, odd. I will try to replicate in Unity. The unity. scene you provided in this thread is missing some scripts and prefabs. It will also be interesting to see what results Dave get with the Orion engine.