Baltman's paper on Verlet integration

[kkw]

Greetings,

Link to the paper http://www.icvonline.com/gdc_dir/files/ ... ration.pdf

I'm new to rigid body dynamics and have been toying around with Rick Baltman's paper on "Verlet Integration and Constraints in a Six Degree of Freedom Rigid Body Physics Simulation" for the past two weeks. I pretty much followed almost all the formulas written in the paper (haven't done the collision part though), and I have the basic working. Notice I bolded the word almost.

When I was doing the angular constraint. I noticed that q_correction resulted from equation 8.6 is not correct. I noticed that if minimum angle is violated, q_correction will rotate the bodies so that both bodies will violate or almost violate the maximum angle (also vice versa). For example, minimum angle of x-axis is -30, current x-axis angle between bodies is -31 (constraint violated), the resulting q_correction rotates both bodies so that their x-axis angle is around +29, about 60 degree rotation (The visual result from this is sudden snap. Imagine two rigid bodies with pin and angular constraint shaped like ">" and on the next frame it becomes "<"). I'm not an expert in quaternions and its intricate operations, I have checked my code over and over again, I am sure my code is an exact translation of the formulas in the paper.

Then I changed equation 8.4 in my code, so that it looks like this (only put x-axis, should be similar for y & z-axis).

Code: Select all

if ( dq_a.x < q_min.x ) { dq_a.x = q_min.x - dq_a.x }
else if ( dq_a.x > q_max.x ) { dq_a.x = q_max.x - dq_a.x }

I calculate the difference between min or max angle with the current angle, instead of setting dq to min or max angle. Also, I changed second part of equation 8.6 from

Code: Select all

q_correction = dq * dq' 

to

Code: Select all

q_correction = dq'

and it correctly work as expected!

So, I was wondering if anyone has ever tried to implement the paper fully ? Did you get the same result as I did ? I just would like to confirm if my fix is correct and that I didn't misunderstood the formulas in the paper. Any help is very much appreciated!

Thank you very much!

[Dirk Gregorius]

Where did you dig this out? I would argue this method is pretty obsolete. Why not using Projected-Gauss-Seidel (PGS) or Sequentiel Impulse (SI) on the velocity level with either Baumgarte stabilization or post-projection? Post projection works on the position level as well after solving for the constraint forces/impulses. Note that PGS and SI are mathematically equivalent, but only differ in the implementation. An example of the former can be found in the ODE and an example of the later here in Bullet. Both use Baumgarte as stabilization method. For a discussion on post-projection search this forum for non-linear Gauss-Seidel (NGS)...


HTH,
-Dirk

[kkw]

Hello again,

@Dirk - Well, I followed your advice . Thanks for the pointer!

For starter, I am taking Box2D(V.1.2.1) solver code and directly translate it to 3D line by line (literally ). Currently, I'm doing the revolute joint with limit in three axis (similar to Generic6DOF in bullet). So far, the joint seems to be very springy, especially at high speed, the constraint is satisfied very slowly over a few seconds (60fps). Using it for ragdoll especially visible, legs and hands bend unnaturally.

I use quaternion for orientation, convert them to euler angle when testing for limit violation and solve each axis separately. Constraint iteration is 10, and only 1 simulation loop with dt of 1/60.

When I run Box2D demo called "motors and limits", the limit is very solid! The 3D code and Box2D(V.1.2.1) code should be analogous, I can't point out what's wrong. I've read a few forum topics regarding stiffness but still at lost, it seems to be the case for many ?

Also, can split impulse be used for joints ? Box2D lite that comes with Erin's GDC2007 presentation only uses it for contact, but not for joint.

Any pointer in reducing the springyness is very appreciated, thanks a bunch!

[Dirk Gregorius]

Please make sure that both the spherical (ball-socket) and revolute (hinge) are working properly without any limits or motors. I assume you use Baumgarte stabilization - that is the bias in Erin's presentation. Do you compute it from the positon constraint and scale it by 0.1/dt -> bias = 0.1 * C / dt? Besides all this you should start with gravity -10 and a fixed timestep of 60 Hz. Also make sure that your contacts and the contact caching work properly. Do you use the slob for the allowed penetration? You should be able to stack 3-5 boxes absolutely stable with mass 1kg and dimension of 1x1x1 m droped onto each other from a small height difference. Does all this work properly? If yes what are your constraint fucntions C for the spherical and revolute joint?

Cheers,
-Dirk

[Dirk Gregorius]

BTW: There is also a Box2D_Lite. It might be easier in the beginning to start with this. Once you have tthis running you can swith to the newer version. I would recommend this...

[Dirk Gregorius]

Regarding split impulse:

Personally I would recommend to either use Baumgarte stabilization or what is sometimes called post-projection. Post-projection is similar to the Jacobsen method (you maybe heard about it) as a second sweep over the constraints, but correcting positions instead of velocities now. In my opinion the spit impulse is only a hack which actually doesn't work for joints. It seems to work good for contacts though. You find an example of the post-projection in the newer Box2D version. Using Baumgarte stabilization is the fastest you can do, while post-projection gives you much better results, especially with joints, but it is slower. I would implement both in the end and start with Baumgarte because of simplicity.

[Erin Catto]

Yes, the split-impulse method is now obsolete. It is a hack with a shakey mathematical foundation. Baumgarte stabilization has a more solid mathematical foundation and is still quite useful in some situtations. However, Baumgarte stabilization has the side effect of changing momentum and makes contacts springy even with zero restitution.

I'm now recommending a velocity correction phase using Projected Gauss-Seidel via impulses with no position stabilization followed by post-projection of position errors using nonlinear Gauss-Seidel via position impulses. IMO, this is the best method known so far.

I would call this the post impulse (PI) method. PI can have performance almost as good as Baumgarte because you can reduce the number of iterations and also you can exit the post-projection iterations early if the position error is below a given tolerance. In my tests, early exits are very common.

[Dirk Gregorius]

Maybe for clarity you should add that a position impulse is basically a mass weighted displacement or rotation (not effecting the momentun)...


BTW: This method is similar to the post-projection e.g. in the Cline thesis and I think also Kenny mentions something like this. One thing I see is that we could add potential energy when moving the bodies. I would argue that this is only a theoretical problem. Can you imagine problems because of this?

[Erin Catto]

Yes post-projection is often mass weighted but it doesn't have to be. The reason many implementations use mass weighting is because they probably already have code/data to handle the dynamics problem that can be re-used to handle position correction.

Unlike Baumgarte stabilization, post-projection does not affect the momentum.

The reasons to use post-projection are:
1) doesn't affect momentum
2) fixes overlap in the current step (not the next step like Baumgarte).
3) lets you handle convergence differently thant the velocity constraints.

On top of this, the post-impulse method has these advantages:
1) you can account for nonlinearity in the position constraints
2) you can treat each constraint differently. For example, with joints you might compute the new nonlinear position error exactly for each iteration, but for contacts you can estimate the new penetration error instead of re-generating the contact points.

I looked at Cline's paper and Kenny's dissertation. They are both stating the same problem, but it looks like the solution is not equivalent to nonlinear Gauss-Seidel. They both don't update the nonlinear position error though the iterations. Basically they are doing one Newton step. Cline uses a pivoting method and I'm guessing Kenny uses PGS.

I think the closest methods would be Jakobsen's stuff or Jan's impulse method. The main difference with Jakobsen's method is that PI has a Lagrange multipliers that can be useful in some situations.

I've seen post projection overshoot contacts, so it's not a silver bullet for inequality constraints. But it seems to do very well with bilateral constraints.

[sbroumley]

How do stacks handle using this method? Since you are updating the positions and not the velocity does it cause the bottom blocks to be squeezed out? Are you going to release a version in Box2D?

thanks

[Erin Catto]

Stacking and joints work much better with post-impulses than with Baumgarte.

I'm currently using this algorithm in Box2D (http://www.box2d.org). Unfortunately I have completely removed the old Baumgarte stuff so you can't do any side-by-side comparisons. I did the comparison while prototyping and found post-impulses to work better in all cases.

[kkw]

Whoa, thanks for all the replies!

Please make sure that both the spherical (ball-socket) and revolute (hinge) are working properly without any limits or motors. I assume you use Baumgarte stabilization - that is the bias in Erin's presentation.

Yes, revolute joint works properly without limits or motors. Not using Baumgarte, but using what Erin would call Post impulse (PI). That is what Box2D(v.1.2.1) code has. As I have said, the Box2D code and what I have should be similar, I convert Box2D code to 3D while trying to understand the math behind. So, basically I should get the same result, unless I'm doing something wrong (which is I'm very sure about) ... I don't have collision & contact yet, because I want to quickly test joint limit before going to collision & contact.

PI can have performance almost as good as Baumgarte because you can reduce the number of iterations and also you can exit the post-projection iterations early if the position error is below a given tolerance. In my tests, early exits are very common.

I have a question about the code in Box2D regarding the early exits.

Code: Select all

bool jointsOkay = true;
for (int i = 0; i < m_jointCount; ++i)
{
	jointsOkay = m_joints[i]->SolvePositionConstraints();
}
if (contactsOkay && jointsOkay)
{
	break;
}

In the above code (taken directly from Box2D, b2Island.cpp), let's say we have 2 joints. First joint returns false, then the second one returns true. This would result with an early exit for all the joints, but then the error in the first joint is still too large. Is this intended ? If I remove the break, I notice the joints tend to satisfy a bit faster when limit is violated, but still not as fast as I would like.

[Erin Catto]

Thanks for letting me know. That is obviously a bug. It's a bit surprising that it didn't cause a poor effect. I guess one iteration on joints was enough in many cases.

I uploaded a fix to the SVN repository.

If you find anymore issues, please let me know at http://www.box2d.org/forum

[kkw]

Thanks for letting me know. That is obviously a bug. It's a bit surprising that it didn't cause a poor effect. I guess one iteration on joints was enough in many cases.

The least I can do for an already great piece of software .