Implicit Time Integration

To simulate cloth realistically, explicit integration often requires tiny time steps to remain stable. Instead, we employ implicit integration, where the state at the next time step depends on forces calculated at that future state. The semi-implicit Euler equation is defined as:

$$v_{new} = v_{old} + \Delta t M^{-1} f(x_{new})$$
$$x_{new} = x_{old} + \Delta t v_{new}$$

Combining these yields a non-linear system of equations where $x_{new}$ appears on both sides inside the force function $f$:

$$x_{new} = x_{old} + \Delta t v_{old} + \Delta t^2 M^{-1} f(x_{new})$$

Energy Minimization Formulation

Rather than finding the root of the momentum equation directly, we can reframe the problem as minimizing the total energy of the system. We define an objective function $F(x)$:

$$F(x) = \frac{1}{2\Delta t^2} \|x - x_{old} - \Delta t v_{old}\|_M^2 + E_{potential}(x)$$

Where $\|x\|_M^2 = x^T M x$. The minimum of this energy function corresponds to the solution of the implicit integration step. Therefore, finding the optimal position $x_{new}$ is equivalent to solving $\nabla F(x_{new}) = 0$.

Newton-Raphson Optimization

To solve the minimization problem, we use the Newton-Raphson method. We approximate the function locally using a second-order Taylor expansion. Given a current guess $x^{(k)}$, we compute the step $\Delta x$ by solving:

$$H(x^{(k)}) \Delta x = -\nabla F(x^{(k)})$$
$$x^{(k+1)} = x^{(k)} + \Delta x$$

Here, $H(x)$ is the Hessian matrix (second derivatives), and $\nabla F(x)$ is the gradient. For the mass-spring system, these components are derived as follows:

$$\nabla F(x) = \frac{1}{\Delta t^2} M (x - x_{pred}) - f(x)$$
$$H(x) = \frac{1}{\Delta t^2} M + \frac{\partial f}{\partial x}$$

Where $x_{pred} = x_{old} + \Delta t v_{old}$. Constructing the exact Hessian for a cloth mesh is computationally expensive. A common approximation in computer graphics is to use a diagonalized Hessian, treating springs as independent entities, which simplifies the update rule significantly.

Mesh and Edge Initialization

The following C# code sets up a grid mesh and extracts structural edges for the simulation. It avoids duplicate edges by sorting and filtering.

void InitializeSimulation()
{
    MeshFilter filter = GetComponent<MeshFilter>();
    Mesh mesh = filter.mesh;

    int gridSize = 21;
    int vertCount = gridSize * gridSize;
    Vector3[] vertices = new Vector3[vertCount];
    Vector2[] uvs = new Vector2[vertCount];

    // Create Grid
    for (int y = 0; y < gridSize; y++)
    {
        for (int x = 0; x < gridSize; x++)
        {
            int idx = y * gridSize + x;
            float fx = (float)x / (gridSize - 1);
            float fy = (float)y / (gridSize - 1);
            vertices[idx] = new Vector3((fx - 0.5f) * 10, 0, (fy - 0.5f) * 10);
            uvs[idx] = new Vector2(fx, fy);
        }
    }

    // Create Triangles
    int trisCount = (gridSize - 1) * (gridSize - 1) * 6;
    int[] triangles = new int[trisCount];
    int t = 0;
    for (int y = 0; y < gridSize - 1; y++)
    {
        for (int x = 0; x < gridSize - 1; x++)
        {
            int tl = y * gridSize + x;
            int tr = tl + 1;
            int bl = (y + 1) * gridSize + x;
            int br = bl + 1;

            triangles[t++] = tl; triangles[t++] = bl; triangles[t++] = tr;
            triangles[t++] = tr; triangles[t++] = bl; triangles[t++] = br;
        }
    }

    mesh.vertices = vertices;
    mesh.triangles = triangles;
    mesh.uv = uvs;
    mesh.RecalculateNormals();

    // Extract Unique Edges
    List<(int, int)> rawEdges = new List<(int, int)>();
    for (int i = 0; i < triangles.Length; i += 3)
    {
        int a = triangles[i];
        int b = triangles[i + 1];
        int c = triangles[i + 2];
        rawEdges.Add(TupleSort(a, b));
        rawEdges.Add(TupleSort(b, c));
        rawEdges.Add(TupleSort(c, a));
    }

    // Sort edges to remove duplicates
    rawEdges.Sort((e1, e2) => 
    {
        if (e1.Item1 != e2.Item1) return e1.Item1.CompareTo(e2.Item1);
        return e1.Item2.CompareTo(e2.Item2);
    });

    List<(int, int)> uniqueEdges = new List<(int, int)>();
    for (int i = 0; i < rawEdges.Count; i++)
    {
        if (i == 0 || rawEdges[i] != rawEdges[i - 1])
        {
            uniqueEdges.Add(rawEdges[i]);
        }
    }

    // Store Edge Data
    springIndices = new int[uniqueEdges.Count * 2];
    restLengths = new float[uniqueEdges.Count];

    for (int i = 0; i < uniqueEdges.Count; i++)
    {
        int idxA = uniqueEdges[i].Item1;
        int idxB = uniqueEdges[i].Item2;
        springIndices[i * 2] = idxA;
        springIndices[i * 2 + 1] = idxB;
        restLengths[i] = Vector3.Distance(vertices[idxA], vertices[idxB]);
    }

    // Initialize velocities
    velocities = new Vector3[vertCount];
}

(int, int) TupleSort(int a, int b)
{
    return a < b ? (a, b) : (b, a);
}

Simulation Loop

The update loop performs the implicit integration. It uses a simplified solver where the Hessian is approximated as a diagonal matrix, avoiding the need for a global linear system solver.

void UpdateStep()
{
    Mesh mesh = GetComponent<MeshFilter>().mesh;
    Vector3[] currentPos = mesh.vertices;
    Vector3[] previousPos = (Vector3[])currentPos.Clone();
    Vector3[] predictedPos = new Vector3[currentPos.Length];
    Vector3[] gradient = new Vector3[currentPos.Length];

    float dt = Time.deltaTime;
    float invDt = 1.0f / dt;
    float invDtSq = invDt * invDt;

    // Predict position and apply damping
    for (int i = 0; i < currentPos.Length; i++)
    {
        // Pin two corners (e.g., top left and top right)
        if (IsPinned(i)) continue;

        velocities[i] *= dampingFactor;
        predictedPos[i] = currentPos[i] + velocities[i] * dt;
        currentPos[i] = predictedPos[i]; // Initial guess for Newton solver
    }

    // Precompute the inverse of the diagonal Hessian approximation
    // H_approx = (1/dt^2 * M) + 4*k
    float diagonalStiffness = invDtSq * mass + 4.0f * stiffness;
    float invHessian = 1.0f / diagonalStiffness;

    // Newton Iteration
    for (int iter = 0; iter < maxIterations; iter++)
    {
        CalculateEnergyGradient(currentPos, predictedPos, dt, gradient);

        for (int i = 0; i < currentPos.Length; i++)
        {
            if (IsPinned(i)) continue;
            // Newton update: x = x - H^-1 * Grad
            currentPos[i] -= invHessian * gradient[i];
        }
    }

    // Update velocities
    for (int i = 0; i < currentPos.Length; i++)
    {
        if (IsPinned(i)) continue;
        velocities[i] = invDt * (currentPos[i] - previousPos[i]);
    }

    mesh.vertices = currentPos;
    ResolveCollisions(currentPos, velocities, dt);
    mesh.RecalculateNormals();
}

bool IsPinned(int index)
{
    // Example: Pin first and last vertex of the first row
    return index == 0 || index == 20;
}

Gradient Calculation

The gradient of the objective function consists of the inertial term (difference between current and predicted position) and the internal spring forces.

void CalculateEnergyGradient(Vector3[] pos, Vector3[] pred, float dt, Vector3[] outGradient)
{
    float invDtSq = 1.0f / (dt * dt);

    // 1. Inertial Force and Gravity
    for (int i = 0; i < pos.Length; i++)
    {
        Vector3 delta = pos[i] - pred[i];
        outGradient[i] = invDtSq * mass * delta - mass * gravity;
    }

    // 2. Spring Forces
    for (int i = 0; i < springIndices.Length; i += 2)
    {
        int idxA = springIndices[i];
        int idxB = springIndices[i + 1];
        int edgeIdx = i / 2;

        Vector3 pA = pos[idxA];
        Vector3 pB = pos[idxB];
        Vector3 diff = pA - pB;
        float dist = diff.magnitude;

        if (dist < 1e-6f) continue; // Avoid division by zero

        // Hooke's Law: F = k * (current_len - rest_len)
        float stretch = dist - restLengths[edgeIdx];
        Vector3 force = (stiffness * stretch / dist) * diff;

        // Apply forces (Newton's 3rd Law)
        outGradient[idxA] += force;
        outGradient[idxB] -= force;
    }
}

Collision Handling

We use a Signed Distance Field (SDF) approach for simple spherical collision. When a vertex penetrates the collider, its position is projected to the surface, and its velocity is adjusted to reflect the impact.

void ResolveCollisions(Vector3[] pos, Vector3[] vel, float dt)
{
    Vector3 sphereCenter = Vector3.zero; // Assuming collider at origin
    float colliderRadius = 2.7f;

    for (int i = 0; i < pos.Length; i++)
    {
        if (IsPinned(i)) continue;

        Vector3 toVertex = pos[i] - sphereCenter;
        float dist = toVertex.magnitude;

        if (dist < colliderRadius)
        {
            Vector3 normal = toVertex / dist; // Normalized direction
            Vector3 surfacePoint = sphereCenter + normal * colliderRadius;

            // Positional Correction
            Vector3 correction = surfacePoint - pos[i];
            pos[i] = surfacePoint;

            // Velocity Update (Impulse)
            vel[i] += correction / dt;
        }
    }
}