Position-Based Dynamics (PBD)

In a cloth simulation, each particle is influenced by gravity when external forces like collisions or spring constraints are ignored. The velocity and position are updated accordingly:

float dt = 0.0333f;
float damping = 0.99f;
Vector3[] positions; // X
Vector3[] velocities; // V
Vector3 g = new Vector3(0.0f, -9.8f, 0.0f);
int[] edges; // E
float[] restLengths; // L

for (int i = 0; i < positions.Length; i++)
{
    if (i == 0 || i == 20) continue; // Fixed points
    velocities[i] += g * dt;
    velocities[i] *= damping;
    positions[i] += velocities[i] * dt;
}

Traditional spring models based on Hooke’s law require increasing stiffness to simulate realistic resistance under large deformations. However, this leads to numerical instability in explicit solvers and increased computational cost in implicit ones. Position-Based Dynamics (PBD) addresses this by enforcing constraints directly on positions rather then modeling forces.

A Jacobi-style constraint solver iteratively adjusts vertex positions to satisfy distance constraints:

void ApplyStrainLimiting()
{
    Mesh mesh = GetComponent<MeshFilter>().mesh;
    Vector3[] verts = mesh.vertices;
    Vector3[] corrected = new Vector3[verts.Length];
    int[] count = new int[verts.Length];

    for (int i = 0; i < verts.Length; i++)
    {
        corrected[i] = Vector3.zero;
        count[i] = 0;
    }

    for (int e = 0; e < restLengths.Length; e++)
    {
        int i0 = edges[e * 2 + 0];
        int i1 = edges[e * 2 + 1];
        Vector3 dir = verts[i0] - verts[i1];
        float currentLen = dir.magnitude;
        float error = (currentLen - restLengths[e]) * 0.5f;
        Vector3 correction = error * dir.normalized;

        corrected[i0] += verts[i0] - correction;
        corrected[i1] += verts[i1] + correction;
        count[i0]++;
        count[i1]++;
    }

    for (int i = 0; i < verts.Length; i++)
    {
        if (i == 0 || i == 20) continue;
        float weight = count[i] + 0.2f;
        Vector3 newPos = (corrected[i] + 0.2f * verts[i]) / weight;
        velocities[i] += (newPos - verts[i]) / dt;
        verts[i] = newPos;
    }

    mesh.vertices = verts;
}

Collision Handling with a Sphere

After PBD constraint resolution, per-vertex collision detection against a sphere is performed within the same frame:

void HandleSphereCollision()
{
    Mesh mesh = GetComponent<MeshFilter>().mesh;
    Vector3[] X = mesh.vertices;
    Vector3 center = sphereTransform.position;
    float radius = sphereRadius;

    for (int i = 0; i < X.Length; i++)
    {
        if (i == 0 || i == 20) continue;
        Vector3 toPoint = X[i] - center;
        float dist = toPoint.magnitude;
        if (dist > radius) continue;

        Vector3 surfacePoint = center + toPoint.normalized * radius;
        Vector3 normal = (surfacePoint - center).normalized;
        Vector3 vNormal = Vector3.Dot(velocities[i], normal) * normal;
        Vector3 vTangent = velocities[i] - vNormal;

        velocities[i] = vTangent + vNormal + (surfacePoint - X[i]) / dt;
        X[i] = surfacePoint;
    }

    mesh.vertices = X;
}

Implicit Integration Approach

An alternative formulation treats the time step as an optimization problem: minimize an energy function $ F(\mathbf{x}) $. When the Hessian of $ F $ is positive definite, a unique minimum exists.

The gradient combines inertial, gravitational, and elastic terms:

void ComputeGradient(Vector3[] x, Vector3[] xPrev, float dt, Vector3[] grad)
{
    float invDt2 = 1.0f / (dt * dt);
    for (int i = 0; i < x.Length; i++)
    {
        grad[i] = mass * (invDt2 * (x[i] - xPrev[i]) - gravity);
    }

    for (int e = 0; e < restLengths.Length; e++)
    {
        int i = edges[e * 2 + 0];
        int j = edges[e * 2 + 1];
        Vector3 diff = x[i] - x[j];
        Vector3 dir = diff.normalized;
        Vector3 force = springK * (diff - restLengths[e] * dir);
        grad[i] += force;
        grad[j] -= force;
    }
}

Using a diagonal approximation of the Hessian (a "magic" constant), the update becomes:

for (int iter = 0; iter < 32; iter++)
{
    ComputeGradient(x, xPrev, dt, grad);
    for (int i = 0; i < x.Length; i++)
    {
        if (i == 0 || i == 20) continue;
        float hessInv = 1.0f / (mass / (dt * dt) + 4.0f * springK);
        x[i] -= hessInv * grad[i];
    }
}

// Update velocities
for (int i = 0; i < velocities.Length; i++)
{
    velocities[i] = (x[i] - xPrev[i]) / dt;
}

This naive fixed-iteration approach is slow to converge. Acceleration can be achieved using Chebyshev-accelerated Jacobi iteration:

Vector3[] lastX = new Vector3[x.Length];
Array.Copy(x, lastX, x.Length);
float w = 1.0f, rho = 0.9f; // spectral radius estimate

for (int iter = 0; iter < 32; iter++)
{
    ComputeGradient(x, xPrev, dt, grad);
    bool converged = true;

    if (iter == 0) w = 1.0f;
    else if (iter == 1) w = 2.0f / (2.0f - rho * rho);
    else w = 4.0f / (4.0f - rho * rho * w);

    for (int i = 0; i < x.Length; i++)
    {
        if (i == 0 || i == 20) continue;
        float hessInv = 1.0f / (mass / (dt * dt) + 4.0f * springK);
        Vector3 old = x[i];
        x[i] -= hessInv * grad[i];
        x[i] = w * x[i] + (1 - w) * lastX[i];
        lastX[i] = old;

        if ((x[i] - old).magnitude > 1e-4f)
            converged = false;
    }

    if (converged) break;
}

Collision handling for the implicit method follows the same geometric projection strategy as in PBD.