2d - How to make a stack stable? Need help for an explicit resting contact scheme (2-dimensional)

Previously, I struggle with the sequential impulse-based method I developed. Thanks to jedediah referring me to this paper, I managed to rebuild the codes and implement the simultaneous impulse based method with Projected-Gauss-Seidel (PGS) iterative solver as described by Erin Catto (mentioned in the reference of the paper as [Catt05]).

So here's how it currently is:

1. 
   

   The simulation handles 2-dimensional rotating convex polygons.

   
   


2. 
   

   Detection is using separating-axis test, with a SKIN, meaning closest points between two polygons is detected and determined if their distance is less than SKIN.

   
   



3. 
   

   To resolve collision, simultaneous impulse-based method is used. It is solved using iterative solver (PGS-solver) as in Erin Catto's paper. Error-correction is implemented using Baumgarte's stabilization (you can refer to either paper for this) using J V = beta/dt*overlap, J is the Jacobian for the constraints, V the matrix containing the velocities of the bodies, beta an error-correction parameter that is better be < 1, dt the time-step taken by the engine, and overlap, the overlap between the bodies (true overlap, so SKIN is ignored).

   
   

However, it is still less stable than I expected :s I tried to stack hexagons (or squares, doesn't really matter), and even with only 4 to 5 of them, they would swing! Also note that I am not looking for a sleeping scheme. But I would settle if you have any explicit scheme to handle resting contacts.

That said, I would be more than happy if you have a way of treating it generally (as continuous collision, instead of explicitly as a special state).

Ideas I have tried: Using simultaneous position based error correction as described in the paper in section 5.3.2, turned out to be worse than the current scheme.

If you want to know the parameters I used:

• 
  
  

  Hexagons, side 50 (pixels)

  
  


• 
  

  gravity 2400 (pixels/sec^2)

  
  


• 
  

  time-step 1/60 (sec)

  
  


• 
  

  beta 0.1

  
  
  


• 
  

  restitution 0 to 0.2

  
  


• 
  

  coeff. of friction 0.2

  
  


• 
  

  PGS iteration 10

  
  



• 
  

  initial separation 10 (pixels)

  
  


• 
  

  mass 1 (unit is irrelevant for now, i modified velocity directly<-impulse method)

  
  


• 
  

  inertia 1/1000

  
  

Thanks in advance! I really appreciate any help from you guys!! :)

EDIT

In response to Cholesky's comment about warm starting the solver and Baumgarte: I forgot to mention! I do save the contact history and the impulse determined in this time step to be used as initial guess in the next time step.

As for the Baumgarte, here's what actually happens in the code. Collision is detected when the bodies' closest distance is less than SKIN, meaning they are actually still separated. If at this moment, I used the PGS solver without Baumgarte, restitution of 0 alone would be able to stop the bodies, separated by a distance of ~SKIN, in mid-air! So this isn't right, I want to have the bodies touching each other. So I turn on the Baumgarte, where its role is actually to pull the bodies together! Weird I know, a scheme intended to push the body apart becomes useful for the reverse.

Also, I found that if I increase the number of iteration to 100, stacks become much more stable, though the program becomes so slow.

UPDATE

Since the stack swings left and right, could it be something is wrong with my friction model?

Current friction constraint: relative_tangential_velocity = 0