Socalled penaltybased schemes handle contact by inserting springs at the penetration points. While this is very simple to implement, it has a number of serious drawbacks. For instance, it is hard to choose suitable spring constants such that, on one hand, objects don’t penetrate too much and, on the other hand, the resulting system doesn’t get unstable. In other schemes for simulating physics, collisions are handled by rewinding time (by binary search for instance) to the exact point of collision, handling the collision analytically from there and then restarting the simulation – this is not very practical from a realtime point of view since the code could potentially run very slowly when there are a lot of collisions.
Here, we use yet another strategy. Offending points are simply projected out of the obstacle. By projection, loosely speaking, we mean moving the point as little as possible until it is free of the obstacle. Normally, this means moving the point perpendicularly out towards the collision surface.
Let’s examine an example. Assume that our world is the inside of the cube (0,0,0)(1000,1000,1000) and assume also that the particles’ restitution coefficient is zero (that is, particles do not bounce off surfaces when colliding). To keep all positions inside the valid interval, the corresponding projection code would be:
// Implements particles in a box
void ParticleSystem::SatisfyConstraints() {
for(int i=0; i
Vector3& x = m_x[i];
x = vmin(vmax(x, Vector3(0,0,0)),
Vector3(1000,1000,1000));
}
}
(vmax operates on vectors taking the componentwise maximum whereas vmin takes the componentwise minimum.) This keeps all particle positions inside the cube and handles both collisions and resting contact. The beauty of the Verlet integration scheme is that the corresponding changes in velocity will be handled automatically. In the following calls to TimeStep(), the velocity is automatically regulated to contain no component in the normal direction of the surface (corresponding to a restitution coefficient of zero). See Figure 1.


Figure 1: Ten timesteps and two particles. 
Try it out – there is no need to directly cancel the velocity in the normal direction. While the above might seem somewhat trivial when looking at particles, the strength of the Verlet integration scheme is now beginning to shine through and should really become apparent when introducing constraints and coupled rigid bodies in a moment.
A common model for cloth consists of a simple system of interconnected springs and particles. However, it is not always trivial to solve the corresponding system of differential equations. It suffers from some of the same problems as the penaltybased systems: Strong springs leads to stiff systems of equations that lead to instability if only simple integration techniques are used, or at least bad performance – which leads to pain. Conversely, weak springs lead to elastically looking cloth.
However, an interesting thing happens if we let the stiffness of the springs go to infinity: The system suddenly becomes solvable in a stable way with a very simple and fast approach. But before we continue talking about cloth, let’s revisit the previous example. The cube considered above can be thought of as a collection of unilateral (inequality) constraints (one for each side of the cube) on the particle positions that should be satisfied at all times:
(C1)
In the example, constraints were satisfied (that is, particles are kept inside the cube) by simply modifying offending positions by projecting the particles onto the cube surface. To satisfy (C1), we use the following pseudocode
// Pseudocode to satisfy (C1)
for i=1,2,3
set xi=min{max{xi, 0}, 1000}
One may think of this process as inserting infinitely stiff springs between the particle and the penetration surface – springs that are exactly so strong and suitably damped that instantly they will attain their rest length zero.
We now extend the experiment to model a stick of length 100. We do this by setting up two individual particles (with positions x1 and x2) and then require them to be a distance of 100 apart. Expressed mathematically, we get the following bilateral (equality) constraint:
Although the particles might be correctly placed initially, after one integration step the separation distance between them might have become invalid. In order to obtain the correct distance once again, we move the particles by projecting them onto the set of solutions described by (C2). This is done by pushing the particles directly away from each other or by pulling them closer together (depending on whether the erroneous distance is too small or too large). See Figure 2.


Figure 2: Fixing an invalid distance by moving particles. 
The pseudocode for satisfying the constraint (C2) is
// Pseudocode to satisfy (C2)
delta = x2x1;
deltalength = sqrt(delta*delta);
diff = (deltalengthrestlength)/deltalength;
x1 = delta*0.5*diff;
x2 += delta*0.5*diff;
Note that delta is a vector so delta*delta is actually a dot product. With restlength=100 the above pseudocode will push apart or pull together the particles such that they once more attain the correct distance of 100 between them. Again we may think of the situation as if a very stiff spring with rest length 100 has been inserted between the particles such that they are instantly placed correctly.
Now assume that we still want the particles to satisfy the cube constraints. By satisfying the stick constraint, however, we may have invalidated one or more of the cube constraints by pushing a particle out of the cube. This situation can be remedied by immediately projecting the offending particle position back onto the cube surface once more – but then we end up invalidating the stick constraint again.
Really, what we should do is solve for all constraints at once, both (C1) and (C2). This would be a matter of solving a system of equations. However, we choose to proceed indirectly by local iteration. We simply repeat the two pieces of pseudocode a number of times after each other in the hope that the result is useful. This yields the following code:
// Implements simulation of a stick in a box
void ParticleSystem::SatisfyConstraints() {
for(int j=0; j
// First satisfy (C1)
for(int i=0; i
Vector3& x = m_x[i];
x = vmin(vmax(x, Vector3(0,0,0)),
Vector3(1000,1000,1000));
}
// Then satisfy (C2)
Vector3& x1 = m_x[0];
Vector3& x2 = m_x[1];
Vector3 delta = x2x1;
float deltalength = sqrt(delta*delta);
float diff = (deltalengthrestlength)/deltalength;
x1 = delta*0.5*diff;
x2 += delta*0.5*diff;
}
}
(Initialization of the two particles has been omitted.) While this approach of pure repetition might appear somewhat naïve, it turns out that it actually converges to the solution that we are looking for! The method is called relaxation (or Jacobi or GaussSeidel iteration depending on how you do it exactly, see [Press]). It works by consecutively satisfying various local constraints and then repeating; if the conditions are right, this will converge to a global configuration that satisfies all constraints at the same time. It is useful in many other situations where several interdependent constraints have to be satisfied at the same time.
The number of necessary iterations varies depending on the physical system simulated and the amount of motion. It can be made adaptive by measuring the change from last iteration. If we stop the iterations early, the result might not end up being quite valid but because of the Verlet scheme, in next frame it will probably be better, next frame even more so etc. This means that stopping early will not ruin everything although the resulting animation might appear somewhat sloppier.