Move over applet to start simulation

To bounce the particles off the edges we simply detect collisions against the edges and reverse the velocity in the offending axis, thus:

  if (P[i].x < -1 && V[i].x < 0)
  {
    V[i].x = -V[i].x;
  }
  if (P[i].y > 1 && V[i].y > 0)
  {
    V[i].y = -V[i].y;
  };
  if (P[i].x > 1 && V[i].x > 0)
  {
    V[i].x = -V[i].x;
  }
  if (P[i].y < -1 && V[i].y < 0)
  {
    V[i].y = -V[i].y;
  }

The first condition of each case checks for intersection of the particle in space and the second checks that the particle is actually heading towards the edge. The second condition is important otherwise we would do the collision response again next frame erronusly should the particles move more than one pixel per frame. This second condition carries though all of rigid body dynamics and is what separates inequality contraints (contacts) and equality constraints (joints).

This is the kind of example code you will no doubt have written a long time ago when you first started getting interested in simulation. I'm including it here because I think there is a natural progression from this easy example to more complex physics code.

At this point i should probably talk about the ordering of integration and collision detection i'm going to be using here. The traditional method is to add forces (i.e. gravity) to velocity and then add velocity to position and then do collision detection and response. The trouble with this is that this method will always over-estimate the time of contact and will allow penetrations to creep in. A better method is to add forces to velocity (as before) but instead of adding velocity to the position we take a copy of the position and add the velocity to that, then do the collision detection against this new temporary position, resolving the velocity to post contact velocity. Finally we integrate by adding the new post collision velocity onto the original position. This method will actually under-estimate the time of contact, which means that objects will hit the surface slightly before they should have, but the stability it adds to the simulation far out weighs any inaccuracy in the detection. This method is that used in the paper Nonconvex Rigid Bodies with Stacking by Guendelman et al. 

Whats going on in this simple simulation?

Well, without knowing it, we have assumed a physical matieral type for our particles including the coefficient of restituition and followed Newton's principle of conservation of momentum . We have also chosen that the world have infinite mass such that it doesn't feel a collision response when the particles impact on it. By choosing to reflect particles on impact we are also preseving their momentum (even though we have ignored their mass in our calculations), indicating that we have chosen a coefficient of restitution of 1, i.e. perfectly elastic - like a super-ball.

What we have here is actually a highly optimised special case physics simulation. Optimised how, you ask? Well, let me explain:

If we were to code the above "properly", not taking any short cuts, we would have to assume the following. Our environment is defined by its boundaries which are four axis aligned planes (actually, they're lines because we are in 2d). Each one has a normal (which points inwards) and a distance to the origin. They look like this:

So, now we have to detect our collisions as we did before, but we can't take any shortcuts this time. In order to detect if our particles have penetrated the planes we must perform a dot product and add the plane's distance to origin:

for all particles i
{
  for all planes j
  {
    distance = P[i].x*Planes[j].a + P[i].y*Planes[j].b + Planes[j].c;

     if (distance < 0)
    {
      // collision responce
    }
  }
}

Now, if we take a closer look at that distance check above, including each plane's coefficents

plane0dist = P[i].x*1  + P[i].y*0  + 1;
plane1dist = P[i].x*0  + P[i].y*-1 + 1;
plane2dist = P[i].x*-1 + P[i].y*0  + 1;
plane3dist = P[i].x*0  + P[i].y*1  + 1;

If we simplify that down a bit:

plane0dist =   P[i].x + 1;
plane1dist =  -P[i].y + 1;
plane2dist =   -P[i].x + 1;
plane3dist =   P[i].y + 1;

Re-arranging slightly we get these plane tests conditions:

if (P[i].x  < -1)
if (-P[i].y < -1) = if (P[i].y > 1)
if (-P[i].x < -1) = if (P[i].x > 1)
if (P[i].y  < -1)

Which are exactly the same as those which we used in our first simple example. The second condition for each plane must also be satisfied so we only do the collision response once per collision. This can be done by performing a dot product between the particle's velocity and the normal for each plane;

if (P[i].x < -1 && V[i]•N[i] < 0)

We call this the normal velocity since its the particle velocity in the direction of the normal. If this results in a value less than zero then the particle is moving towards the plane and we allow the collision. Once the collision is resolved, the normal velocity will be >= 0 depending on the coefficient of restitution. I could perform the same analysis on the second condition as i did the first to prove that this general solution is the same as the first specific version, but i will leave that as an excercise for the reader.

This shows our two approaches are physically the same.

Ok, so what do we do now for the collision response?

We need a solution which achieves the same result as the original program but is more general. We can do using the reflection vector from the law of reflection. The reflection vector is calculated thus:

R = V - 2*N*(V•N)

Where V is the velocity vector and N the surface normal. We need to do a similar comparison with this operator as we did with the collision checks:

Which, you should be able to see, are the exact same resulting velocities after collision with each plane in the first simple example.

You will have noticed that mass has been ignored in our calculations so far, this is because it just drops right out of the equations when you are dealing with collision against a body of infinite mass (i.e. our simple world).

So what does this new version look like so far in pseudo-code?

Move over applet to start simulation

Ok, so this is all very well but its not a very real world simulation. There is no gravity and the coefficient of restistution is far too high - in the real world, most everyday things have a near plastic coefficient of restitution, i.e near 0.

So, adding gravity is pretty simple, we just subtract gravity from our particle's velocities before we do any collision, remembering to account for the fact that gravity is an acceleration:

V[i] += G*f

In this case G is the vector {0, -9.8, 0} and f is the number of seconds per frame as before.

We would also like to use a different coefficient of restitution as well because our particles might not be made of super-rubber. Recall the reflection vector equation that we used earlier:

R = V - 2*N*(V•N)

We can actually re-write this equation to include the coefficient of restitution:

R = V - (1+e)*N*(V•N)

Where e is the coefficient of restitution which varies from 0 (totally plastic) to 1 (totally elastic). You can see that these two equations are equivelent; if we set e to 1, we get the original equation. The particles in this applet have coefficients of restitution varying from 0 to 1.

Click left mouse over applet to move particles

The only thing we can now do to this particle simulation to make it more realistic (without going to a full rigid body simulation with rotation) is to handle masses in our collisions. To accomplish this, we will have to use circles instead of particles, since its quite unlikely that two particles will ever collide what with them having no actual size! I'm presuming you have read through the section on collision detection with circles .

In order to deal with the masses of our particles we need to do some thinking regarding the reflection equation that we already have. What we really need is some kind of mass ratio to apply a bias to the equation so that lighter particles get pushed out of the way by heavier particles, but we need it to conserve momentum as well.

ratio_a = M_b / (M_a + M_b)
ratio_b = M_a / (M_a + M_b)

The above two ratios will do the trick. Here is an example:

M_a = 1.0
M_b = 0.5

ratio_a = 0.5 / (1.0 + 0.5) = 1/3
ratio_b = 1.0 / (1.0 + 0.5) = 2/3

Obviously 1/3 + 2/3 = 1, so you can see that this conserves momentum correctly, and since body b is half the mass of body a, it should experience twice the reaction upon collision, which it does.

So we can plug these two ratios into our reflection vector equation to obtain two new equations for the reflected velocity of each body after collision:

V_a = V_a - (1+e)*N*(V_a • N)*(M_b / (M_a+M_b))
V_b = V_b - (1+e)*-N*(V_b • -N)*(M_a / (M_a+M_b))

Our collision code gives us the normal pointing from body b towards body a, so we must reverse the normal to calculate body b's reflected velocity. You will have noticed that there are a lot of duplicated calculations in the two equations above. We would like to simplify this down a bit. What we can do is use the relative velocity of the two objects in question:

V_r = V_a - V_b

Now that we have only one velocity to deal with the equation gets simpiler.

I = (1+e)*N*(V_r • N)

V_a = -I*(M_b / (M_a+M_b))
V_b = +I*(M_a / (M_a+M_b))

What we are doing here is to treat one object as being stationary relative to the other moving object. But there is still more we can do here to simplify these equations:

I = (1+e)*N*(V_r • N)*(M_a*M_b)/( M_a+M_b)

V_a - = I /  M_a
V_b + = I / M_b

By combining the ratios we can save some more calculation and reduce the final step to a simple divide though by the mass of the objects in question to produce the final velocities for those objects. You can see this works because

(M_a*M_b)/( M_a+M_b) / M_a = M_b/( M_a+M_b)

and

(M_a*M_b)/( M_a+M_b) / M_b = M_a / (M_a+M_b)

Finally we can show that (from fractions)

(M_a+M_b)/( M_a*M_b) = 1/M_a + 1/M_b

Which means we can re-write the final linear impulse equation

I = (1+e)*N*(V_r • N) / (1/M_a + 1/M_b)

This is handy because we can just store 1/M_a and 1/M_b inside the definition of our circle rigid bodies so we don't have to re-calculate. It also means that we now have a way to represent infinitly massive objects, such as our world (by storing the inverse mass as 0).

V_a - = I * 1/M_a
V_b + = I * 1/M_b

Drag circle with left mouse

The above applet shows just one rigid body circle colliding with its environment. There is no penetration correction happening here at all, its simply the physics half step mentioned above which prevents it; You can probably also see the artifact of this which is the fact that it bounces slightly before it actually hits the surface.

Already here we are actually dealing with a multi-contact type situation which requires some kind of specialised 'solver' to compute the correct impulses to resolve collision without errors; The circle can come into contact with two planes in the world at once.

Problems

I'm going to carry on here in a logical chain of progression and just let the problems happen so i can show the things a first time physics programmer will have to deal with. The next step is to increase the number of circles in the simulation.

Drag circles with left mouse

You should be able to see that we've run into a big problem here. With two rigid body circles inside a small box world we run straight into the main problem with rigid body contact: multiple simultaneous contact resolution.

Many papers hide this problem with collision by making the assumption that collisions happen when objects are colliding at speed. They turn low speed collisions into resting contact to be solved differently. However, it would be nice to deal with this correctly. For clarity, its worth remembering that i'm not distinguishing collision from contact (as many authors do), since i'm only considering coefficents of restitutions of 0.

Previous work; definitions

Before i talk about this, i feel i should introduce some concepts needed to understand the problem.

Constraint

Loosely, this is something which enforces a condition between two rigid bodies. So, when there is a collision we create a constraint which enforces the rule that the colliding rigid body may not move through the object it collides with but also that it do this by only pushing. We already derrived the code for this above (not taking into account rotation yet):

I = (1+e)*N*(V_r • N) / (1/M_a + 1/M_b)

All the constraints i'm going to work with are impulse velocity level constraints. Forces and accellerations will not come into it.

Inequality constraint

This is a constraint which only acts in one direction. So a collision constraint is an inequality constraint because it is only allowed to push and never pull. If it were allowed to pull, it would stick the rigid body to the object it collided with. The way we enforce this inequality is simply the check i mentioned at the beginning:

Equality constraint

Virtually every other type of constraint is an equality constraint. This constraint is allowed to pull and push equally. An example is a point to point constraint, where two rigid bodies are connected by a single point - one point on each body is forced to lie in the same place. With a hinge constraint the whole hinge axis is constrained to lie in the same place on both rigid bodies.

Previous work; methods for solving multiple simultaneous collisions

Now, the thing which makes solving multiple simultaneous collisions so difficult is that we are essentially trying to solve a whole system of quadratic equations and with the added constraint that we never 'pull' and only 'push' (inequality constaints). It would be easier if we were just solving a system of equality constraints and there are proven linear time methods to do this .

LCP

LCPs are a popular method for solving the multi contact problem. Originally formulated at the force/accelleration level there were problems with friction which could lead to an unfeasible system. Later reforumulated by Stuart and Trinkle at the impulse/momentum level these problems were solved. However, solving an LCP is a polynomial time problem.

Baraff uses an LCP in his paper of 1989 Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies

Quadratic programming

This is a way to find the solution to a system of simultaneous quadratic equations under a set of constraints. One example would be a set of circular rigid bodies which are constrained to stay 3 metres away from each other but allowed to fall as rigid bodies otherwise. The constraint is then the squared distance function between each body. Quadratic programming is a polynomial time problem. Baraff uses quadratic programming in his siggraph notes from 1997.

Penatly methods

Springs are inserted at points of contact to prevent interpenetration. The springs all work together to form a system of equations which 'solve' the problem. The trouble is that they can cause instability and are difficult to set up. Generally, they are frowned upon in simulation circles.

Jacobsen extended penalty methods and made them stable by using infintely stiff springs. His approach was interesting because he actually re-defined the concept of a rigid body as the set of supporting distance constraints which actually formed the rigid body itself by connected a system of particles. This results in a kind of wire mesh of edges for each body. His system was inescapably 'bouncy' because of the defomation in the wire stucture of edge constraints which occured at the point of collision and it suffered from structural instability - edge constraints could be flipped and objects turned inside out. However, it is especially easy to get something up and running quickly because of the simple concepts and the itterative nature of his system will help with the impulse based method i'm describing on this page.

Impulse methods

Mirtich originally designed the impulse based rigid body simulator for his thesis. Impulse based simulators are noted for their physical accuracy. Typically these use a series of impulses applied in sequence; The trouble being that this is not going to find the correct solution to a system of simultaneous collisions and it leads to artifacts where objects appear to jitter when they should be at rest. This problem is compounded when a stack of objects is made. The problem of simultaneous impulsive collisions was only really addressed fully in the paper by Guendelman et al Nonconvex Rigid Bodies with Stacking

They use a suite of techniques to help the sequential impulse solver reach the correct solution quickly:

• Integration half step as mentioned above which predicts collisions and prevents penetrations
• Contact graph for objects
• Contacts within an object sorted by relative normal velocity / penetration distance
• Contacts solved itteratively using restitution values from -1 to 0
• Shock propagation step

Its this paper which the method on this webpage is most fully based on.

Solving the problem

Integration half step

I've talked about this above.

Contact graph

Guendelman et al use a Directed Acyclic Graph or DAG for their objects, but i opted to go for something more simple and cheaper to construct. Each object A keeps a list of objects B which are in contact and 'above' A - this is easy to construct at collision detection time. Then at run-time this list is recursively traversed and each object is prepended to the global contact graph of all objects. The following is psudocode for building the graph:

BuildContactGraph(CRigidBody obRB)
{
   if (!obRB.m_bVisisted)
   {
      obRB.m_bVisited = true;
      for (i=0; i<obRB.m_iNumRBsAbove; i++)
      {
          BuildContactGraph(obRB.m_aobObjectsAbove[i]);
      }

      // add this rigid body to the start of the list
      gobContactGraphList.PrependNode(obRB);
   }
}
BuildGlobalContactGraph(void)
{
   // reset the contact graph
   gobContactGraphList.Reset();
   for (i=0; i<giNumRigidBodies; i++)
   {
      BuildContactGraph(gaobRB[i]);
   }
}

The resulting list contains all objects and when processed in order will result in objects getting visited from the ground upwards. It might help you to dry run this routine with a 2d triangular stack of blocks, something like this:

      0
  0 0
 0 0 0
0 0 0 0

Imagine that the pyramid is stacked with no gaps, so that each block rests on each block below it. Start from any block you like, find the touching block(s) above until there are no more (as per the recursive part) and number the block with a 1, then blocks just touching and immedaitely below (2, 3 etc) in the recursion until you run out. Then pick any other non occupied block you like and repeat until none are left, remembering that you can't visit a block once it has a number in. Then read the numbers backwards from the highest - this is the contact graph ordering. Obviously contacts against the world should be processed before any others.