Cloth Simulation via Verlet Integration

1. The Mass-Spring Model

Cloth is discretised as a grid of particles (point masses) connected by springs. Three spring types capture cloth mechanics:

2. Verlet Integration

Standard Euler integration accumulates error and can become unstable for stiff springs. Verlet integration is symplectic — it conserves energy well and is unconditionally stable for small enough steps:

The key insight is that velocity is derived implicitly from the position history — there is no separate velocity variable to accumulate errors. This makes Verlet uniquely stable for constraint-based simulations.

3. Constraint Types

A distance constraint between two particles p₁ and p₂ with rest length d₀ projects both particles outward/inward to restore their distance to d₀:

Constraints are inextensibility constraints: they enforce a maximum distance (cloth doesn't stretch beyond rest length). Compression is usually allowed — cloth can wrinkle but not go negative-length.

4. Constraint Relaxation (Jakobsen Method)

With many interleaved constraints, satisfying one disturbs others. Jakobsen's approach: iterate over all constraints multiple times per frame. Each pass brings the system closer to satisfying all constraints simultaneously:

This is equivalent to Gauss-Seidel iterative solving of the constraint system. It converges quadratically near the solution and is robust to over-constrained or near-degenerate configurations that would crash matrix solvers.

5. Sphere and Floor Collisions

Collision response for a particle against a sphere of centre C and radius R: if the particle is inside the sphere, project it to the surface:

This is the elegance of Verlet: collision response is just position correction. No velocity modification needed — the integrator infers the new velocity from the corrected positions.

6. Self-Intersection

Preventing cloth from passing through itself is expensive. The naive approach (test every triangle pair) is O(n² × m²). Production solutions:

• Spatial hashing: Map particles to a hash grid; only test nearby particles and triangles within the same or adjacent cells. O(n) expected.
• Repulsion forces: Add a soft repulsion between nearby particle pairs — prevents most self-intersection without explicit resolution.
• Continuous collision detection (CCD): Test swept triangles (over the timestep) for intersection. Expensive but required for fast-moving cloth (flags in wind).
• Strain limiting (robust): Limit the per-step position change to a fraction of the particle diameter — cloth cannot teleport through itself.

For interactive demos, repulsion forces + increased constraint iterations handle 90% of visible self-intersection artefacts acceptably.

7. Full JavaScript Implementation

// Cloth simulation — Verlet integration + constraint relaxation
class Particle {
  constructor(x, y, z) {
    this.pos  = {x, y, z};
    this.prev = {x, y, z};
    this.acc  = {x: 0, y: 0, z: 0};
    this.pinned = false;
  }
  integrate(dt, damping = 0.99) {
    if (this.pinned) return;
    const {pos, prev, acc} = this;
    const vx = (pos.x - prev.x) * damping;
    const vy = (pos.y - prev.y) * damping;
    const vz = (pos.z - prev.z) * damping;
    prev.x = pos.x; prev.y = pos.y; prev.z = pos.z;
    pos.x += vx + acc.x * dt * dt;
    pos.y += vy + acc.y * dt * dt;
    pos.z += vz + acc.z * dt * dt;
    acc.x = 0; acc.y = 0; acc.z = 0;
  }
}

class Constraint {
  constructor(a, b) {
    this.a = a; this.b = b;
    const dx = b.pos.x-a.pos.x, dy = b.pos.y-a.pos.y, dz = b.pos.z-a.pos.z;
    this.rest = Math.sqrt(dx*dx + dy*dy + dz*dz);
    this.torn  = false;
  }
  satisfy(tearThreshold = Infinity) {
    if (this.torn) return;
    const {a, b} = this;
    const dx = b.pos.x - a.pos.x;
    const dy = b.pos.y - a.pos.y;
    const dz = b.pos.z - a.pos.z;
    const dist = Math.sqrt(dx*dx + dy*dy + dz*dz) || 0.0001;
    if (dist > this.rest * tearThreshold) { this.torn = true; return; }
    const diff = (dist - this.rest) / dist;
    if (!a.pinned) { a.pos.x += dx * diff * 0.5; a.pos.y += dy * diff * 0.5; a.pos.z += dz * diff * 0.5; }
    if (!b.pinned) { b.pos.x -= dx * diff * 0.5; b.pos.y -= dy * diff * 0.5; b.pos.z -= dz * diff * 0.5; }
  }
}

class Cloth {
  constructor(cols, rows, spacing) {
    this.cols = cols; this.rows = rows;
    this.particles   = [];
    this.constraints  = [];
    for (let r = 0; r < rows; r++)
      for (let c = 0; c < cols; c++) {
        const p = new Particle(c * spacing, -r * spacing, 0);
        if (r === 0) p.pinned = true; // pin top row
        this.particles.push(p);
      }
    const at = (r, c) => this.particles[r * cols + c];
    for (let r = 0; r < rows; r++)
      for (let c = 0; c < cols; c++) {
        if (c < cols-1) this.constraints.push(new Constraint(at(r,c), at(r,c+1)));  // struct H
        if (r < rows-1) this.constraints.push(new Constraint(at(r,c), at(r+1,c)));  // struct V
        if (c < cols-1 && r < rows-1) {                                         // shear
          this.constraints.push(new Constraint(at(r,c), at(r+1,c+1)));
          this.constraints.push(new Constraint(at(r,c+1), at(r+1,c)));
        }
        if (c < cols-2) this.constraints.push(new Constraint(at(r,c), at(r,c+2)));  // bend H
        if (r < rows-2) this.constraints.push(new Constraint(at(r,c), at(r+2,c)));  // bend V
      }
  }

  update(dt, gravity = -980, iters = 5, sphere = null) {
    // Apply gravity
    for (const p of this.particles) { if (!p.pinned) p.acc.y += gravity; }
    // Integrate
    for (const p of this.particles) p.integrate(dt);
    // Constraint relaxation
    for (let i = 0; i < iters; i++)
      for (const c of this.constraints) c.satisfy();
    // Sphere collision
    if (sphere) for (const p of this.particles) {
      const dx = p.pos.x-sphere.x, dy = p.pos.y-sphere.y, dz = p.pos.z-sphere.z;
      const len = Math.sqrt(dx*dx+dy*dy+dz*dz);
      if (len < sphere.r) {
        p.pos.x = sphere.x + dx/len*sphere.r;
        p.pos.y = sphere.y + dy/len*sphere.r;
        p.pos.z = sphere.z + dz/len*sphere.r;
      }
    }
  }
}

// Usage: 20×20 cloth hung from top row
const cloth = new Cloth(20, 20, 15);
const sphere = {x: 140, y: -150, z: 0, r: 80};
function loop() {
  cloth.update(0.016, -980, 5, sphere);
  // ... render via canvas or Three.js
  requestAnimationFrame(loop);
}
loop();

8. Extensions and Further Reading

• Wind forces: Apply a force perpendicular to each cloth triangle proportional to the dot product of the surface normal with a wind vector.
• Implicit integration: Baraff & Witkin (1998) — solves a sparse linear system per step, allowing 10× larger timesteps at the cost of one linear solve. Used in production VFX (Houdini, Maya Nucleus).
• Position-Based Dynamics (PBD): Müller et al. (2007) — extends the Jakobsen framework to rigid bodies, fluids and deformables in a unified projective approach. Used in NVIDIA PhysX, Unreal Chaos.
• XPBD: Extended PBD (2016) decouples constraint stiffness from iteration count, enabling compliant constraints.
• Finite element cloth: Shell elements with Saint-Venant–Kirchhoff or neo-Hookean material models. Accurate but expensive — used for offline simulation of garments.