In this project we build a simple physically-based real-time simulation for rigid bodies. Its been a long while since I last worked on anything physics-related and I thought it is time to revisit the subject and restore some of my former self-confidence for this type of work.

Aside from the basics (i.e. position/velocity computations and collision detection/resolution), we are also interested in making our simulation decently scalable to higher quantities of entities, say, up to a few thousand. To do that, we will be using some spatial partitioning techniques.

Table of contents


Rigid body basics

For simplicity's sake, rigid bodies in our simulation are disc-shaped, and so in addition to a radius \(r\), they each possess a mass \(m\), position \(\boldsymbol{x}\), velocity \(\boldsymbol{v}\), and experience an applied force \(\boldsymbol{F}\).

We advance the simulation in time steps with some reasonably small duration \(dt\), and use \(\boldsymbol{F}\) to derive new position and velocity \(\boldsymbol{x_1}\), \(\boldsymbol{v_1}\) from their previous values \(\boldsymbol{x_0}\), \(\boldsymbol{v_0}\) using the fundamental equations for motion:

$$ \boldsymbol{a} = \frac{\boldsymbol{F}}{m} $$ $$ \boldsymbol{x_1} = \boldsymbol{v_0}dt + \frac{1}{2}\boldsymbol{a}(dt)^2 $$ $$ \boldsymbol{v_1} = \boldsymbol{v_0} + \boldsymbol{a}dt $$

Handling collisions

There are three steps to the handling of collisions: (1) detection, (2) separation, and (3) resolution.


Our initial objective is to determine whether two bodies \(a = (r_1, m_1, \boldsymbol{x_1}, \boldsymbol{v_1})\) and \(b = (r_2, m_2, \boldsymbol{x_2}, \boldsymbol{v_2})\) collide in the next \(dt\) seconds. First we check if the bodies in their current positions intersect, and if so report the pair as a collision:

$$ |\boldsymbol{x_1} - \boldsymbol{x_2}| \le (r_1 + r_2) \text{ iff the bodies intersect.} $$

Otherwise, we must find out whether the two collide in the future. Let \(D(t)\) be the squared distance between the centers of \(a\) and \(b\) at time \(t\):

$$ D(t) = |(\boldsymbol{x_1} + \boldsymbol{v_1}t) - (\boldsymbol{x_2} + \boldsymbol{v_2}t)|^2 $$

We are interested in finding out whether \(D(t) = (r_1 + r_2)^2\) on the interval \([0, dt]\). Using the following shorthands:

$$ \boldsymbol{dx} = \boldsymbol{x_1} - \boldsymbol{x_2} $$ $$ \boldsymbol{dv} = \boldsymbol{v_1} - \boldsymbol{v_2} $$ $$ R = r_1 + r_2 $$

...the equation can be expressed like so:

$$ (\boldsymbol{dx} + \boldsymbol{dv}t) \cdot (\boldsymbol{dx} + \boldsymbol{dv}t) - R^2 = 0 $$

Opening up the left side yields a quadratic equation \(at^2 + bt + c = 0\). Solve it to find two roots \(t_0\) and \(t_1\). All that is left to do is check whether either of them is in \([0, dt]\). If both of them are, then the time of collision is the smaller of the two: \(t_c = \text{min}(t_0, t_1)\).


Before we proceed with computation of the effects of the collision on \(a\) and \(b\), we first take a moment to bring them to the point of collision. In the case where they are already colliding and penetrating each other, we first separate them along the axis of collision. Otherwise, if the collision happens at time \(t_c \in [0, dt]\) then we first advance both up to that point.


Now that we have a pair of colliding bodies at the time of their collision, we look to compute their new velocities immediately proceeding it. However, before we do that we must determine the factor of kinetic energy preserved, i.e. the collision's coefficient of restitution \(e\).

For our purposes, we compute \(e\) as the average of two elasticity properties \(e_1\), \(e_2\) of the colliding bodies (these are constant, just like a body's radius and mass):

$$ e = \frac{1}{2}(e_1 + e_2) $$

Using \(e\), we compute the new velocities for both bodies according to the motion equation for inelastic collisions. When \(e = 1\), the collision is said to be perfectly elastic, and when \(e = 0\), it is said to be perfectly inelastic. Figure 1 illustrates collision resolution for different values of \(e\).

Figure 1: Collision resolution for different values of the coefficient of restitution.
\(e = 1\)
\(e = 0\)

Spatial partitioning

Now that we know how to determine whether two bodies collide and how to resolve such a collision, one question remains: How should we go about finding all future-colliding pairs? The naive option of course is to check all pair combinations. However, this quickly becomes unfeasible when the number of bodies in the simulation \(n\) grows. For \(n = 1000\), we already need to check \(10^6\) pairs for a single time step, and there are many time steps per second! Real-time simulation is not going to happen unless we find a way to trim the number of pairs to check.

One way to do this is to divide the simulation space into partitions, and assume that bodies in separate partitions are too far away from each other to collide in the near future. This way we need only check for colliding pairs in those partitions of space containing more than one body.

There are many types of spatial partitions, but for our purposes a simple quadtree is sufficient. The quadtree divides space into four quadrants, with each quadrant dividing its own partition into four more, and so on up to some predetermined depth. Now, we need only check for collision those groups of entities that share a leaf in the quadtree. Figure 2 illustrates this.

Figure 2: Using a quadtree to partition space allows us to check for collision only for groups of bodies that share a leaf (highlighted squares).

Source code

You can view the source code and documentation of the project on GitHub.

See also