Mass-Spring Cloth Simulation with Verlet in GLSL

This article details the mathematical formulation and GLSL implementation of Verlet numerical integration for mass-spring cloth simulations. It covers the core integration equations, force accumulation across structural, shear, and flexion springs, and how to execute the entire simulation pipeline efficiently on the GPU using GLSL compute shaders.

The Verlet Integration Formulation

Verlet integration is a second-order numerical method that computes future particle positions based on current and previous positions, eliminating the need to store velocity explicitly. This implicit velocity representation provides natural numerical stability, making it ideal for real-time cloth systems.

The standard Verlet formula is derived from the Taylor expansion of position \(\mathbf{x}(t)\):

\[\mathbf{x}(t + \Delta t) = 2\mathbf{x}(t) - \mathbf{x}(t - \Delta t) + \mathbf{a}(t)\Delta t^2 + \mathcal{O}(\Delta t^4)\]

To account for air resistance or energy dissipation, a damping factor \(d \in [0, 1]\) is introduced to damp the implicit velocity term \((\mathbf{x}(t) - \mathbf{x}(t - \Delta t))\):

\[\mathbf{x}_{\text{next}} = \mathbf{x}_{\text{curr}} + (1 - d)(\mathbf{x}_{\text{curr}} - \mathbf{x}_{\text{prev}}) + \frac{\mathbf{F}_{\text{total}}}{m}\Delta t^2\]

Where:

\(\mathbf{x}_{\text{next}}\) is the position at \(t + \Delta t\).
\(\mathbf{x}_{\text{curr}}\) is the position at \(t\).
\(\mathbf{x}_{\text{prev}}\) is the position at \(t - \Delta t\).
\(\mathbf{F}_{\text{total}}\) is the net force acting on the particle.
\(m\) is the particle mass.
\(\Delta t\) is the simulation timestep.

Mass-Spring System Dynamics

A cloth grid consists of particles arranged in a 2D network connected by three types of springs:

Structural Springs: Connect particle \((i, j)\) to immediate orthogonal neighbors \((i \pm 1, j)\) and \((i, j \pm 1)\) to resist stretching and compression.
Shear Springs: Connect particle \((i, j)\) diagonally to \((i \pm 1, j \pm 1)\) to resist in-plane shearing.
Flexion (Bending) Springs: Connect particle \((i, j)\) to secondary orthogonal neighbors \((i \pm 2, j)\) and \((i, j \pm 2)\) to resist out-of-plane bending.

For any pair of connected particles \(\mathbf{x}_i\) and \(\mathbf{x}_j\), the restoring spring force \(\mathbf{F}_{ij}\) acting on particle \(i\) is calculated using Hooke's Law:

\[\mathbf{F}_{ij} = -k \left( \|\mathbf{x}_i - \mathbf{x}_j\| - L_0 \right) \frac{\mathbf{x}_i - \mathbf{x}_j}{\|\mathbf{x}_i - \mathbf{x}_j\|}\]

Where \(k\) is the spring stiffness coefficient and \(L_0\) is the rest length between the two particles.

The total force \(\mathbf{F}_{\text{total}}\) acting on particle \(i\) is the sum of external forces (e.g., gravity, wind) and internal spring forces:

\[\mathbf{F}_{\text{total}} = m\mathbf{g} + \mathbf{F}_{\text{ext}} + \sum_{j \in \mathcal{N}(i)} \mathbf{F}_{ij}\]

GLSL Compute Implementation

On the GPU, particle data is stored in textures or Shader Storage Buffer Objects (SSBOs). When using textures, two buffers (ping-pong) store current and previous positions.

Below is a GLSL compute shader illustrating the force computation and Verlet integration step.

#version 450 core

layout(local_size_x = 16, local_size_y = 16) in;

layout(rgba32f, binding = 0) readonly uniform image2D u_PrevPos;
layout(rgba32f, binding = 1) readonly uniform image2D u_CurrPos;
layout(rgba32f, binding = 2) writeonly uniform image2D u_NextPos;

uniform float u_DeltaTime;
uniform float u_Damping;
uniform float u_Mass;
uniform vec3  u_Gravity;
uniform ivec2 u_GridSize;

// Spring parameters: (offset_x, offset_y, rest_length, stiffness)
struct Spring {
    ivec2 offset;
    float restLength;
    float stiffness;
};

// Computes force exerted on particle at 'coord' by a neighbor at 'coord + offset'
vec3 computeSpringForce(ivec2 coord, vec3 pos, ivec2 offset, float restLength, float stiffness) {
    ivec2 neighborCoord = coord + offset;
    
    // Boundary check
    if (neighborCoord.x < 0 || neighborCoord.x >= u_GridSize.x ||
        neighborCoord.y < 0 || neighborCoord.y >= u_GridSize.y) {
        return vec3(0.0);
    }

    vec3 neighborPos = imageLoad(u_CurrPos, neighborCoord).xyz;
    vec3 delta = pos - neighborPos;
    float dist = length(delta);
    
    if (dist < 1e-6) return vec3(0.0);

    // Hooke's Law: F = -k * (x - L0) * normalize(delta)
    return -stiffness * (dist - restLength) * (delta / dist);
}

void main() {
    ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
    if (coord.x >= u_GridSize.x || coord.y >= u_GridSize.y) return;

    vec4 currData = imageLoad(u_CurrPos, coord);
    vec3 currPos = currData.xyz;
    float isPinned = currData.w; // 1.0 if fixed in space, 0.0 if dynamic

    // Pinned particles do not integrate
    if (isPinned > 0.5) {
        imageStore(u_NextPos, coord, vec4(currPos, isPinned));
        return;
    }

    vec3 prevPos = imageLoad(u_PrevPos, coord).xyz;

    // Accumulate external forces
    vec3 totalForce = u_Mass * u_Gravity;

    // Define neighbor offsets for Structural, Shear, and Bending springs
    const float d = 1.0; // Base cell dimension
    const float diag = 1.41421356 * d;

    // Structural springs
    totalForce += computeSpringForce(coord, currPos, ivec2( 1,  0), d, 1500.0);
    totalForce += computeSpringForce(coord, currPos, ivec2(-1,  0), d, 1500.0);
    totalForce += computeSpringForce(coord, currPos, ivec2( 0,  1), d, 1500.0);
    totalForce += computeSpringForce(coord, currPos, ivec2( 0, -1), d, 1500.0);

    // Shear springs
    totalForce += computeSpringForce(coord, currPos, ivec2( 1,  1), diag, 1000.0);
    totalForce += computeSpringForce(coord, currPos, ivec2(-1,  1), diag, 1000.0);
    totalForce += computeSpringForce(coord, currPos, ivec2( 1, -1), diag, 1000.0);
    totalForce += computeSpringForce(coord, currPos, ivec2(-1, -1), diag, 1000.0);

    // Bending springs
    totalForce += computeSpringForce(coord, currPos, ivec2( 2,  0), 2.0 * d, 500.0);
    totalForce += computeSpringForce(coord, currPos, ivec2(-2,  0), 2.0 * d, 500.0);
    totalForce += computeSpringForce(coord, currPos, ivec2( 0,  2), 2.0 * d, 500.0);
    totalForce += computeSpringForce(coord, currPos, ivec2( 0, -2), 2.0 * d, 500.0);

    // Acceleration
    vec3 acceleration = totalForce / u_Mass;

    // Verlet Integration step
    vec3 velocity = (currPos - prevPos) * (1.0 - u_Damping);
    vec3 nextPos = currPos + velocity + acceleration * (u_DeltaTime * u_DeltaTime);

    // Output newly integrated position
    imageStore(u_NextPos, coord, vec4(nextPos, 0.0));
}

State Cycling and Buffer Swapping

To advance the simulation each frame:

Run the compute shader with u_PrevPos containing \(\mathbf{x}(t - \Delta t)\), u_CurrPos containing \(\mathbf{x}(t)\), and u_NextPos bound to a writeable target for \(\mathbf{x}(t + \Delta t)\).
Insert a memory barrier using glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT).
Shift buffer pointers so that \(\mathbf{x}(t)\) becomes the new u_PrevPos, and \(\mathbf{x}(t + \Delta t)\) becomes the new u_CurrPos.