Chris Hecker's impulse calculation versus block solver

Please don't post Bullet support questions here, use the above forums instead.
Post Reply
mikeshafer
Posts: 49
Joined: Fri Feb 24, 2012 10:45 am

Chris Hecker's impulse calculation versus block solver

Post by mikeshafer »

Hi there,

I was hoping someone could answer a simple question for me. Chris Heckers GD article is widely referenced as a starting point for building your own physics engine. I know I looked at the article and it helped me a great deal when I was first starting out. Now I got a GJK/EPA collision engine coupled with dynamics that resemble bullet's block solver. It works well. Now, I would like to look back and say Chris Hecker's method is like bullet, only not block-solver, but something baffles me. In that article, we take the cross product twice to what should be the Jacobian, but I just checked my block solver and we only take the cross product once. Is this the reason why the block solver method performs better? I think there are a multitude of reasons why my engine performs better without the "jiggle", including help I have received from this forum and bullet's methodology. But yeah, that double cross-product doesn't seem to be anywhere. Anyone can shed some light on this?

Mike
Dirk Gregorius
Posts: 861
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA

Re: Chris Hecker's impulse calculation versus block solver

Post by Dirk Gregorius »

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.
mikeshafer
Posts: 49
Joined: Fri Feb 24, 2012 10:45 am

Re: Chris Hecker's impulse calculation versus block solver

Post by mikeshafer »

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
Dirk Gregorius
Posts: 861
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA

Re: Chris Hecker's impulse calculation versus block solver

Post by Dirk Gregorius »

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 )
Dirk Gregorius
Posts: 861
Joined: Sun Jul 03, 2005 4:06 pm
Location: Kirkland, WA

Re: Chris Hecker's impulse calculation versus block solver

Post by Dirk Gregorius »

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! :)
mikeshafer
Posts: 49
Joined: Fri Feb 24, 2012 10:45 am

Re: Chris Hecker's impulse calculation versus block solver

Post by mikeshafer »

Ah I see. The following is my constraint code:

Code: Select all

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
mikeshafer
Posts: 49
Joined: Fri Feb 24, 2012 10:45 am

Re: Chris Hecker's impulse calculation versus block solver

Post by mikeshafer »

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
bone
Posts: 231
Joined: Tue Feb 20, 2007 4:56 pm

Re: Chris Hecker's impulse calculation versus block solver

Post by bone »

Sorry, which equation number of Hecker's is wrong? Thanks.
mikeshafer
Posts: 49
Joined: Fri Feb 24, 2012 10:45 am

Re: Chris Hecker's impulse calculation versus block solver

Post by mikeshafer »

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

Mike
bone
Posts: 231
Joined: Tue Feb 20, 2007 4:56 pm

Re: Chris Hecker's impulse calculation versus block solver

Post by bone »

Ahh, okay, thanks. I misunderstood.
mikeshafer
Posts: 49
Joined: Fri Feb 24, 2012 10:45 am

Re: Chris Hecker's impulse calculation versus block solver

Post by mikeshafer »

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
Post Reply