Dynamical Simulation with Finite Elements

The finite element method is a very versatile method for any kind of boundary value problem. It was first extensively used in engineering as a tool for structural analysis, but it can also be used for heat flow, fluid and electrodynamical simulation to name some fields.

The underlying principle (as the name implies) is the discretization of the domain of interest into simple connected elements like triangles, rectangles, tetrahedrons or dozens more possible types. Then the physical property (force, velocity, temperature…) is evaluated at the nodes, or vertices if you look at an isolated element, of the resulting mesh. Between vertices, for example inside the volume of a tetrahedron the variable is interpolated.

Now in structural dynamics the elements represent actual physical elements with a mass and a volume (or area in 2D). In comparison to rigid bodies, usually more than one element makes up a body. Consider for example a box:
It can be described by two triangular finite elements, or four smaller rectangular elements, or any other mesh covering the box. As a rigid body the box would be described as a … box. So why use these elements if a ‘one element’ formulation is sufficient? Because finite elements offer a number of physical properties which just cannot be described efficiently by a rigid body. For a real time simulation the most important properties are elastic behavior, permanent deformations and breakability.
But how are these elements actually formulated? To illustrate how finite elements are described consider a system of two masses connected by a spring.

A spring element with nodes 1 and 2

This is a two node one dimensional element.
The spring’s elongation is linearly proportional to the applied force:


where δ ist the elongation, k is the spring constant which describes the spring’s resiliance to deformations and u stands for the displacement from the initial positions of the masses. In equilibrium the two forces at the nodes are opposite to each other.


The forces can be expressed by the previously found relation:


And these two formulas can be written as a matrix equation:


The matrix is called the stiffness matrix K and it contains the material and geometric characteristics of the element. In this simple case it is just a 2×2 matrix but the more complicated an element becomes the size of the matrix can easily reach 12×12 or even bigger sizes.
After finding the description for one lonely spring the next step is to formulate the connection between several elements. So another spring is connected at the right node of the first spring. It has the same stiffness matrix K and subsequently a similar matrix equation, but with a new force and displacement variable for the third node.


The two stiffness matrices must be combined to form a global stiffness matrix. For this simple example it looks like this:


Where the upper index denotes the element the force stems from and the lower index the node number.
This basically describes the whole procedure behind a finite element formulation: the stiffness matrix is formulated, the mesh generated and the global stiffness matrix is computed. These steps are all done before the actual simulation in which the matrix equation is solved for u, the global displacements. This is done by calculating the inverse of K,


which depending on the size of K can be an expensive computation and generally is one reason (the other being high memory requirements) why the finite element method is rarely used for real time physics engines which should be as fast as possible.
But none the less for few elements it is still fast enough to be useful. Hence after all this introduction it is time to present an implementation. The following demo is a dynamical simulation of a bar consisting of 16 triangular elements (click on flash & press space to start the simulation):

This seems at first sight fun to look at but useless, however this finite element method should be integratable into motor2, which would give the best of both worlds. A fast physics engine for rigid bodies and the possibility to simulate deformable (and possibly breakable) objects.
Besides, this is a perfect candidate for using PixelBender since it involves heavy math crunching on big matrices, which doesn’t make sense for the current Box2D iterative impulse solver. Meanwhile Michael is working hard to release the next motor2 version soon, and once it’s done we hope to expand motor2 with finite element capabilities.

Feature-based geometrical algorithms

Performing collision queries over a large number of objects is computationally expensive and therefore it’s important to use clever algorithms to avoid bottlenecks. After all potentially colliding object pairs have been sorted out during the broad phase, the remaining task is to run some algorithms to determine if and where the objects intersect. Of course there are tons of ways to do this.

motor2 mainly uses a SAT approach for this task, which is fast, easy to implement and quite robust. The only disadvantage is that you end up with an unpleasant quadratic-time complexity when shapes come into contact. In hope to improve this situation I spent some time learning different collision detection algorithms and stumbled across some papers talking about feature-based algorithms which like the name implies operate on the features (vertex, edge, face) of two polygons/polyhedra to determine if they are disjoint.
To say in advance: It’s really hard to find a better replacement for SAT in 2D because things are just so simple in two dimensions (this behaves differently in 3D, where SAT can be really slow for complex polyhedral shapes, since you need to check many more separation axes…).

The Lin-Canny algorithm [1] seems to be the most commonly used of these types of algorithms. The V-Clip, or Voronoi-Clip algorithm [2] is an enhancement of Lin-Canny and according to the author offers superior numerical robustness and can handle penetrations and degenerate situations as well.
I’m not quite sure if the well-known GJK algorithm (Gilbert-Johnson-Keerthi) falls under this category as well – from the geometrical point of view it’s simplex-based and looks at the minkowski difference.

Anyway, while I found many Lin-Canny implementations, V-Clip is a bit rare and after adopting the algorithm to 2D (fun!) I was wondering why so few people are using it until Erin Catto, the brain behind Box2D, pointed out that the algorithm is patented – I simply missed that (always read the small print!).

Here’s the result:

The V-CLIP algorithm in 2D, shaded areas represent voronoi regions of the nearest features, the solid line the minimum positive distance

Let me shortly explain how the algorithm works: The basic theorem states that if X and Y are a pair of features on two separated shapes, then if X is contained in the Voronoi Region of Y and Y is contained in the Voronoi Region of X, X and Y are the closest features and thus must contain the closest points between those features. Knowing the closest points, we can then compute the distance between them and declare a collision when the distance falls below some small value.

In a physics simulation we are usually integrating over a small time step so we can assume that objects only move by a small amount (called spatial and temporal coherence). This also means the new pair of closest features is probably near the last one from the previous time step, so an almost constant query time is achieved when re-starting the algorithm using the last result.

Talking about performance, V-CLIP is roughly 2-3 times faster than SAT and even better – in most situations it performs independently of the number of polygon vertices – which makes sense since you are only looking at the “difference” from the last state rather than the whole picture. Unluckily, the last problem I’m facing preventing it from being useful is that you have to do some extra work in creating a good contact manifold, so at the end it’s actually slower than SAT. But I have not yet given up and have some ideas to get the job done faster.

If you are interested in the code, it’s included in motor2 since version 0.9.
I have now looked more closely at the open Lin-Canny algorithm which is very similar to the Voronoi Clip algorithm, especially in 2D where I don’t need some extra cases introduced by the paper. For now I have dumped the V-Clip classes from my engine (they weren’t used anyway) and will revise the algorithm in the future.

There also other interesting applications, for example it allows you to compute the time of impact (TOI) for conservative advancement.

[1] A Fast Algorithm for Incremental Distance Calculation
[2] V-Clip: Fast and Robust Polyhedral Collision Detection