Advanced Character Physics

• A so-called Verlet integration scheme.
• Handling collisions and penetrations by projection.
• A simple constraint solver using relaxation.
• A nice square root approximation that gives a solid speed-up.
• 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 “velocity-less” 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₁, x₂, x₃).

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 time-stepping 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 velocity-less 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 well-known demo effect for creating ripples in water uses a similar approach.) It works due to the fact that 2x-x^*=x+(x-x^*) and x-x^* 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<NUM_PARTICLES; i++) {
      Vector3& x = m_x[i];
      Vector3 temp = x;
      Vector3& oldx = m_oldx[i];
      Vector3& a = m_a[i];
      x += x-oldx+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<NUM_PARTICLES; i++) m_a[i] = m_vGravity;
}
// 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.

________________________________________________________