Position Based Dynamics (PBD): Cloth, Soft Bodies & Constraints

1. Why Not Just Use Forces?

Traditional force-based simulation (mass-spring systems, FEM) computes forces at the current state, then integrates Newton's second law to update positions and velocities. The problem: stiff springs require tiny timesteps to stay stable, and cloth has very stiff springs (otherwise it stretches unrealistically).

A cloth panel with 10 000 particles connected by springs with stiffness k = 10⁶ N/m requires a timestep Δt < √(m/k) ≈ 10⁻⁴ s to remain stable with explicit integration — 600 substeps per 60 fps frame. Implicit integration would require solving a large sparse linear system each substep.

The key insight of PBD: geometry constraints (distances, volumes, angles) are always satisfiable regardless of timestep size. Instead of modelling springs as forces, enforce the geometric constraint directly by moving particles to the nearest configuration that satisfies it. No stiffness constant → no stability condition → any timestep works.

PBD was introduced by Müller et al. in their 2007 SIGGRAPH paper "Position Based Dynamics" and has since become the dominant method for real-time cloth, hair, and soft-body simulation in games and film visual effects toolchains.

2. The PBD Algorithm Loop

The complete PBD simulation loop per frame is:

// PBD main loop (Müller et al. 2007) for each particle i: v_i += Δt \cdot (f_ext / m_i) // apply external forces (gravity, wind) p_i = x_i + Δt \cdot v_i // predict new position generate collision constraints(x_i, p_i) // detect overlaps for iter = 1 to solverIterations: // constraint projection loop for each constraint C_j: projectConstraint(C_j, p_1...p_n) // move predicted positions for each particle i: v_i = (p_i - x_i) / Δt // derive velocity from position change x_i = p_i // commit positions apply velocity damping apply friction / restitution from collision constraints

The crucial property: the inner loop only modifies the predicted positions p, not the actual positions x. Velocity is recovered at the end as a finite difference. This means unstable oscillations cannot build up — any over-correction in p is absorbed as a velocity change at the end of the frame, not fed back into the next constraint projection.

The number of solverIterations controls the quality-performance tradeoff: 1 iteration is fast but constraints are not fully satisfied; 10-20 iterations give near-rigid behaviour for structural constraints.

3. Constraint Types

A constraint is a function C(p₁, …, pₙ) = 0 (equality) or ≥ 0 (inequality) that the predicted positions must satisfy. Common constraint types for cloth and soft bodies:

Distance Constraint

Enforces a fixed rest length between two particles — the cloth "edge" constraint:

C_dist(p₁, p₂) = |p₁ - p₂| - d_rest = 0 d_rest = initial edge length (rest configuration) Violated when: current distance \neq d_rest

Bending Constraint

Resists bending along an edge shared by two triangles. The dihedral angle between the two triangles should match the rest angle:

C_bend(p₁, p₂, p₃, p₄) = arccos(n₁ \cdot n₂) - φ_rest = 0 n₁, n₂ = unit normals of the two triangles sharing the edge (p₃, p₄) φ_rest = dihedral angle at rest Alternatively: use the quadratic bending constraint formulation C = (p₁ - p₄)\cdot(p₂ - p₄) - cos(φ_rest)\cdot|\dots|\cdot|\dots| (avoids arccos, cheaper to differentiate)

Volume Conservation Constraint

For soft bodies (deformable meshes), the total volume of a tetrahedron should be preserved:

C_vol(p₁, p₂, p₃, p₄) = (1/6) \cdot (p₂-p₁)\cdot((p₃-p₁)\times(p₄-p₁)) - V_rest = 0 V_rest = rest volume of the tetrahedron

Collision Constraints

Generated dynamically when a particle penetrates a surface. The inequality constraint pushes the particle back to the surface:

C_col(p) = (p - q) \cdot n \geq 0 q = closest point on the collision surface n = surface normal at q Satisfied when: particle is on or above the surface

Shape-Matching Constraint

For rigid or near-rigid bodies: compute the optimal rigid transformation (rotation + translation) that maps rest positions to current positions, then pull all particles toward their transformed rest positions. This preserves the overall shape while allowing some elastic deformation.

4. Constraint Projection: The Core Maths

Given a violated constraint C(p₁,…,pₙ) ≠ 0, how do we minimally move the particles to satisfy it? PBD solves a constrained minimisation: find corrections Δpᵢ that minimise ½Σ(1/wᵢ)|Δpᵢ|² subject to C = 0, where wᵢ = 1/mᵢ is the inverse mass.

// General constraint projection formula Correction for particle i: Δpᵢ = -wᵢ \cdot [C(p₁\dotspₙ) / Σⱼ wⱼ |\nablapⱼ C|²] \cdot \nablapᵢ C Where: wᵢ = 1/mᵢ (inverse mass — zero for pinned/static particles) \nablapᵢ C = gradient of constraint function w.r.t. pᵢ C(p) = current constraint value (violation amount) Intuition: Move each particle along the gradient of C, scaled by its inverse mass (lighter particles move more) and by the violation amount (small C \to small correction)

Distance Constraint Projection (Closed Form)

// For C(p₁, p₂) = |p₁ - p₂| - d_rest = 0 n = (p₁ - p₂) / |p₁ - p₂| // unit direction Δp₁ = -w₁/(w₁+w₂) \cdot (|p₁-p₂| - d_rest) \cdot n Δp₂ = +w₂/(w₁+w₂) \cdot (|p₁-p₂| - d_rest) \cdot n p₁ += Δp₁ p₂ += Δp₂ If both particles have equal mass (w₁=w₂=1/m): Each moves half the violation distance toward the other. They meet exactly at the midpoint of the over/under-stretched segment.

This closed-form solution is why distance constraint projection is so cheap — just a few float operations per constraint, no matrix factorisation, no iterative solver. The same simplicity applies to almost all PBD constraint types: the gradient ∇C is analytic, and the projection formula has a closed-form solution.

Gauss-Seidel ordering: Constraints are projected sequentially (not simultaneously), and later constraints see the corrections made by earlier ones — this is a Gauss-Seidel iteration scheme. The order matters for convergence; random shuffling each iteration step helps avoid systematic bias and improves convergence for overdetermined systems.

5. Cloth Simulation with PBD

A cloth mesh is a regular or irregular grid of particles connected by three sets of constraints:

// Cloth constraint structure for a grid mesh Structural constraints: edges between adjacent particles (resist stretching/compression) \leftarrow most stiffness, most iterations needed Shear constraints: diagonal edges across each quad (resist shearing deformation — cloth "wrinkles") Bending constraints: edges connecting next-nearest neighbours (resist bending around each edge between two triangles) \leftarrow smallest stiffness, optional for thin fabrics Typical particle count: 400 (20\times20) for game cloth 10000 (100\times100) for film Pinned particles: set wᵢ = 0 (infinite mass) \to fixed attachment points (collar, belt loops)

Wind and self-collision are handled as additional forces and constraints respectively. Self-collision is the most expensive part: a spatial hash or BVH finds nearby particle pairs, and a distance inequality constraint prevents interpenetration.

Cloth Hanging from Two Pins: Minimal Implementation

// JavaScript cloth skeleton (~60 lines core) const N = 20; // grid size NxN const particles = []; const constraints = []; // Create particles for (let y = 0; y < N; y++) for (let x = 0; x < N; x++) { particles.push({ pos: [x*0.05, -y*0.05, 0], prev: [x*0.05, -y*0.05, 0], invMass: (y===0 && (x===0 || x===N-1)) ? 0 : 1 // pin top corners }); } // Add distance constraints for structural edges for (let y = 0; y < N; y++) for (let x = 0; x < N; x++) { if (x+1 < N) addDistConstraint(idx(x,y), idx(x+1,y)); if (y+1 < N) addDistConstraint(idx(x,y), idx(x,y+1)); } function step(dt) { // 1. Predict (Verlet-style: v implicit in pos-prev difference) for (const p of particles) { if (p.invMass === 0) continue; const v = sub(p.pos, p.prev); p.prev = [...p.pos]; p.pos = add(p.pos, add(v, scale([0, -9.8*dt*dt, 0], 1))); // gravity } // 2. Project constraints (10 iterations) for (let i = 0; i < 10; i++) for (const c of constraints) projectDistance(c); // velocity is implicitly (pos - prev)/dt — no explicit velocity storage }

Note the Verlet-style integration: velocity is implicit in the difference pos - prev. This is the "XPBD" or "PBD with Verlet" variant commonly used in practice — it avoids storing explicit velocities and gives automatic Rayleigh damping proportional to velocity.

6. Stiffness, Iterations, and Convergence

In classic PBD, stiffness is controlled only by the number of solver iterations and a per-constraint stiffness parameter k ∈ [0, 1] that partially scales the correction:

// Stiffness scaling in classic PBD Δpᵢ_effective = k \cdot Δpᵢ_projected k = 1.0 \to full correction each iteration (stiff / near-rigid) k = 0.3 \to partial correction (soft / elastic) Problem: k is iteration-count-dependent! With 1 iteration: k=0.5 gives a certain effective stiffness With 10 iterations: k=0.5 gives MUCH stiffer behaviour (because corrections accumulate) \to Stiffness in classic PBD is not timestep- or iteration-independent. A cloth tuned for 10 iterations badly misbehaves at 5 iterations.

This iteration-count dependency is the central practical problem of classic PBD. Doubling iterations makes cloth stiffer; halving iterations (due to a performance budget change) makes cloth floppier. Artists re-tune the simulation every time the solver budget changes — a painful workflow.

Convergence is also slow for stiff constraints in large meshes — the Gauss-Seidel-like projection propagates constraint corrections at the speed of O(1) constraint per iteration, so a wave of deformation needs O(N) iterations to propagate across N particles. Techniques to mitigate this:

Hierarchical / multigrid PBD: Project constraints at multiple resolution levels simultaneously.
Chebyshev acceleration: Over-relax corrections with an iteration-dependent coefficient to accelerate convergence.
Constraint grouping: Independent constraint sets (colourable as non-adjacent graph edges) can be projected in parallel — practical for GPU PBD.

7. XPBD: Extended PBD with Compliance

XPBD (Macklin et al. 2016) fixes the iteration-count dependency by introducing a compliance parameter α that is the reciprocal of physical stiffness and modifies the projection formula with a Lagrange multiplier:

// XPBD constraint projection α̃ = α / Δt² // scaled compliance (α = 1/stiffness) // α=0 \to rigid, α\to\infty \to no constraint Lagrange multiplier update: Δλ = -(C + α̃\cdotλ) / (Σᵢ wᵢ|\nablapᵢ C|² + α̃) Particle correction: Δpᵢ = wᵢ \cdot Δλ \cdot \nablapᵢ C λ initialised to 0 at the start of each timestep. λ accumulates across iterations (not reset each iteration!) Result: - Stiffness is now controlled by α (physical units: m/N) - Behaviour is timestep-independent: same α gives same result at Δt=1/30 s or Δt=1/240 s - Behaviour is iteration-independent: more iterations \to better convergence to the same physical solution, not stiffer behaviour

XPBD is the successor to PBD used in all modern physics engines. Unreal Engine's Chaos Physics, NVIDIA PhysX 5, and Unity's ECS cloth all implement XPBD. The compliance parameter α maps directly to physical material properties: α ≈ 10⁻⁸ m/N for near-rigid bone, α ≈ 10⁻⁴ m/N for stiff fabric, α ≈ 10⁻¹ m/N for soft elastic tissue.

Positional vs. Velocity-Level Constraints

XPBD also introduces velocity-level constraint correction at the end of the solver loop (damping), so energy can be correctly dissipated without relying on artificial position-level damping that breaks the XPBD energy model. This enables physically accurate material damping as a separate α_damp parameter.

8. Applications and Industry Usage

PBD and XPBD are the dominant real-time simulation methods across games, film, and medical simulation:

Games

Cloth: Capes, banners, hair, ropes in virtually every AAA game since 2010. Often 1000–5000 particles per cloth piece, 5–15 solver iterations, substeps: 1–3 per frame.
Rigid body stacking: XPBD with high stiffness replaces impulse-based solvers for pile stability with similar quality but unified solver code.
Destruction: Fracture pieces modelled as PBD distance constraints between fragments — breaking constraint when force exceeds a threshold.

Film & VFX

Houdini Vellum: Production cloth, hair, and grain solver — entirely XPBD-based. Used on dozens of Oscar-winning films. Can run 200 000+ particle simulations with substep counts of 20–100.
Maya nCloth: Earlier PBD-derived cloth solver in Maya, widely used in TV production.

Medical & Scientific

Surgical simulators: PBD soft tissue deformation for laparoscopic training — real-time, stable even with large deformations from surgical tools.
Fluid-like granular materials: PBD density constraints simulate sand, mud, and non-Newtonian fluids. Macklin's "Position Based Fluids" (2013) adds a fluid density constraint alongside SPH kernel functions.

GPU PBD: Cloth constraints can be coloured as a graph (non-adjacent constraints share no particles and can project simultaneously). With 4-8 colour passes per iteration, an entire cloth mesh can be projected on the GPU in parallel using compute shaders. NVIDIA FleX and Chaos Cloth use this approach to simulate 100 000+ particles at 60 fps on modern GPUs.

PBD's blend of simplicity, stability, and physical plausibility — even if not strict accuracy — has made it the foundation of real-time physics simulation for over 15 years. XPBD's compliance model now bridges the gap toward physically accurate behaviour, making the method competitive with FEM for many applications where millisecond-scale simulation budgets apply.