Position-Based Dynamics and the Projection It Iterates

The standard way to simulate a cloth is to integrate forces and then correct the result. Position-based dynamics inverts the emphasis: it treats the constraints, not the forces, as the primary object, integrates positions forward with a crude predictor, and then moves the points so the constraints hold again. The method is not a physical model in the sense that a force-based solver is, but it is stable where force solvers are not, and its tuning knob, the number of projection iterations, is a stiffness in disguise. Understanding it means understanding what is being given up for that stability.

Predict, then correct, then correct again

The loop is three stages per time step. First the positions are updated with a cheap explicit step, typically semi-implicit Euler, using only the external forces such as gravity. Second, each constraint is processed in turn: a distance constraint moves its two particles along the line between them until the distance is correct; a collision constraint pushes a particle back out of a surface. Third, the velocities are reconstructed from the position change, by dividing the displacement by the step, and the loop repeats.

The reconstruction of velocity from displacement is the method’s signature. There is no stored velocity in the constraint projection; the displacement that satisfied a constraint is converted back into a velocity so the next step’s predictor has something to integrate. That conversion loses information about how the velocity was changing, and it is the source of the method’s main physical artifact, a damping that depends on how many iterations the projection ran.

/* One position-based distance constraint on particles a and b.
   w is inverse mass; 0 inverse mass means the particle is fixed.
   Moves both points toward the target length by the share of their mass. */
static void project_distance(float *a, float *b, float wa, float wb,
                             float rest, float *dx) {
    float ab[3] = { b[0]-a[0], b[1]-a[1], b[2]-a[2] };
    float len = sqrtf(ab[0]*ab[0] + ab[1]*ab[1] + ab[2]*ab[2]);
    if (len < 1e-9f) len = 1e-9f;             /* avoid a zero direction */
    float err = len - rest;                   /* how far the constraint is violated */
    float n[3] = { ab[0]/len, ab[1]/len, ab[2]/len };
    float wsum = wa + wb;
    if (wsum < 1e-9f) return;                 /* both fixed: nothing to move */
    float s = err / wsum;
    a[0] -= wa * s * n[0]; a[1] -= wa * s * n[1]; a[2] -= wa * s * n[2];
    b[0] += wb * s * n[0]; b[1] += wb * s * n[1]; b[2] += wb * s * n[2];
    dx[0] = a[0]; dx[1] = a[1]; dx[2] = a[2];  /* caller accumulates for velocity */
}

The projection moves each particle in proportion to its inverse mass, so a fixed particle with zero inverse mass does not move and a heavy particle moves less than a light one. That is the same mass weighting a force solver applies, but here it is applied as a positional share, and it is what keeps the total momentum of the pair unchanged under the correction. The share is the physical content of the projection, and everything else in the method is bookkeeping.

Constraints are equations, not forces

A distance constraint is the equation that the distance between two points equals a rest length. A bend constraint fixes the angle at a hinge between two cloth triangles. A volume constraint fixes the signed volume of a tetrahedron, which keeps a soft body from collapsing. Each is a scalar equation in the particle positions, and the solver’s job is to move the positions to satisfy all of them, which is a system of as many equations as constraints.

Force-based methods turn each constraint into a spring and solve the spring system, which is why they hit the stability limit of the stiffest spring. Position-based dynamics never forms a spring: the constraint is enforced directly, and the “stiffness” of the constraint is the number of times the projection iterates over it in a step. One iteration of a distance constraint reduces its error to zero for that pair, but the move may have broken a neighbouring constraint, which is why the solver sweeps all constraints several times. The number of sweeps is the stiffness knob, and raising it makes the cloth stiffer without any change to a force constant, because there is no force constant to change.

The trade is that the stiffness is only approximately enforced. After a finite number of sweeps the constraints hold only approximately, and the residual error is visible as stretch in the cloth, sag in a rope, or a soft body that does not keep its volume. Running more iterations reduces the residual but does not remove it, and the cost is linear in the iteration count, which is the honest price of the stiffness.

Why Gauss-Seidel order matters

Constraints are processed one at a time, each using the positions the previous constraint left behind. That is Gauss-Seidel iteration, and it has two consequences. First, the order of constraints changes the result, because each constraint sees a different state; a solver that processes constraints in a different order produces a different pose for the same input, which is a reproducibility problem across machines and compilers. Second, the Gauss-Seidel order converges faster than a Jacobi pass that reads a snapshot and writes all corrections at once, because each constraint benefits immediately from its predecessor’s work.

The reproducibility problem is real and underappreciated. Two builds that iterate constraints in a different order, or that parallelize the sweep and merge corrections, produce poses that differ by the residual error. The differences are small but they do not cancel, and a deterministic replay that requires bit-identical results must fix the order and the parallelization scheme. The standard fix is a graph-colouring: constraints that do not share a particle are independent and can be processed in any order, so a solver partitions the constraints into colour groups and processes the groups in a fixed order, which restores determinism while keeping most of the parallelism.

Colouring is also where the method’s parallel limits show. The number of independent constraints per sweep is bounded by the graph structure, and a highly connected constraint graph, a cloth with many bend constraints sharing edges, has few colours and little parallelism. The sequential Gauss-Seidel sweep is the fallback, and it is the reason a position-based solver on a many-core machine does not scale linearly with cores: the constraint graph, not the machine, limits the parallelism.

Stiffness as a continuum from force to constraint

Position-based dynamics and force-based methods are not competitors so much as two ends of a knob. At one end a very soft spring integrated explicitly behaves like a force-based cloth; at the other end the same constraint projected many times behaves like a rigid constraint. The middle, a spring projected a few times, is the useful region for cloth, hair, and soft bodies, and it is exactly the region where neither a pure force solver nor a hard constraint solver works well.

The continuum is why the method’s damping is hard to pin down. A projection that moves a particle back to a constraint removes energy, and the amount removed depends on how far the predictor overshot, which depends on the step, the stiffness of the neighbours, and the iteration count. The effective damping is not a constant the author sets; it emerges from those three. Authors who need a specific bounce set the iteration count and step to hit a feel, then leave them fixed, because changing either changes the feel in a way that is hard to predict from the equations.

Restitution and friction are added as modifications to the projection. A restitution coefficient scales how much of the displacement is kept after a collision, and a friction coefficient scales a tangential displacement; both are applied inside the projection and inherit its non-physical status. The result is controllable and visually convincing, and it is not a rigorous contact model, which is the same bargain the rest of the method makes.

Where the method is the wrong tool

Position-based dynamics gives up the things a force solver guarantees. It does not conserve energy, it does not conserve angular momentum except in special cases, and its time step is coupled to its stiffness in a way that makes the physical properties depend on the frame rate. A cloth simulated at 60 Hz with a given iteration count is stiffer than the same cloth at 30 Hz, because the projection runs more often per second; that frame-rate dependence is absent in a properly solved implicit method.

The loss of angular momentum matters for anything that must spin. A soft body or a rope that rotates does not preserve its spin exactly, and over many seconds the spin decays or grows depending on the projection order, which reads as a slowly changing tumble. Simulations where rotation fidelity matters, gyroscopes, spinning wheels, any rigid body with real rotational dynamics, are better served by a rigid-body solver with a proper angular velocity integrator than by position-based projection.

The method also does not give a clean way to express some constraints that forces express naturally. A constraint that a body stays at a constant temperature, or that a fluid stays incompressible, has a positional form that is awkward to project and loses information each step. Position-based dynamics is a good fit for the constraints that are naturally positional, distances, contacts, volumes, and a poor fit for those that are naturally energetic.

For the cases it fits, it is the workhorse. The cloth and rope in a real-time character are almost always position-based, because they must never explode, must look acceptable under any user input, and must run at a fixed cost, and those are the properties the method guarantees by construction. The sequential-impulse contact solver that keeps a character on the ground is the same idea applied to rigid bodies, and the two are siblings in the same family of projection methods.