Advanced Character Physics

Cloth Simulation

The fact that a stick constraint can be thought of as a really hard spring should make apparent its usefulness for cloth simulation as sketched in the beginning of this section. Assume, for example, that a hexagonal mesh of triangles describing the cloth has been constructed. For each vertex a particle is initialized and for each edge a stick constraint between the two corresponding particles is initialized (with the constraint’s “rest length” simply being the initial distance between the two vertices).

The function HandleConstraints() then uses relaxation over all constraints. The relaxation loop could be iterated several times. However, to obtain nicely looking animation, actually for most pieces of cloth only one iteration is necessary! This means that the time usage in the cloth simulation depends mostly on the N square root operations and the N divisions performed (where N denotes the number of edges in the cloth mesh). As we shall see, a clever trick makes it possible to reduce this to N divisions per frame update – this is really fast and one might argue that it probably can’t get much faster.

// Implements cloth simulation
struct Constraint {
   int particleA, particleB;
   float restlength;
};
// Assume that an array of constraints, m_constraints, exists
   void ParticleSystem::SatisfyConstraints() {
      for(int j=0; j<NUM_ITERATIONS; j++) {
         for(int i=0; i<NUM_CONSTRAINTS; i++) {
         Constraint& c = m_constraints[i];
         Vector3& x1 = m_x[c.particleA];
         Vector3& x2 = m_x[c.particleB];
         Vector3 delta = x2-x1;
         float deltalength = sqrt(delta*delta);
         float diff=(deltalength-c.restlength)/deltalength;
         x1 -= delta*0.5*diff;
         x2 += delta*0.5*diff;
      }

      // Constrain one particle of the cloth to origo
      m_x[0] = Vector3(0,0,0);
   }
}

We now discuss how to get rid of the square root operation. If the constraints are all satisfied (which they should be at least almost), we already know what the result of the square root operation in a particular constraint expression ought to be, namely the rest length r of the corresponding stick. We can use this fact to approximate the square root function. Mathematically, what we do is approximate the square root function by its 1st order Taylor-expansion at a neighborhood of the rest length r (this is equivalent to one Newton-Raphson iteration with initial guess r). After some rewriting, we obtain the following pseudo-code:

// Pseudo-code for satisfying (C2) using sqrt approximation
delta = x2-x1;
delta*=restlength*restlength/(delta*delta+restlength*restlength)-0.5;
x1 -= delta;
x2 += delta;

Notice that if the distance is already correct (that is, if |delta|=restlength), then one gets delta=(0,0,0) and no change is going to happen.

Per constraint we now use zero square roots, one division only, and the squared value restlength*restlength can even be precalculated! The usage of time consuming operations is now down to N divisions per frame (and the corresponding memory accesses) – it can’t be done much faster than that and the result even looks quite nice. Actually, in Hitman, the overall speed of the cloth simulation was limited mostly by how many triangles it was possible to push through the rendering system.

The constraints are not guaranteed to be satisfied after one iteration only, but because of the Verlet integration scheme, the system will quickly converge to the correct state over some frames. In fact, using only one iteration and approximating the square root removes the stiffness that appears otherwise when the sticks are perfectly stiff.

By placing support sticks between strategically chosen couples of vertices sharing a neighbor, the cloth algorithm can be extended to simulate plants. Again, in Hitman only one pass through the relaxation loop was enough (in fact, the low number gave the plants exactly the right amount of bending behavior).

The code and the equations covered in this section assume that all particles have identical mass. Of course, it is possible to model particles with different masses, the equations only get a little more complex.

To satisfy (C2) while respecting particle masses, use the following code:

// Pseudo-code to satisfy (C2)
delta = x2-x1;
deltalength = sqrt(delta*delta);
diff = (deltalength-restlength)
      /(deltalength*(invmass1+invmass2));
x1 -= invmass1*delta*diff;
x2 += invmass2*delta*diff;

Here invmass1 and invmass2 are the numerical inverses of the two masses. If we want a particle to be immovable, simply set invmass=0 for that particle (corresponding to an infinite mass). Of course in the above case, the square root can also be approximated for a speed-up.

Rigid Bodies

The equations governing motion of rigid bodies were discovered long before the invention of modern computers. To be able to say anything useful at that time, mathematicians needed the ability to manipulate expressions symbolically. In the theory of rigid bodies, this lead to useful notions and tools such as inertia tensors, angular momentum, torque, quaternions for representing orientations etc. However, with the current ability to process huge amounts of data numerically, it has become feasible and in some cases even advantageous to break down calculations to simpler elements when running a simulation. In the case of 3D rigid bodies, this could mean modeling a rigid body by four particles and six constraints (giving the correct amount of degrees of freedom, 4x3-6 = 6). This simplifies a lot of aspects and it’s exactly what we will do in the following.

Consider a tetrahedron and place a particle at each of the four vertices. In addition, for each of the six edges on the tetrahedron create a distance constraint like the stick constraint discussed in the previous section. This is actually enough to simulate a rigid body. The tetrahedron can be let loose inside the cube world from earlier and the Verlet integrator will let it move correctly. The function SatisfyConstraints() should take care of two things: 1) That particles are kept inside the cube (like previously), and 2) That the six distance constraints are satisfied. Again, this can be done using the relaxation approach; 3 or 4 iterations should be enough with optional square root approximation.

Now clearly, in general rigid bodies do not behave like tetrahedrons collision-wise (although they might do so kinetically). There is also another problem: Presently, collision detection between the rigid body and the world exterior is on a vertex-only basis, that is, if a vertex is found to be outside the world it is projected inside again. This works fine as long as the inside of the world is convex. If the world were non-convex then the tetrahedron and the world exterior could actually penetrate each other without any of the tetrahedron vertices being in an illegal region (see Figure 3 where the triangle represents the 2D analogue of the tetrahedron). This problem is handled in the following.

In the first case, we simply project x1 out of the invalid region like earlier (in the direction of q) and that’s it (x2 is not touched). In the second case, p is still nearest to x1 and one might reason that consequently x1 should be moved more than x2. Actually, since p=0.75 x1 + 0.25 x2, we will choose to move x1 by an amount of 0.75 each time we move x2 by an amount of 0.25. In other words, the new particle positions x1’ and x2’ are given by the expressions:

  (*)

where l is some unknown value. The new position of p after moving both particles is p’=c1 x1’+ c2 x2’.

Recall that we want p’=q, i.e., we should choose l exactly such that p’ ends up coinciding with q. Since we move the particles only in the direction of D, also p moves in the direction of D and consequently the solution to the equation p’=q can be found by solving:

  (**)

for l. Expanding the left-hand side yields:

which together with the right-hand side of (**) gives

Plugging l into (*) gives us the new positions of the particles for which p’ coincide with q.

Figure 5 shows the situation after moving the particles. We have no object penetration but now the stick length constraint has been violated. To fix this, we do yet another iteration of the relaxation loop (or several) and we’re finished.

and the new positions are then given by:

Here, we have collided a single rigid body with an immovable world. The above method generalizes to handle collisions of several rigid bodies. The collisions are processed for one pair of bodies at a time. Instead of moving only p, in this case both p and q are moved towards each other.

Again, after adjusting the particle positions such that they satisfy the non-penetration constraints, the six distance constraints that make up the rigid body should be taken care of and so on. With this method, the tetrahedron can even be imbedded inside another object that can be used instead of the tetrahedron itself to handle collisions. In Figure 6, the tetrahedron is embedded inside a cube.

First, the cube needs to be ‘fastened’ to the tetrahedron in some way. One approach would be choosing the system mass midpoint 0.25*(x1+x2+x3+x4) as the cube’s position and then derive an orientation matrix by examining the current positions of the particles. When a collision/penetration is found, the collision point p (which in this case will be placed on the cube) is then treated exactly as above and the positions of the particles are updated accordingly. As an optimization, it is possible to precompute the values of c1-c4 for all vertices of the cube. If the penetration point p is a vertex, the values for c1-c4 can be looked up and used directly. Otherwise, p lies on the interior of a surface triangle or one of its edges and the values of c1-c4 can then be interpolated from the precomputed values of the corresponding triangle vertices.

By rearranging the positions of the particles that make up the tetrahedron, the physical properties can be changed accordingly (mathematically, the inertia tensor changes as the positions and masses of the particles are changed).

Other arrangements of particles and constraints than a tetrahedron are possible such as placing the particles in the pattern of a coordinate system basis, i.e. at (0,0,0), (1,0,0), (0,1,0), (0,0,1). Let a, b, and c be the vectors from particle 1 to particles 2, 3, and 4, respectively. Constrain the particles’ positions by requiring vectors a, b, and c to have length 1 and the angle between each of the three pairs of vectors to be 90 degrees (the corresponding dot products should be zero). (Notice, that this again gives four particles and six constraints.)