Iterative n-body Collision Problem

An essential part of a complete rigid-body physics engine

notes Math physics school:cal-poly=csc-572 project:planet-collision-572

For my CSC 572 final project, I implemented a 3D rigid-body physics engine. Naturally, the problem of n-body collisions arises and we need to figure out what to do about them. Here are notes I wrote about my implementation.


Particles are spheres. They essentially contain the following information:

  • mass mm
  • radius rr
  • moment of inertia
  • position vector ss
  • velocity vector vv
  • rotation matrix RR
  • angular velocity vector ω\omega (right-hand rule, scaled to angular velocity in radians)

Contacts essentially contain the following information:

  • a reference to the two particles, A and B. A and B both contain references to the contact, as well.
  • the position where it happened PP
  • the contact normal nn , which is a vector from B to A, with length equal to the penetration depth
  • a state value that says if this contact is stable or not

Every loop, we calculate our list of contacts. This, along with gravity calculation, are the most computationally expensive operations simply due to it being O(n2)O(n^2) operations, so they are offloaded onto the GPU.

Calculating our list of contacts

This is a O(n2)O(n^2) operation that just checks if two particles are intersecting. If they are, it pushes a new contact onto the list contacts. We are currently working on parallelizing it on the GPU.

n-body contact solving

Once we have our list of contacts contacts, to take care of the n-body case, we essentially do the following in pseudo code:

repeat contacts.size() times:
    for contact in contacts:

where contact.solve_momentum() takes care of the 2-body case described here.

Solving contacts multiple times takes care of cases like the Newton’s Cradle. This may seem inefficient in theory because if there’s nn particles, there might be n3n 3 calls for solve_momentum(). However, in practice, most particles spend most of their time floating around freely so momentum calculations can be safely done on the CPU, and this step ends up only taking <10% of total processing time.


Suppose you have a bunch of balls in a line like this:

Velocity  |  > ....
Balls     |  o oooo

where >, <, and . are right, left, and zero respectively. We expect this situation to end up like this:

.... >
oooo o

So, by running it iteratively, we can propagate the momentum through the chain:

   > ....
   o oooo

1.  >....

2.  .>...

3.  ..>..

4.  ...>.

5.  ....>

    .... >
    oooo o

The number of iterations could be tuned, or calculated in a smarter way (I was thinking, maybe calculate contacts within contiguous groups that are all touching each other to reduce the number of iterations) but it works well enough for our purposes.

2-body problem


If two balls are leaving each other, or if they have low relative velocity/the contact is marked as “stable,” then there is no momentum or friction to be applied.

The 3D momentum problem gets reduced to a 1D momentum problem along the collision normal nn . So, we use vnv \cdot n as our velocity and plug it into the conservation of momentum equation with restitution derived here.


From the final momenta, we derive the impulse II . Multiplying II by some arbitrary coefficient approximates the normal force FnF_n , giving us the magnitude of friction f=μFn|f| = \mu F_n where μ\mu is the coefficient of friction.

Next, we need to calculate the relative surface velocity of the particles. Suppose we were only looking at a single particle. The contact point is essentially undergoing two movements:

  • it is moving at vv , the particle’s translational velocity
  • it is moving along the particle’s surface due to rotation, or ω×(Ps)\omega \times (P - s)

So, the surface velocity at the contact point of a single particle is v+ω×(Ps)v + \omega \times (P - s) . The relative surface velocity Δvsurf\Delta v_{surf} is just the difference between the two particles’ surface velocities at the contact point.

However, friction can only happen tangentially to the normal, so the direction of friction is actually in the direction of Δvsurfprojn(Δvsurf)\Delta v*{surf} - proj_n(\Delta v*{surf}) (its orthogonal component). Normalize this, and multiply by the magnitude derived earlier, and we have our final friction.