PBD Distance Constraints in GLSL Compute Shaders

Position-Based Dynamics (PBD) simulates physical systems by directly projecting vertex positions to satisfy geometric constraints rather than integrating forces. Implementing iterative PBD distance constraints in a GLSL compute shader requires translating the mathematical projection onto the GPU's massively parallel architecture. This article details the core mathematical formulation of distance constraints, the race condition challenges inherent to parallel constraint solving, the use of graph coloring to prevent data hazards, and the complete compute shader pipeline required to iteratively converge on a stable simulation.

Mathematical Formulation of Distance Constraints

A distance constraint between two particles with positions \(\mathbf{p}_1\) and \(\mathbf{p}_2\) and a rest length \(d\) is defined by the scalar constraint function:

\[C(\mathbf{p}_1, \mathbf{p}_2) = |\mathbf{p}_1 - \mathbf{p}_2| - d\]

To satisfy \(C(\mathbf{p}_1, \mathbf{p}_2) = 0\), position corrections \(\Delta \mathbf{p}_1\) and \(\Delta \mathbf{p}_2\) are applied along the constraint gradient \(\mathbf{n} = \frac{\mathbf{p}_1 - \mathbf{p}_2}{|\mathbf{p}_1 - \mathbf{p}_2|}\), weighted by the inverse masses \(w_1 = 1/m_1\) and \(w_2 = 1/m_2\):

\[\Delta \mathbf{p}_1 = -\frac{w_1}{w_1 + w_2} (|\mathbf{p}_1 - \mathbf{p}_2| - d) \frac{\mathbf{p}_1 - \mathbf{p}_2}{|\mathbf{p}_1 - \mathbf{p}_2|}\]

\[\Delta \mathbf{p}_2 = +\frac{w_2}{w_1 + w_2} (|\mathbf{p}_1 - \mathbf{p}_2| - d) \frac{\mathbf{p}_1 - \mathbf{p}_2}{|\mathbf{p}_1 - \mathbf{p}_2|}\]

If a particle is immovable, its inverse mass \(w\) is set to zero, receiving no positional displacement.

The Parallel Hazard: Gauss-Seidel vs. Jacobi

On a CPU, PBD typically uses a Gauss-Seidel solver, resolving constraints sequentially so subsequent constraints immediately read updated positions. On a GPU, multiple threads executing constraints simultaneously on shared vertices create write-after-read and write-after-write hazards.

Two strategies resolve this on the GPU:

  1. Jacobi Approach (Accumulation): Each constraint computes displacements \(\Delta \mathbf{p}\) and writes them atomically into accumulation buffers. A secondary shader pass divides the accumulated displacement by the number of connected constraints and updates positions. While trivially parallel, this method converges slowly and can exhibit damping or unnatural oscillations.
  2. Gauss-Seidel via Graph Coloring (Constraint Partitioning): Constraints are pre-processed into independent groups (batches) where no two constraints in the same group share a vertex. Threads execute constraints within a single group concurrently with zero race conditions. A synchronization barrier is placed between groups, preserving fast Gauss-Seidel convergence rates.

GLSL Compute Shader Implementation

The Gauss-Seidel method with graph coloring is standard for high-performance physics. The data is held in Shader Storage Buffer Objects (SSBOs).

1. Data Structures and Buffers

#version 450

layout(local_size_x = 256) in;

struct Constraint {
    uint indexA;
    uint indexB;
    float restLength;
    float padding;
};

layout(std430, binding = 0) buffer PositionBuffer {
    vec4 positions[]; // w component stores inverse mass (invMass)
};

layout(std430, binding = 1) readonly buffer ConstraintBuffer {
    Constraint constraints[];
};

uniform uint u_BatchOffset;
uniform uint u_BatchCount;

2. Constraint Resolution Kernel

Each invocation resolves one constraint from the active batch:

void main() {
    uint id = gl_GlobalInvocationID.x;
    if (id >= u_BatchCount) return;

    Constraint c = constraints[u_BatchOffset + id];
    
    vec4 pA = positions[c.indexA];
    vec4 pB = positions[c.indexB];

    vec3 delta = pA.xyz - pB.xyz;
    float currentDist = length(delta);
    
    if (currentDist < 1e-6) return; // Prevent division by zero

    float wA = pA.w;
    float wB = pB.w;
    float wSum = wA + wB;

    if (wSum == 0.0) return;

    float error = currentDist - c.restLength;
    vec3 correction = (error / wSum) * (delta / currentDist);

    positions[c.indexA].xyz -= correction * wA;
    positions[c.indexB].xyz += correction * wB;
}

Execution Loop and Synchronization

The host application drives the simulation step by coordinating compute dispatches and memory barriers:

  1. Predict Positions: A compute shader integrates external forces (e.g., gravity) and predicts candidate positions: \(\mathbf{p}^* = \mathbf{x} + \mathbf{v}\Delta t + \mathbf{f}_{ext}\Delta t^2\).
  2. Iterative Constraint Loop: For a fixed number of solver iterations (typically 10 to 30):
    • Loop through each constraint partition (color batch):
      • Set u_BatchOffset and u_BatchCount.
      • Dispatch the constraint-solving compute shader.
      • Issue an execution and memory barrier: glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT) to ensure that all position writes are visible before the next dependent batch reads them.
  3. Velocity Update and Position Commit: A final compute pass computes the new velocities \(\mathbf{v} = (\mathbf{p}^* - \mathbf{x}) / \Delta t\) and updates the previous positions \(\mathbf{x} = \mathbf{p}^*\).