Seq. impulse solver: 'wrong' velocities for SI?

[kalesony]

Hello .

I have noticed that btSequentialImpulseConstraintSolver::solveGroupCacheFriendlySetup(.) is using the previous-step corrected velocities (i.e. those stored in btRigidBody objects) to determine relative velocity and thus velocity error to be fed into SI. In Erin Catto's GDC slides describing SI the correction impulse is determined basing on the current tentative velocities (i.e. predict unconstrained velocities -> correct unconstrained velocities -> obtain constrained velocities -> update config).

What is more - contacts (btSequentialImpulseConstraintSolver::setupContactConstraint(.)) are treated in mixed fashion since the target relative velocity (restitution variable) is evaluated using old post-correction velocity but the current source/pre-correction velocity (rel_vel variable) is evaluated using the current tentative velocities as expected/described by Erin. I can see the thing being based on the very interesting thread http://bulletphysics.org/Bullet/phpBB3/ ... 146#p27918.

On one hand - both joint and contacts seem to be working "alright" using their respective approaches, but on the other: is there any particular reason for using "old" post-correction velocities for evaluating joints' current pre-correction velocity (as opposed to contacts)?

Thanks,
kalesony

[Erwin Coumans]

Your point is not entirely clear.

As far as I remember, both joints (non-contact) and contact constraints are handled the same. Penetration (error) only depends on position, so you can ignore it.

Originally, we applied the gravity/forces before going into the solver, but that was not optimal (causing jitter).
Now the forces are processed inside the solver, in solveGroupCacheFriendlySetup:

Code: Select all

			solverBody.m_linearVelocity += body->getTotalForce()*body->getInvMass()*infoGlobal.m_timeStep;
			solverBody.m_angularVelocity += (body->getTotalTorque()-gyroForce)*body->getInvInertiaTensorWorld()*infoGlobal.m_timeStep;

I recall that the btRigidBody velocity and btSolverBody velocities are not the same, because of the handling of the forces.

I discussed with Erin to process forces inside the solver, instead of applying them beforehand, and he agreed that is a better approach.

[kalesony]

Hi Erwin, thanks for taking your time to reply.

Now - I am also not sure I understood you well.

What I meant was the fact that joints/non-contact relative velocity is determined this way (btSequentialImpulseConstraintSolver::solveGroupCacheFriendlySetup(.)):

Code: Select all

btScalar vel1Dotn = solverConstraint.m_contactNormal.dot(rbA.getLinearVelocity()) + solverConstraint.m_relpos1CrossNormal.dot(rbA.getAngularVelocity());
btScalar vel2Dotn = -solverConstraint.m_contactNormal.dot(rbB.getLinearVelocity()) + solverConstraint.m_relpos2CrossNormal.dot(rbB.getAngularVelocity());

rel_vel = vel1Dotn+vel2Dotn;

i.e. velocity stored in btRigidBody objects is used which is the corrected velocity from the previous step, since as far as I could see, btRigidBody object velocities don't get updated until the solver is done.

Relative velocity for contacts is determined differently (btSequentialImpulseConstraintSolver::setupContactConstraint(.)):

Code: Select all

btScalar vel1Dotn = solverConstraint.m_contactNormal.dot(body0?solverBodyA.m_linearVelocity:btVector3(0,0,0)) 
	+ solverConstraint.m_relpos1CrossNormal.dot(body0?solverBodyA.m_angularVelocity:btVector3(0,0,0));
btScalar vel2Dotn = -solverConstraint.m_contactNormal.dot(body1?solverBodyB.m_linearVelocity:btVector3(0,0,0)) 
	+ solverConstraint.m_relpos2CrossNormal.dot(body1?solverBodyB.m_angularVelocity:btVector3(0,0,0));

rel_vel = vel1Dotn+vel2Dotn;

i.e. velocity stored in btSolverBody object is used which is the current step velocity before correction (updated precisely in the code snippet you pasted).

Erin's slides (e.g. Modeling and Solving Constraints, 2009, pages 47-56) suggest the 2nd option but it could well be that you have established an alternative version which I am simply unaware of and this what made me create this thread .

I discussed with Erin to process forces inside the solver, instead of applying them beforehand, and he agreed that is a better approach.

I am afraid I don't understand what you meant - inside the solver meaning during SI? But I can't see external forces used there - all I can see are velocity and position errors being corrected with impulses established using effective masses. Have I missed something essential?

Last but not least: I know that there has been the anti-jitter contact handling alteration about a year ago (I have linked the discussion in the original post), which explains the "funny" way of evaluating contact constraints target velocity for SI, but it does not seem to really explain the apparent inconsistency (which probably exists only in my imagination ) described above.

Cheers ,
kalesony

[Erwin Coumans]

Ah, I just noticed the rel_vel is computed a second time in setupContactConstraint (as local variable), which is likely a copy/paste bug:
the lines 571-575 should be removed:
https://code.google.com/p/bullet/source ... r=2547#536
That would make the contacts/joints consistent, and use the v_old to compute the relative velocity for both cases.

I think Erin applied the external forces (gravity) to the rigid bodies, before solving the constraints.
This could created the issue that an object resting on a table would still bounce up due to restitution, because of the gravity application.
rel_vel_old=0 but rel_vel_gravity = -9.81*time_step

Here is the commit that fixed the 'restitution' issue:
https://code.google.com/p/bullet/source ... cpp?r=2547

I agree the code is not clear, especially the difference in velocity in btRigidBody and btSolverBody. I hope to address this at some stage.

Thanks for the feedback,
Erwin

[kalesony]

Great, thanks for clearing that up .

[Edit] Argh, I am still uneasy about this thing. I re-read the anti-jitter post (http://bulletphysics.org/Bullet/phpBB3/ ... 146#p27918) and it clearly states that:

Solve constraints. The old velocity is used to compute restitution forces. the constraints however are used to constrain the new tentative velocity. A final resultant velocity is obtained and written in the object.

, which is surprisingly similar the several lines you have described as a copy-paste bug :O.

-
kalesony

[kalesony]

Ad. contact and anti-jitter fix: I did a bit of testing - using the "old" relative velocity, i.e.

Code: Select all

vel1 = rb0? rb0->getVelocityInLocalPoint(rel_pos1) : btVector3(0,0,0);
vel2 = rb1? rb1->getVelocityInLocalPoint(rel_pos2) : btVector3(0,0,0);
vel  = vel1 - vel2;
old_rel_vel = cp.m_normalWorldOnB.dot(vel);

both for evaluating the target/desired relative velocity, i.e.:

Code: Select all

restitution =  restitutionCurve(old_rel_vel, cp.m_combinedRestitution);

and the velocity error injected into the solver, i.e.:

Code: Select all

btScalar velocityError = restitution - old_rel_vel;// * damping;

leads to a roughly constant "steady-state" error/penetration, which is NOT the case if you use the new tentative/gravity-based velocity for velocity error evaluation, i.e.

Code: Select all

btScalar vel1Dotn = solverConstraint.m_contactNormal.dot(rb0?bodyA->m_linearVelocity:btVector3(0,0,0)) 
	+ solverConstraint.m_relpos1CrossNormal.dot(rb0?bodyA->m_angularVelocity:btVector3(0,0,0));
btScalar vel2Dotn = -solverConstraint.m_contactNormal.dot(rb1?bodyB->m_linearVelocity:btVector3(0,0,0)) 
	+ solverConstraint.m_relpos2CrossNormal.dot(rb1?bodyB->m_angularVelocity:btVector3(0,0,0));
new_rel_vel = vel1Dotn+vel2Dotn;
btScalar	velocityError = restitution - new_rel_vel;// * damping;

I set the reference parameter rel_vel to old_rel_vel in both cases.

If you take a look at the original solver patch-code provided by Laurent Coulon - http://bulletphysics.org/Bullet/phpBB3/ ... php?id=932 - you will see that it contains this old_rel_vel/new_rel_vel mixed approach.

-
kalesony