Rigid Body Simulator

This rigid body simulator was written in late 2014 to compete in the final competition of the "Physically Based Simulation" course at ETH Zurich, where it won 1st place in the Jury Award and also received the Audience Award for "Best Project". It was used to simulate the short film "Planky with a chance of Meatballs", shown above, which was produced in collaboration with Antoine Milliez.

The main features of the simulator are robust resting contacts without jittering, friction, stable stacking and flexible rods. The main purpose of the program was to simulate large structures made from wodden planks, made popular by toys such as KAPLA or KEVA, and robustly simulate their destruction at the hands of oversized food products such as spaghetti and meatballs. It is mainly intended for offline simulation, but can achieve realtime rates for low numbers of objects. It fully exploits SIMD hardware and makes use of multithreading for performance.

Rigid body simulation can be challenging in many aspects. One of the most difficult properties to achieve is resting contacts and stable stacking. Naive approaches to simulation have great difficulties to model objects resting on a flat surface and are prone to "jittering", that is, erroneous bouncing of objects that should be still. Similarly, simulating large resting stacks of objects that do not fall apart over time can be challenging to model, in part due to the aforementioned jittering that causes stacks to slide apart.

Many different ways exist to deal with these problems. For this project, I chose to implement a constraint-based approach to rigid body physics, mainly based on the paper Iterative Dynamics with Temporal Coherence by Erin Catto. The main idea is to model objectives such as sliding contacts, friction etc. as scalar valued constraint functions that are zero only on valid simulation states. By taking a differential view and restricting rigid body movement along the gradient of these constraint functions, we can evolve the simulation while staying in a valid state (e.g. objects don't interpenetrate, objects slow down due to friction, etc.).

Without going into detail, the main problems when writing an implementation of this method is primarily in identifying full contact manifolds between all objects as well as writing a fast, robust solver for the Linear Complementarity Problem (LCP).

Computing contact manifolds is mainly an exercise in computational geometry, and generally I would advise to use existing libraries for this task. For learning purposes, I implemented it myself, restricted to certain shapes (boxes, spheres, capsules). Even so, robust manifold generation in presence of floating point round-off error can be frustratingly difficult, and much time was spent debugging corner cases. To accelerate intersection tests, I also implemented a bounding volume hierarchy (BVH).

Solving the LCP is also surprisingly difficult. On top of solving a linear system of equations, the program must also honor so-called complementarity conditions that further restrict the solution. I researched various different approaches to this problem and implemented the algorithms Conjugate Projected Gradient, Projected Steepest Descent, Projected Gauss-Seidel, Subspace Minimization with Conjugate Gradient, Minimum Residual and Symmetric LQ solvers, as well as Projected Gauss-Seidel with line search using the Armijo rule. Ultimately, Gauss-Seidel with line search performed best for my specific problem, and it was the algorithm used to produce all the videos on this page.

Even though I was mainly targeting offline simulation, performance was still a big concern during development. Unfortunately, vectorization and parallelization are not straight-forward in this case, since the LCP solver suffers from numerous dependencies between computation steps, making it prone to race conditions. To achieve full SIMD and core utilization, I had to generate a scheduling plan each frame that could schedule solver operations across SIMD registers and threads without race conditions. Ultimately I was able to reduce the problem to a more abstract edge coloring problem on a multigraph, where edges represent constraints and vertices rigid bodies. Edges with the same color can be solved in parallel, whereas edges of different color generally cannot.

In addition, the program employs contact caching to boost the performance of the solver, and sends object to sleep that haven't moved in a while, reducing their impact on the performance. Using these tricks, the solver can achieve real-time performance for low numbers of objects (at most 1000).

To be able to simulate spaghetti, I also implemented torsion-, twisting-, bending- and ball-joint constraints that allowed the program to link multiple bodies together and restrict their relative motion to model a flexible rod that wants to return to a straight resting position. By modifying the softness of the constraints, different rod stiffnesses can be simulated. The video below shows an array of rods, with stiffness increasing from the front (very flexible) to the back (relatively stiff).