This article explains the basic elements of an approach to physicallybased modeling which is well suited for interactive use. It is simple, fast, and quite stable, and in its basic version the method does not require knowledge of advanced mathematical subjects (although it is based on a solid mathematical foundation). It allows for simulation of both cloth; soft and rigid bodies; and even articulated or constrained bodies using both forward and inverse kinematics.
The algorithms were developed for IO Interactive’s game Hitman: Codename 47. There, among other things, the physics system was responsible for the movement of cloth, plants, rigid bodies, and for making dead human bodies fall in unique ways depending on where they were hit, fully interacting with the environment (resulting in the press oxymoron “lifelike death animations”).The article also deals with subtleties like penetration test optimization and friction handling.
The use of physicallybased modeling to produce nicelooking animation has been considered for some time and many of the existing techniques are fairly sophisticated. Different approaches have been proposed in the literature [Baraff, Mirtich, Witkin, and others] and much effort has been put into the construction of algorithms that are accurate and reliable. Actually, precise simulation methods for physics and dynamics have been known for quite some time from engineering. However, for games and interactive use, accuracy is really not the primary concern (although it’s certainly nice to have) – rather, here the important goals are believability (the programmer can cheat as much as he wants if the player still feels immersed) and speed of execution (only a certain time per frame will be allocated to the physics engine). In the case of physics simulation, the word believability also covers stability; a method is no good if objects seem to drift through obstacles or vibrate when they should be lying still, or if cloth particles tend to “blow up”.
The methods demonstrated in this paper were created in an attempt to reach these goals. The algorithms were developed and implemented by the author for use in IO Interactive’s computer game Hitman: Codename 47, and have all been integrated in IO’s inhouse game engine Glacier. The methods proved to be quite simple to implement (compared to other schemes at least) and have high performance.
The algorithm is iterative such that, from a certain point, it can be stopped at any time. This gives us a very useful time/accuracy tradeoff: If a small source of inaccuracy is accepted, the code can be allowed to run faster; this error margin can even be adjusted adaptively at runtime. In some cases, the method is as much as an order of magnitude faster than other existing methods. It also handles both collision and resting contact in the same framework and nicely copes with stacked boxes and other situations that stress a physics engine.
In overview, the success of the method comes from the right combination of several techniques that all benefit from each other:

A socalled Verlet integration scheme.

Handling collisions and penetrations by projection.

A simple constraint solver using relaxation.

A nice square root approximation that gives a solid speedup.

Modeling rigid bodies as particles with constraints.

An optimized collision engine with the ability to calculate penetration depths.
Each of the above subjects will be explained shortly. In writing this document, the author has tried to make it accessible to the widest possible audience without losing vital information necessary for implementation. This means that technical mathematical explanations and notions are kept to a minimum if not crucial to understanding the subject. The goal is demonstrating the possibility of implementing quite advanced and stable physics simulations without dealing with loads of mathematical intricacies.
The content is organized as follows. First, in Section 2, a “velocityless” representation of a particle system will be described. It has several advantages, stability most notably and the fact that constraints are simple to implement. Section 3 describes how collision handling takes place. Then, in Section 4, the particle system is extended with constraints allowing us to model cloth. Section 5 explains how to set up a suitably constrained particle system in order to emulate a rigid body. Next, in Section 6, it is demonstrated how to further extend the system to allow articulated bodies (that is, systems of interconnected rigid bodies with angular and other constraints). Section 7 contains various notes and shares some experience on implementing frictionetc. Finally, in Section 8 a brief conclusion.
In the following, bold typeface indicates vectors. Vector components are indexed by using subscript, i.e., x=(x_{1}, x_{2}, x_{3}).
Verlet integration
The heart of the simulation is a particle system. Typically, in implementations of particle systems, each particle has two main variables: Its position x and its velocity v. Then in the timestepping loop, the new position x’ and velocity v’ are often computed by applying the rules:
where Dt is the time step, and a is the acceleration computed using Newton’s law f=ma (where f is the accumulated force acting on the particle). This is simple Euler integration.
Here, however, we choose a velocityless representation and another integration scheme: Instead of storing each particle’s position and velocity, we store its current position x and its previous position x^{*}. Keeping the time step fixed, the update rule (or integration step) is then:
This is called Verlet integration (see [Verlet]) and is used intensely when simulating molecular dynamics. It is quite stable since the velocity is implicitly given and consequently it is harder for velocity and position to come out of sync. (As a side note, the wellknown demo effect for creating ripples in water uses a similar approach.) It works due to the fact that 2xx^{*}=x+(xx^{*}) and xx^{*} is an approximation of the current velocity (actually, it’s the distance traveled last time step). It is not always very accurate (energy might leave the system, i.e., dissipate) but it’s fast and stable. By lowering the value 2 to something like 1.99 a small amount of drag can also be introduced to the system.
At the end of each step, for each particle the current position x gets stored in the corresponding variable x^{*}. Note that when manipulating many particles, a useful optimization is possible by simply swapping array pointers.
The resulting code would look something like this (the Vector3 class should contain the appropriate member functions and overloaded operators for manipulation of vectors):
// Sample code for physics simulation
class ParticleSystem {
Vector3 m_x[NUM_PARTICLES]; // Current positions
Vector3 m_oldx[NUM_PARTICLES]; // Previous positions
Vector3 m_a[NUM_PARTICLES]; // Force accumulators
Vector3 m_vGravity; // Gravity
float m_fTimeStep;
public:
void TimeStep();
private:
void Verlet();
void SatisfyConstraints();
void AccumulateForces();
// (constructors, initialization etc. omitted)
};
// Verlet integration step
void ParticleSystem::Verlet() {
for(int i=0; i
Vector3& x = m_x[i];
Vector3 temp = x;
Vector3& oldx = m_oldx[i];
Vector3& a = m_a[i];
x += xoldx+a*fTimeStep*fTimeStep;
oldx = temp;
}
}
// This function should accumulate forces for each particle
void ParticleSystem::AccumulateForces()
{
// All particles are influenced by gravity
for(int i=0; i
}
// Here constraints should be satisfied
void ParticleSystem::SatisfyConstraints() {
// Ignore this function for now
}
void ParticleSystem::TimeStep() {
AccumulateForces();
Verlet();
SatisfyConstraints();
}
The above code has been written for clarity, not speed. One optimization would be using arrays of float instead of Vector3 for the state representation. This might also make it easier to implement the system on a vector processor.
This probably doesn’t sound very groundbreaking yet. However, the advantages should become clear soon when we begin to use constraints and switch to rigid bodies. It will then be demonstrated how the above integration scheme leads to increased stability and a decreased amount of computation when compared to other approaches.
Try setting a=(0,0,1), for example, and use the start condition x=(1,0,0), x^{*}=(0,0,0), then do a couple of iterations by hand and see what happens.