This page contains information about the physics engine I'm writing. It's based on the paper "Non-convex rigid bodies with stacking" by Eran Guendelman, Robert Bridson and Ronald Fedkiw. You can find the paper here. Also see my other page here, as it contains source code to an earlier demo that is based on the same ideas.

It also supports the accumulated impulses scheme described by Erin Catto here in the downloads area. Penetration is handled by splitting the impulse, with the penetration-correction part not adding to the momentum.

Scroll down to the bottom of this page for some notes on the implementation. They actually refer to the versions prior to 0.61 - hopefully I'll get a chance to update the notes at some point.

You can download my test executable that uses the library here (9 Nov 2007). Extract everything before running it. It includes the source code for the complete engine licensed under the zlib license.

All screenshots below came from a slightly older version of the demo running at interactive frame rates (on my 2GHz Athlon machine). The current version looks prettier!

You can see a movie of it (version 0.70) running a wall/pyramid (128 boxes along the base) here. This was running rather slow on my machine, but I think it's rather elegant how the wall collapses!

The initial setup - all objects start active (except the ragdolls) then naturally asleep after a second or so. The box stack (behind the "jenga" stack) is 50 blocks high. There are two chains (box and sphere) in the background.

Here are some notes on the implementation that might be useful if you want to try writing something similar:

• The basic order of integration is (well... was, things have changed a bit since I wrote what's below) this:
  1. Get all forces acting on the bodies.
  2. Integrate velocity using the forces. Use the new velocities integrate positions. Detect collisions (ideally a swept test from old positions to new ones... but I currently just do an overlap test with the new positions). Restore positions and velocity to their original values.
  3. Process (see below) all collisions and constraints using impulses, iterating a few times over the list. These impulses update the velocities, but not the positions, of course.
  4. Integrate the velocities (again) using the forces.
  5. Process all collisions and constraints using impulses again, except now treat all collisions as inelastic, iterating a few times over the list.
  6. Do a "shock" step - see below. This allows objects to stack up even when the number of iterations in the previous steps is small.
  7. Integrate the positions using the new velocities.
  Work through this algorithim considering, say, a sphere sitting on a plane, and you'll see why it's good. Step 1 gets gravity acting on the sphere. Step 2 predicts that if gravity had its way, there would be a collision between the sphere and the plane. Step 3 processes the collision using the initial velocity (zero, assuming it's at rest) - nothing to be done. Step 4 adds the effect of gravity to the sphere velocity - now it's moving downward.... but Step 5 applies an inelastic collision between the sphere and the plane, zeroing the velocity. So in step 6 the velocity is (still) zero, and the sphere doesn't move.
• During the collision detection a collision record is generated for each pair of touching collision shapes. This record contains information about the bodies involved, the collision normal, and a list of points involved in the collision.
• There are different ways of processing the collisions. The most obvious is to, for each collision record, iterate over the collision points applying impulses until all points are satisfied (i.e. the two objects are separating at the collision point). However, if you imagine a box sitting on a plane, with 4 collision points at each corner of the box, then applying impulses in sequence will result in the box jittering very slightly, because the impulses at each point are affected by the impulses at the other points. It will also drift. When you have lots of boxes stacked up, this effect accumulates, and they wiggle around even more.
  One way around this is to calculate the impulse at each point, but don't apply these impulses until all points have been considered. Then apply, for example, an averaged impulse at the mean collison point. This reduces jittering a lot, but has the problem that, for a box sitting on a plane, for example, the resulting impulse is not big enough. You can see this if you imaging a really heavy box with a small moment of inertia - the impulse required at each corner to solve the corner's collision is rather small since it can be done by making the box rotate. However, when you add up all the impulses the resultant isn't enough to stop the box as a whole moving!
  I'm still experimenting for how to work around this - one option is to solve each contact-with-multiple-points using a LCP solver (so each use of the solver is with a rather small-sized set of matrices). At the moment I'm experimenting with generating an extra collision point by averaging the others - so for a box on a plane this would generate a collision point at the centre of the box's bottom face. If this point is processed before the corner points, it should generate the whole impulse required to support the box, and the corner points would just "help" to stabilise it. Work in progress...
• Friction - At each collision point the normal impulse has been calculated. Applying friction at each point is treated like a collision - calculate the impulse required to prevent the two points on the two objects diverging. Now compare this tangential friction impulse to staticFriction*normalImpulse. If it's greater, then friction won't be able to stop the objects, so apply dynamicFriction*normalImpulse in the tangential direction. If it's less then apply the originally calculated tangential friction.
• The shock step reduces the effect of only applying a few iterations in the collision processing. This is especially apparant when you get a very light body trapped between the ground and a heavy one. Firstly the bodies are sorted in the direction of gravity to help things. Then we walk through the list of bodies from bottom to top). If the body is touching something immovable (e.g. the static world, or another body marked as immovable), then all collisions between this body and other immovable bodies are processed (slightly differently to previously described - the normal impulse is applied without affecting the angular velocity). After this, the body is itself marked temporarily as immovable, and processing of the body-list continues. At the end the immovable status is restored.
  Processing things this way means that objects can "stack" in any direction - for example if you drive into a wall and squash a box between the car and the wall, the box shouldn't be pushed through the wall.
• An object deactivation algorithm is pretty essential. To see if an active object should be deactivated I use a distance method - i.e. if the object stays in the same place (also orientation), to within some tolerance, for more than a certain time, then it's a candidate for deactivation. Just prior to updating positions I look for such objects and mark them as inactive.
  Working out when to activate objects is actually quite a lot harder - you need to balance between activating things unnecessarily and having objects stay inactive when they should be activated. So, I activate objects when:
  1. After each collision/contact iteration I look at each inactive objects velocity, and if it's greater than a tolerance, I active it. Otherwise, I set that objects velocity to zero.
  2. When an object goes from inactive to active I run collision detection on it (because there might be collisions between it and another inactive object). Then for every collision, I see if the other object would move towards the object that's being woken up - if so then the other object gets woken up too (and this is done recursively - so that removing an object from the bottom of a stack makes the whole stack start to drop down immediately).
• When processing collisions/contacts I add on an extra term that would result in the penetration being resolved over a certain timescale (and I limit the relative velocity that this extra term can generate). The reason I do this instead of physically moving objects is that if you move an object during a penetration-resolution loop, you'd need to do collision detection on it again in case you move it into an object that it wasn't originally touching.
• The basic constraint used to generate hinges etc is a point-on-body-1 to point-on-body-2 constraint. This uses the same impulse equations that are used for collision handling to make the velocities at the points on the two bodies be equal. An extra term is added on to dry to eliminate divergence of position over a timescale (just like in the collision/penetration resolution step above). The other useful constraint is a maximum distance constraint - this tries to keep point-on-body1 no more than a certain distance from point-on-body2 by calculating the velocity that would be required to acheive this over a certain timescale.
• Hinges are built up from the two constraints above - my basic hinge uses a single point-point constraint to keep the objects together. Then two max-distance constraints are set up either side, along the hinge axis - the max-distance bit allows the hinge to be made slightly sloppy (i.e. bend in not quite the right direction), which is good for ragdoll joints. Finally, a max-distance constraint is set up to limit the joint movement. Overall this works pretty well, but I'm sure it's not the fastest way of doing things.
• The car just uses ray casts for the suspension (actually you can set it up to use multiple rays per wheel). Drive forces are applied at each wheel to both the chassis, the wheel, and the object that it touches. It's pretty hacky! I want to set up a system so that rather than using forces, impulses are applied during the collision processing - impulses perpendicular to the tyre would be applied to the chassis (using the chassis' mass), and impulses parallel to the tyre would use the tyre's mass/inertia.