Real-Time Physics Simulation Forum

Re: Chris Hecker's impulse calculation versus block solver

2013-05-14T06:41:37+00:00

Yeah specifically, there is a triple product identity. Turns out taking the cross-product twice is the same thing as taking the cross product and dotting it with itself. Or something. Let me see if I can figure it out:

Chris Hecker's impulse formula: (r x n) x r
Widely used impulse formula in solvers: (r x n) . (r x n) = n . ((r x n) x r) by scalar triple product identity

There it is! The jacobian entry is the normal so it is dotted with Chris Hecker's formula. Ok I don't feel stupid anymore.

Mike

Statistics: Posted by mikeshafer — Tue May 14, 2013 6:41 am

Re: Chris Hecker's impulse calculation versus block solver

2013-05-13T17:13:55+00:00

Ahh, okay, thanks. I misunderstood.

Statistics: Posted by bone — Mon May 13, 2013 5:13 pm

Re: Chris Hecker's impulse calculation versus block solver

2013-05-10T18:07:23+00:00

No no, nothing is wrong. It was my mistake. Thank you.

Mike

Statistics: Posted by mikeshafer — Fri May 10, 2013 6:07 pm

Re: Chris Hecker's impulse calculation versus block solver

2013-05-10T18:05:11+00:00

Sorry, which equation number of Hecker's is wrong? Thanks.

Statistics: Posted by bone — Fri May 10, 2013 6:05 pm

Re: Chris Hecker's impulse calculation versus block solver

2013-05-10T07:13:29+00:00

So it turns out that 2nd cross product is very necessary. It's why my boxes wouldn't stack! Ha! Go Chris Hecker! For some reason, I thought the dividing by d in PGS from Erin's paper was this calculation so I left it out, but I added it and voila, nice stability. Thanks for all your help.

Mike

Statistics: Posted by mikeshafer — Fri May 10, 2013 7:13 am

Re: Chris Hecker's impulse calculation versus block solver

2013-05-10T02:43:53+00:00

Ah I see. The following is my constraint code:

Code:

oid MEHEditorConstraint::ToContact(const MEHVector3 &overlap, const MEHVector3 &radius1, const MEHVector3 &radius2,                               const MEHVector3 &normal, scalar erp, scalar cfm){   m_jacobiVectors[0] = normal;         m_jacobiVectors[1] = radius1;   m_jacobiVectors[1].ToCrossProductWith(normal);      m_jacobiVectors[2] = normal;   m_jacobiVectors[2] *= (scalar)-1.0;   m_jacobiVectors[3] = normal;   m_jacobiVectors[3].ToCrossProductWith(radius2);     m_constraint = overlap.DotProductWith(normal);      m_erp = erp;   m_cfm = cfm;   m_calcDamping = false;   m_lowerBound = 0.0;   m_upperBound = (scalar)1e10;}void MEHEditorConstraint::ToNormal(const MEHVector3 &radius1, const MEHVector3 &radius2, const MEHVector3 &normal,                           MEHEditorRigidBody *bodyA, MEHEditorRigidBody *bodyB, scalar penetrationVelocity,                            scalar erp, scalar cfm){   m_jacobiVectors[0] = normal;      m_jacobiVectors[1] = radius1;   m_jacobiVectors[1].ToCrossProductWith(normal);         m_jacobiVectors[2] = normal;   m_jacobiVectors[2] *= (scalar)-1.0;   m_jacobiVectors[3] = normal;   m_jacobiVectors[3].ToCrossProductWith(radius2);         scalar num1 = 0.0;      num1 = m_jacobiVectors[0].DotProductWith(bodyA->GetVelocity()) + m_jacobiVectors[1].DotProductWith(bodyA->GetAngularVelocity());   scalar num2 = 0.0;      if (bodyB && bodyB->GetMass() != MEHMAXSCALAR)      num2 = m_jacobiVectors[2].DotProductWith(bodyB->GetVelocity()) + m_jacobiVectors[3].DotProductWith(bodyB->GetAngularVelocity());   scalar numerator = (scalar)(-1.0 + 0.0)*(num1 + num2) - (penetrationVelocity > 0.0 ? penetrationVelocity : penetrationVelocity * erp);        m_constraint = numerator;     m_erp = 1.0;   m_cfm = cfm;   m_calcDamping = false;   m_isContact = true;   m_lowerBound = 0.0;   m_upperBound = (scalar)1e10;}

The first constraint is Baumgarte stabilization, but you'll notice we're doing it in the contact constraint too (the 2nd function). The first function takes the overlap magnitude as the error value while the 2nd one incorporates the overlap magnitude divided by elapsedSeconds which is penetration velocity. So you're saying to put the overlap through a min(overlap + value)?

As far as the jitter goes, I thought it was because I was using GJK/EPA and not box clipping like bullet. Then I rehashed bullet's boxbox algorithm to use GJK/EPA and it can stack boxes pretty well compared to mine. So now I just got to track (and learn) what bullet is doing that I'm not. I really should give back to bullet somehow, I wish I could. I have written a Win32 editor that let's you define game objects ("widgets") in XML files and then compile the widgets into a packed data structure to be used in your game. It has a 3D view too. I could pitch it in with bullet if you guys want. Right now I'm working on a rag-doll viewer with a mouse joint. Let me know if this sounds interesting.

Mike

Statistics: Posted by mikeshafer — Fri May 10, 2013 2:43 am

Re: Chris Hecker's impulse calculation versus block solver

2013-05-10T00:11:03+00:00

I don't know why there are two cross products in Chris Hecker's paper. The formula I gave you is correct and we can solve more interesting problems instead of digging into this too much. Don't worry about it!

Statistics: Posted by Dirk Gregorius — Fri May 10, 2013 12:11 am

Re: Chris Hecker's impulse calculation versus block solver

2013-05-10T00:07:45+00:00

I think we refer to 1) as Projected Gauss Seidel (PGS) here in the forum. And 2) is refereed to as Sequential Impulses (SI). Both are equivalent iterative solvers and only differ in the implementation.

You should add gravity to your velocity before you send it to the solver. The contact solver should see the velocity! The jitter you are seeing might due to other reasons. Make sure that you don't correct the full penetration, but keep objects slightly touching. E.g. define some amount of allowed penetration (e.g. 0.5 cm) and only correct penetration if you exceed this number. Something along the lines like this:

float PenetrationToResolve = min( Penetration + cAllowedPenetration, 0 )

Statistics: Posted by Dirk Gregorius — Fri May 10, 2013 12:07 am

Re: Chris Hecker's impulse calculation versus block solver

2013-05-09T04:05:05+00:00

Ah! Well first of all thank you for your reply. I think I was getting lost in my own block solver. Basically, I look at physics calculations in two different approaches:
1) matrix style (block solver) which I got into knowing thanks to Erin's paper
2) impulse addition/subtraction (Paul Firth's speculative contact algorithm is an example of this).

Basically, it's all the same, but I didn't understand it at first. I was specifically talking about Figure 4 in http://chrishecker.com/images/b/bb/Gdmphys4.pdf He takes the cross-product twice and we don't do that with the Jacobian in the "block-solver"--I should just say "in how bullet does it".

While I have your attention, I got another question for you: I moved my gravity calculation out of a pre-solver step and put it into the solver as "F_ext" as Erin describes in Iterative Dynamics. Then I put an if statement for contact points: if there is no contact manifold on the rigid body, then apply gravity to velocity outside of the solver (otherwise nothing moves if it's not colliding yet). This seems to almost work as my boxes stack okay now (they weren't able to before), but the whole stack wobbles continuously. Do you have any thoughts as to what I could do to fix this? BASE_CONTACT_BREAKING_THRESHOLD = 0.02 which is the equivalent of m_contactBreakingThreshold in bullet. If I make this smaller the wobbling is less noticeable, but manifolds aren't as persistent anymore and it looks like it's just wobbling at a higher frequency. If I don't do the gravity application in the solver and keep it as a pre-step, then no box stacking, but no wobbling either. Let me know your thoughts. I would be so appreciative if you could shed some light on this.

Mike

Statistics: Posted by mikeshafer — Thu May 09, 2013 4:05 am

Re: Chris Hecker's impulse calculation versus block solver

2013-05-09T02:26:02+00:00

It doesn't matter how you compute the effective mass. You hopefully get the same result with both computations. Otherwise one is wrong. The inverse effective mass for the non-penetration constraint is:

K = invM1 + invM2 + (r1 x n)^T * invI1 * (r1 x n) + (r2 x n)^T * invI2 * (r2 x n)

The reason the Bullet solver works better is because he clamps against the accumulated impulse. Instead of just applying only positive impulses.

I am not sure I understand what you mean with block solver. A block solver usually means that you solve n-constraints as a block. E.g. you have a hinge which is defined by 5 constraints. Instead of solving 5 scalar equations you solve a 5x5 linear system. The 5x5 inverse effective mass will be symmetric positive definite.

Statistics: Posted by Dirk Gregorius — Thu May 09, 2013 2:26 am