# What's all this "accumulate impulses" stuff

## Recommended Posts

japro    887
Hi,

I'm working on an educational project that involves rigid body physics in 2d. It's basically working and right now I mostly try to get better stacking/overall stability. Whenever I search on the topic I always either end up redownloading the Guendelman paper or with some forum topic that says you are supposed to use accumulated impulses and refers to Erin Cattos GDC 2006 presentation. I don't seem to understand this accumulation business from the power point thingy alone and am unable to find any explanation of the concept anywhere else. I tried to stitch the whole thing together by looking at Box2D lite... I guess I'll just try to show what I did and what my understanding problems are:

Similar to what Box2D lite does, I use a SAT test and calculate those "clipped" points and save them in an arbiter object:
//simplified version....class Arbiter {//....private:  Entity *a, *b;  std::map<CPIndex, ContactPoint, CompareCPIndex> contactPoints;};

CPIndex identifies points by the Entity they belong to and the index of a point in the collision polygon. ContactPoint looks like this:
struct ContactPoint {//.... some constructors...  double normalImpulse;  double tangentImpulse;  coords2d penetrationVector, location;};

and the velocities are updated like this:
void updateVelocity(Arbiter &arb, int dt, double r, double mu){      const double tolerance = 0.01;      const double biasFactor = 0.3;      PhysicalEntity *a = dynamic_cast<PhysicalEntity*>(arb.getA());      PhysicalEntity *b = dynamic_cast<PhysicalEntity*>(arb.getB());	        if(a==0 || b==0) //check for nonphysical objects	return;      //some obvious stuff left out for...            for(Arbiter::iterator i = arb.begin();i!=arb.end();++i)	{	  if(i->second.penetrationVector == coords2d(0,0))	    continue;     /*.....     lots of calculations for relative velocities etc.     .....*/	  double delta = abs(i->second.penetrationVector);	  double vBias = biasFactor/dt * std::max(0., delta - tolerance);	  	  double kn = mAinv + mBinv + iAinv*(dotP(rA, rA) - rnA*rnA) + iBinv*(dotP(rB, rB) - rnB*rnB);	  double kt = mAinv + mBinv + iAinv*(dotP(rA, rA) - rtA*rtA) + iBinv*(dotP(rB, rB) - rtB*rtB);	  	  if(kn<10e-10 || kt<10e-10)	    return;	  	  	  double jn = (dotP(vRel, n)-vBias)/kn;	  double jt = dotP(vRel, t)/kt;	  	  double normalImpulse0 = i->second.normalImpulse;	  i->second.normalImpulse = std::min(normalImpulse0 + jn, 0.);	  jn = i->second.normalImpulse - normalImpulse0;//friction stuff	  double tangentImpulse0 = i->second.tangentImpulse;	  i->second.tangentImpulse = clamp(tangentImpulse0 + jt , mu * i->second.normalImpulse, -mu * i->second.normalImpulse);	  jt = i->second.tangentImpulse - tangentImpulse0;//end of friction stuff	  coords2d j =  r * n * jn + t * jt;	  	  a->setVelocity(vA + j*mAinv);	  b->setVelocity(vB - j*mBinv);	  	  a->setAngularVelocity(avA - dotP(j, coords2d(rA.y,-rA.x)) * iAinv); 	  b->setAngularVelocity(avB + dotP(j, coords2d(rB.y,-rB.x)) * iBinv);	}}

Now the simulation loop is something like:
for(1..n)   predictLocations();   updateCollisions();   updateVelocities(); //elasticendforfor(1..m)   predictLocations();   updateCollisions();   updateVelocities(); //inelasticendforupdatePositions()

So, one obvious difference to the code in Box2D lite is, that I update collisions between the velocity update iterations. If I don't do that I get crazy jitter. So maybe the that is already a sign that I do something terribly wrong. Also this velocity update for some reason results in no friction at all. If I replace the "friction stuff" with:
 	  jt = clamp(jt, mu * jn, -mu * jn);

it works (as in: "it looks reasonable"). But then the tangential impuls isn't used for anything?

Now, I posted all that code because there might be something very obvious that I just don't see because I already stared at it too long. Also I'd really appriciate if someone could explain what the underlying principle of that accumulation stuff is and where the warm starting comes in.

##### Share on other sites
DonDickieD    2751
The accumulated impulse becomes apparent when you study how PGS (in particular the solver from the ODE) solves a velocity level constraint system. Basically you solve J*M^-1*JT * lambda = -J * v. If you solve this system with PGS you notice that the difference to Guendelman (besides some differences in the friction model) is that instead of requiring that each delta impulse must be greater zero only the sum of all delta impulses (called the accumulated impulses) must be greater zero. This allows negative delta impulses as long as the accumulation of all impulses stays greater zero. Hence the name.

Look at Erin Catto earliest GDC presentation and you can also look at the PhD of Kenny Erleben and the master thesis of H. Garstenauer. You should be able to find everything with Google.

HTH,
-Dirk