Sequential impulses have made the system extremely stable during stacking.

There is a small amount of sliding around though which i'm sure is a common problem.

Here' how i apply friction:

* During the warm start i cache all the data related to the contact point. Precompute all the redundant calculations.
- This includes the tangential direction, which is perpendicular to the normal in the direction of the relative velocity. Can be zero or a normalized direction
- Apply warm start impulses. Old normal impulse magnitude is applied in current normal direction. Previous friction impulse magnitude is applied in the current tangential direction. (Sounds wrong, that direction has changed.)

For each iteration
- Apply normal impulse, keep running accumulative impulse positive.
- Apply friction impulse, clamp running accumulative impulse based on friction coefficient times the current accumulated normal impulse (Sounds wrong, that magnitude has not settled yet)
- All of these friction impulses are applied in the precomputed (during warmstart) tangential direction (sounds wrong, the direction is likely changing)

The results are extremely pleasing. How ever i'm sure with more effort they could be better.

Options:
* Keep friction direction around from previous frame for warm starting.
* Handle friction in two dimensional coordinate frame relative to the body.
* Update tangential direction with each iteration

I briefly heard of central friction but haven't found many resources on it.
Maybe there are new best practices for friction?

Thank you

2D or 3D?

3D is implied. Certainly handling friction in more than one direction will make things more stable

You need indeed two friction directions then. It might help that you align the first friction axis with the relative velocity in the tangent plane. It is also important to project the last friction impulse onto the new friction directions for warmstarting. Say you have your two accumulated friction impulses lambda1 and lambda2 and the associated direction vectors tangent1 and tangent2 from the last and current frame. Then you need to do the following to project the impulse:

Vec3 OldImpulse = Lambda1 * OldTangent1 + Lambda2 * OldTangent2;
Lambda1 = dot( Lambda1, NewTangent1 );
Lambda2 = dot( Lambda2, NewTangent2 );

Don't skip friction if the relative velocity in the tangent plane is zero. You create an arbitrary frame then. You can look at dPlaneSpace() in the ODE if you need an example.

HTH,
-Dirk

Awesome, thanks guys.

Alright friction works great now.

Now my question is about "bias" velocities for penetration resolution.

I'm using basic 0.4 * (penetration - slop) / dt resolution.

Do i treat this as a parallel impulse computation to the normal impulse?
So that iterations of multiple contact constraints can adjust a bodies position bias velocity?

Or is it treated in its own unique way? Such as only considering the value from the last iteration?

Then i assume you use the combined velocity to integrate position and then zero the bias velocity. Or is that value also persistent from frame to frame/

0.4 is too aggressive for contacts from my experience. I suggest 0.1. You can look at Erin Cattos first GDC (poster) session for a derivation of this value. To avoid overshooting you can also clamp the bias velocity.

Split-impulses are another alternative. Instead of solving J*M^-1*JT * lambda = -J*v - 0.1*C/dt you split this into two separate equations J*M^-1*JT * lambda = -J*v and J*M^-1*JT * mu = -0.1*C/dt. The result of the second equation is then directly applied to the position and doesn't effect the momentum this way. Split-impulses seem to work reasonable well for contacts, bu have problems with joints since the Jacobian is not updated.

For a third alternative search in Box2D for NGS (Non-Linear Gauss-Seidel). With this approach you solve two LCPs. First you project velocities to satisfy dC/dt = J*v = 0 and after you have updated the positions you project C( x(t+dt) ) = 0. The position projection damps quite a lot. A similar effect can be seen in cloth simulation. If you want better energy conservation you need to look into symmetric projection methods. Sadly this is not so straight forward since symplectic Euler is not a symmetric integrator.

If you really want to go deep into this topic look into the books from Hairer. There are pretty advanced and expensive though.

HTH,
-Dirk

Thanks you both have confirmed my theory that i keep the penetration velocity separate and add them together during integration, then zero. I thought this process was (strangely) called bias-velocity. Bias as in separate from velocity and biased to position resolution.

I did indeed want information on split-impulses. Dirk I dont fully understand your equations. However i'm sure of what they're describing.

How does one penetration impulse interact with the others so that a global solution can converge?
Is this merely a parallel computation to the regular normal impulse?

Where i have now:
float newVelocityImpulse = ComputeVelImpulse( contact );
float lastVImpulse = contact.LastVImpulse;
contact.LastVImpulse = Max( 0.f, contact.LastVImpulse + newVelocityImpulse );
float delta = contact.LastVImpulse - lastVImpulse;
ApplyVelImpulse( contact, delta );

I would now have?
float newVelocityImpulse = ComputeVelImpulse( contact );
float lastVImpulse = contact.LastVImpulse;
contact.LastVImpulse = Max( 0.f, contact.LastVImpulse + newVelocityImpulse );
float delta = contact.LastVImpulse - lastVImpulse;
ApplyVelImpulse( contact, delta );

float newPositionImpulse = ComputePosImpulse( contact );
float lastPImpulse = contact.LastPImpulse;
contact.LastPImpulse = Max( 0.f, contact.LastPImpulse + newPositionImpulse );
delta = contact.LastPImpulse - lastPImpulse ;
ApplyPosImpulse( contact, delta );

I assume i wouldn't want to warm start these values. Warm starting with the two axis friction did not seem to work well.

For friction I am storing the world space total friction impulse from the last solution. When i compute two new arbitrary friction direction in the new frame i dot them with this vector do initialize vec2 LastImpulse for friction.

Applying that impulse at the beginning of solving seemed to not jive well, but works very well without.