SimCloth3 simulates cloth as a particle system where each vertex of the cloth mesh represents one particle. The particles are connected with forces, which depend on the currently chosen force model (see below). Each particle has a certain position and velocity. The task of SimCloth3 then is to compute the particle positions and velocities for each frame of the animation.
SimCloth3 uses an implicit Euler integration method, the details of which are described below.
Why implicit integration?
Recently, implicit integration methods have become more and more popular for dynamics simulations, since they have one major advantage over explicit methods: they are much more stable and can handle the large forces necessary for producing stiff (and thus realistic) cloth while keeping the time step large. The most often used integration method is the implicit Euler method, which is also used by SimCloth3.
Implementing an implicit integration method involves the solution of a system of non-linear differential equations for computing the positions and velocities of cloth particles.
In many cases either the velocities are solved for and the positions derived from them, or the other way round. However this leads to some problems with collision response schemes. Robust collision response requires us to have control over both particle positions and velocities in one time step. Therefore SimCloth3 solves for both the positions and velocities of particles. As it turned out, this is not much slower than computing only one of these quantities and deriving the other.
Linear vs non-linear methods for solving the equations
Often for speed reasons the non-linear system of equations is reduced to a linear one by approximating the forces acting on particles with a first-order Taylor expansion. It is then easy to use the conjugate gradients method for solving the linear system. While the resulting method is quite fast, it has some disadvantages. Firstly, it requires the computation of the derivatives of the forces acting on particles. This is easy for simple force models, but becomes increasingly difficult for more complex forces and some approximations have to be used to reduce the amount of computations. Secondly, the fact the forces are approximated reduces the realism of the simulation for large time steps and can lead to unnaturally damped or even incorrect simulation. And thirdly, using the conjugate gradient method imposes some restrictions on the forces and the resulting linear system (the matrix of the system needs to be symmetric and positive-definite). Fourthly, building the linear system is somewhat complicated and requires careful book-keeping of the interaction between particles.
It is difficult to implement all these requirements properly. Futhermore, adding more complex force models is almost impossible as the exact force derivatives become too costly to compute.
SimCloth3 avoids all these issues by solving the non-linear system directly with a non-linear conjugate gradients method. This approach is very easy to implement, stable, reasonably fast on today's hardware and produces very realistic results.
Implementation
Here are some details about the implicit integrator used by SimCloth3. Suppose that we have N particles with positions p[0], p[1], ..., p[N-1] and velocities v[0], v[1], ..., v[N-1] at the current time. We want to find the positions p'[0], p'[1], ..., p'[N-1] and velocities v'[0], v'[1], ..., v'[N-1] of the particles after a time step of dt.
The explicit Euler method is defined as:
v'[i] = v[i] + im[i] * f[i] * dt for i=0, ..., N-1 and
p'[i] = p[i] + v[i] * dt for i=0, ..., N-1
Where f[i] is the force acting on the i-th particle at the current time and im[i] is the inverse mass of the particle. Since the new positions and velocities depend only on the current state of the cloth, the implementation of this method is straightforward. Unfortunately, the method is not very stable and the cloth tends to "blow up" when large forces are used.
The implicit Euler method, on the other hand is defined as:
v'[i] = v[i] + im[i] * f'[i] * dt for i=0, ..., N-1 and
p'[i] = p[i] + v'[i] * dt for i=0, ..., N-1
where f'[i] is the force acting on the i-th particle after a time step of dt. Since we don't know p'[i], v'[i] and f'[i], we have a system of equations (non-linear in general) which needs to be solved. The non-linear method used is pretty much the same as described in [1], with some simplifying assumptions.
The algorithm takes three parameters:
- the number M of variables
- a function residual() which returns an M-dimensional vector - the direction of steepest descent for the current value of the variables
- a function add(d, f) which adds the given M-dimensional vector d multiplied by the scalar f to the current value of the variables
r=residual();
add(r, 1);
while (1) {
// Check convergence
maxt=max(r[i]*r[i], i=0..M-1);
if (maxt<precision) break;
// Choose a direction
t=sum(r[i]*r[i], i=0..M-1);
if (s==0) { d=r; oldt=t; }
else {
b=t/oldt;
oldt=t;
d[i]=r[i]+b*d[i], i=0..M-1;
}
// Compute a suitable distance along that direction
t0=sum(r[i]*d[i], i=0..M-1);
add(d, 1);
r=residual();
t1=sum(r[i]*d[i], i=0..M-1);
k=-t0/(t1-t0);
add(d, k-1);
}
Note that the algorithm does not need explicitly the variables or the system of equations that need to be solved. This is quite a useful property, since we don't need to construct the system in memory (as opposed to linear methods, which typically construct the linear matrix explicitly and thus take a lot of additional RAM).
In order to apply this algorithm to the problem of finding the positions and velocities of particles at the next moment of time, we need to determine M, an initial starting solution that will be refined, and the functions residual() and add(). Since we want to solve for both positions and velocities, M=2*N. In each M-dimensional vector x[0], x[1], ..., x[M-1] variables with even indexes represent positions and variables with odd indexes represent velocities:
x[i*2] = p'[i/2] for i=0, ..., N-1 and
x[i*2+1] = v'[i/2] for i=0, ..., N-1
The starting point is simply the positions and velocities at the current moment of time:
p'[i] = p[i] for i=0, ..., N-1 and
v'[i] = v[i] for i=0, ..., N-1
The function residual() then returns an M-dimensional vector r such that:
r[i*2] = (p[i] + v'[i]*dt) - p'[i] for i=0, ..., N-1 and
r[i*2+1] = ((v[i] + im[i]*f'[i]*dt) - v'[i])*k for i=0, ..., N-1
where k is a constant scaling factor for the velocity conditions. The introduction of this factor is necessary since in general positions and velocities have different magnitudes. The constant allows us to control the target precision of velocities relative to positions and allows the system to be solved for fewer iterations of the non-linear CG algorithm. Note that we need the forces f' for the current positions p' and velocities v' of the particles.
The function add() is straightforward, it takes one M-dimensional vector d and a scalar f, multiplies them and adds the result to the current solution:
p'[i] = p'[i]+d[i*2]*f for i=0, ..., N-1 and
v'[i] = v'[i]+d[i*2+1]*(f/k) for i=0, ..., N-1
Note that the velocity changes must be divided by the same scaling factor k as used in residual().
The algorithm is very simple to implement and requires minimal additional storage. However, the system that we are solving is non-linear and therefore convergence is not guaranteed. We solve this problem by tracking the convergence rate and reducing the time step if the solution diverges. We also adaptively determine the timestep based on the number of iterations required to solve the system, increasing the time step if there are only a few iterations and decreasing the time step if we get too many iterations. This strategy seems to produce stable results even for very large initial time steps.
Stretch forces
Linear springs
Stretch/shear models
The trispring
Bend forces
Linear springs
Face angle
Plane projection