Accumulating Impulses, how do I store this data?

[arcoraven]

I'm implementing an iterative SI rigid body solver in 3D similar to Guendelman's paper. Simulation results are working nicely, though I'm getting some weird jumping/unexpected large energy gains from thin objects (like long, thin boxes or cylinders with very small height). I was hoping implemented accumulated impulses would resolve this.

First, let me know if I understand the "clamping the accumulated impulse" technique correctly (no warm-starting just yet):
At the start of each iteration in the collision stage, I run broad and narrow phase tests to identify a possible collision. I create a Collision object that stores collision data between the two objects. Since I know which points on each object are involved in the collision, I create an array of collision points, which are all the points which should have an impulse applied to it. For each point, if this is the first time in this timestep that this point is involved in a collision, set the point's accumulated impulse to zero. I compute a correctional impulse needed and clamp the accumulated impulse, resulting in a modified accumulated impulse and (maybe) modified correctional impulse (this part I adapted from Box2D's implementation). In subsequent iterations, I repeat the above but use the stored accumulated impulse value.

My question then is this:
Accumulated impulses are attached to each collision point. In my sim, my collision points are simply represented by an X, Y, Z coordinate in world space. How do I save an accumulated impulse for each collision point? In one iteration, a box may collide on an edge with points p0 and p1, but in the next iteration, it may be colliding on just one vertex p1. If I saved a vector<double> of accumulated impulses to correspond to the points on the first iteration, it would be unusable on the second iteration. How am I supposed to save accumulated impulses across several iterations, if collision points come and go?

Thanks, this forum has been an invaluable resource for me.

[Dirk Gregorius]

Did you slowly ramp the impulses up when you tried the Guendelman paper? The paper basically suggests to go from 0.1 to 1.0 in 10 iterations.

The accumulated impulse is stored with each contact point. If a point changes dramatically, e.g. if you bounce from one vertex to another in between two frames you have to clear the impulse. The whole art for robust contact point creation is to create coherent points (if possible) so you can match impulses and associate the last impulses with the new contact points. Look at Box2D or even simpler at Box2D Lite to see how this contact matching works. If you take an incremental manifold approach like Bullet the matching is actually straight forward since you store the accumulated impulse with the persistent points.

HTH,
-Dirk

[arcoraven]

Thanks Dirk, I ended up throwing away the accumulated impulse in the iteration if the number of contact points changed between iterations within a single timestep. Weird jumpiness was resolved for most parts, but still evident in some places. I'll work on it more some other time.

Rather than creating another topic, I was wondering if you (or anyone else) can help me on my shock propagation implementation. I created a contact graph (kind of, I just have a list of the objects sorted by ascending y-coordinate for now) and am using this list order for the contact stage (10 iterations) and the shock propogation step (1 iteration). Should I be using it for the collision stage as well?

Anyway, my pseudocode looks like this:

Code: Select all

level = 0
for each collision C in the collisionList (which was created in the order of the contact graph list, and thus should be resolved in order)
	if(rigidBody[C.i].contactGraphLevel <= level && rigidBody[C.j].contactGraphLevel <= level) {
		if(rigidBody[C.i].contactGraphLevel < level) {
			rigidBody[C.i].mass = 0
			rigidBody[C.i].massInv = 0
		}
		if(rigidBody[C.j].contactGraphLevel < level) {
			rigidBody[C.j].mass = 0
			rigidBody[C.j].massInv = 0
		}
		resolveCollision(C, restitution = 0) // perfectly inelastic
		C.valid = false // make sure I don't resolve this collision again
	} else {	// done with all objects on current level, move on to an upper level
		level++
		c--		// re-process this collision with the new level value
	}

Any help is greatly appreciated. Thanks.

[Dirk Gregorius]

I cannot help with shock-propagation, but I have seen plenty of posts about it on GameDev.net. Maybe search the forum there and also here.

Anyway, I strongly recommend to get the accumulated impulse working. If you see strange bouncing this is very likely a bug in your implementation. You don't need to warmstart in the beginning and maintain the accumulated impulse over several frames. Just initialize the impulse with zero before you enter the iterative solver loop. I am sure you know, but just for completeness a sample how to clamp the impulse:


float delta_impulse = -effective_mass * velocity_error;

float impulse = contact_point.mImpulse;
contact_point.mImpulse = Max( impulse + delta_impulse, 0 );
delta_impulse = contact_point.mImpulse - impulse;

Apply( delta_impulse );



HTH,
-Dirk

[arcoraven]

Regarding shock propagation, I've decided to scrap it after reading several threads like this one:
http://www.bulletphysics.org/Bullet/php ... ?f=9&t=368
My stack of 5 boxes is unstable at anything less than 40 contact iterations at 24Hz. It appears the impulses that I'm applying for each collision isn't strong enough, and the boxes sink into one another. I'll look into fixing my normal impulses first to see if that helps.

Yes, your code is pretty similar to what I have. I'm confused about this line, however:

Code: Select all

float delta_impulse = -effective_mass * velocity_error;

I followed Catto's Box2D Lite closely and don't see that mentioned in the Arbiter code. What are those two terms on the right side? My delta_impulse is the computed impulse that's seen everywhere (J = -(1+e)*vn / blah blah). Is this equivalent to what you had in mind?

A simplified pseudocode for my collision/contact steps looks something like this:

Code: Select all

create a contact graph: a list of the objects sorted along the -Y (gravity) axis
for each iteration (5 for collision, 10 for contact)
     create an empty collision list L
     identify all collisions in the current timestep, getting collision points, collision normal, and penetration depth
          note:  if a collision does not exist in L, populate its data and set accumulated impulse values at each point to zero
                  if the collision does exist in L, only update the collision details and leave accumulated impulse values alone
     for each collision C
          for each collision point P
               if the objects are non-separating at P (relative velocity < 0)
                    calculate the impulse J
                    clamp to the accumulated impulse at J (with the middle 3 lines in your code)
                    apply J to object 1, apply -J to object 2

Since I'm starting the collision list empty each time step, I'm assuming that means I'm not warm starting.
Perhaps my problem has something to do with the way I'm storing the contact points. I was a little confused by your earlier post:

The whole art for robust contact point creation is to create coherent points (if possible) so you can match impulses and associate the last impulses with the new contact points.

Do you mean that instead of storing the 3D coordinates of the contact points, I should be storing indices of contact points? That way I can match accumulated impulses to points? E.g. accumulatedImpulse[6] is the accumulated impulse for objectA.vertex[6]? I thought about this, but how would I store accumulated impulse values for points that I compute "per-collision", such as the intersection of two objects' edges? (It's hard to explain...I hope I'm getting my point across)

Last question: in Box2D Catto applies the normal impulse before even computing the friction impulse. I'm computing both together and applying the vector sum at the end. Is that okay?

Again, thanks a lot for your help, Dirk.

[Dirk Gregorius]

I might find time later this day to answer more of your questions. Here two quick thoughts:

1) Drop the Guendelman scheme and just iterate over the contacts applying sequential impulses (see Bullet or Box2D). Also use a frequency of 60 Hz for now.
2) I think you have to compute the friction impulses separately as in Box2D. Otherwise I dont see how you could associate the delta impulses with the correct directions for accumulation.

In general: Start as simple as possible with proven algorithms and then change from there.


Cheers,
-Dirk

[arcoraven]

I took your advice and decided to start simple.
I scrapped all the Guendelman features (using predicted positions, separate contact step with -1...0 restitution, shock propagation, etc.). I'm following Box2D's time stepping method very closely now, and think I have a good way to caching 3d contact points.

Results are mixed. At 60Hz and 10 iterations, I'm getting good speeds, but I'm getting much more jitter than before. There's still some fairly frequent odd behavior as well, so I'm certain my code isn't 100% correct yet anyway.

A few questions about Box2D that I have:
- What are the variables involved in line 121 of Arbiter.cpp:

Code: Select all

c->bias = -k_biasFactor * inv_dt * Min(0.0f, c->separation + k_allowedPenetration);

It looks like Box2D doesn't use restitution, and thus the numerator looks different than the impulse equation I'm used to. Is it okay to continue using -(1+restitutionCoefficient) * vDotN for the numerator of the impulse equation?
- What method is recommended to remove jitter and deactivate/freeze objects in an SI implementation? I can't compile Box2D Lite on the machine I'm on, but browsing through the code I don't see any deactivation/freezing mechanism.

Thanks again.

[Dirk Gregorius]

Deactivation/Freezing/Sleeping is used to save CPU time and not to hide problems of your solver.

I don't remember whether Box2D *Lite* can handle restitution. Box2D can for sure. Actually this is quite simple and can be added later (just 5-10 lines of code). Get things working without restituion first.

The line of code you quoted basically allows the shapes to slightly sink into each other. This helps to make contact points more persistent over several frames. Otherwise you would loose the contact points and more importantly the associated accumulated impulse. This is important for stable stacking. All engines I am aware of use some kind of "skin".

Regarding the variables:
k_biasFactor is the Baumgarte stabilization factor. Read e.g. Erin's first GDC presentation for an explanation. Usually 0.1 - 0.2
k_allowedPenetration is the described "skin" I mentioned above. Usually 0.01 - 0.02 (1-2 cm)

I compiled some very basic demos which you can use as reference and orientation for your own solver. The demos can be downloaded here: www.dondickied.de/download/demos.rar

Controls review:
- ALT + Mouse camera (similar to Maya)
- RMB grap body
- R reset scene
- P pause / unpause
- SPACE single step if paused
- Arrow up/down iteration through scenes

There is also a small GUI. Continuous simulation is not working yet (work in progress). Also notice the friction artifacts in the pyramid if you disable warmstarting. I use a different friction model sometime refered to as central friction, but you should get similar if not even better results with per point friction.

HTH,
-Dirk

[arcoraven]

Thanks for the explanation, and your demos. I will test it when I get home to my windows machine.
I'm using the bias and slop factors now, and everything looks good. I tried adding restitution, but didn't get the result I wanted. I'm guessing it's not as easy as just adding (restitutionCoefficient * -vDotN) to the bias, is it?

Stacking boxes works nicely if I crank up the iterations. A 6-box stack runs at 60fps with 60+ iterations. Accumulated impulses are warm-starting are implemented (though the latter hasn't been thoroughly tested...)

I read one of Catto's posts (either on the bulletphysics or box2d forums) that he uses a perfectly inelastic restitution of 0 for low-velocity collisions (not uncommon in other physics models). I set my bias, then, to something like:

Code: Select all

bias = (vDotN > LOW_VELOCITY_THRESHOLD ? 0.2 * timestep_inverse * min(0.0, 0.01 - separation) : 0.0);

(note: my separation is positive for interpenetrating objects)
So that I would use a bias of 0 for low velocity contacts and the typical bias for all other contacts.

The problem is that in my implementation now, I'm allowing objects to sink into one another for several timesteps (as you suggested). This means that ALL collisions, such as dropping a cube onto a plane, will undergo these "low velocity contacts" and be treated as inelastic. The end result is that all collisions end up looking inelastic. Am I going about implementing this feature the right way? On the bright side, stacking is much more stable, and I'm able to stack more with fewer iterations, so I have high hopes for it.

One last bit I'm confused on (sorry for all the questions!): am I supposed to skip contact points that are already non-separating? Many papers mention that if vDotN > 0, skip and continue to the next contact point. But in Box2D (if I understand correctly) I don't see this anywhere. Or perhaps we need this for accumulated impulses, because correctional impulses can negative as long as the accumulated impulse is positive.

Thanks again. I'll be sure to check out your demo.

[Dirk Gregorius]

One last bit I'm confused on (sorry for all the questions!): am I supposed to skip contact points that are already non-separating? Many papers mention that if vDotN > 0, skip and continue to the next contact point. But in Box2D (if I understand correctly) I don't see this anywhere. Or perhaps we need this for accumulated impulses, because correctional impulses can negative as long as the accumulated impulse is positive.

No, this is the whole idea of the accumulated impulse. Instead of requiring that each *individual* impulse is greater zero you require now that the *sum* of all impulses is greater zero. This way you allow negative delta impulses as long as the sum remains positive. This is actually the major difference to the impulse-based method from e.g. Mirtich or Guendelman. The code example I posted earlier takes care of exactly this clamping.

I'm using the bias and slop factors now, and everything looks good. I tried adding restitution, but didn't get the result I wanted. I'm guessing it's not as easy as just adding (restitutionCoefficient * -vDotN) to the bias, is it?

Here is what I use:

Code: Select all

// Penetration recovery and restitution
const plReal kMinPenetration = 0.01f;
out.mBias = k_cor * plMath::Min( separation + kMinPenetration, 0.0f );
out.mEffectiveMass = 1.0f / K;
out.mLambda = lambda;

plVector3 v_rel;
plComputeRelativeVelocityAt( v_rel, pBody1, r1, pBody2, r2 );

plReal vn_rel = plVec3Dot( v_rel, n );

const plReal kMinVelocity = 1.0f;
if ( vn_rel < -kMinVelocity )
	{
	plReal ve = restitution * vn_rel;
	if ( ve < out.mBias )
		{
		out.mBias = ve;
		}
	}

Stacking boxes works nicely if I crank up the iterations. A 6-box stack runs at 60fps with 60+ iterations. Accumulated impulses are warm-starting are implemented (though the latter hasn't been thoroughly tested...)

Five boxes should be stable between 4-10 iterations. I am not 100% sure, but you can check this in my demos.