Verlet Physics 3D — Mass-Spring Cloth

Math essence

Verlet integration: x_new = 2·x - x_prev + a·dt². Damping: x_new = x + (1-d)·(x - x_prev) + a·dt². Constraint relaxation: for each spring with rest length L and stiffness k, move each end by ±(k/2)·(current_length - L)·normalized_direction. Bending constraints: a chain (a, b, c) gets rest_angle preservation via a perpendicular-offset projection. NaN/Infinity guard on every step. Floor plane: y >= 0 with reflection (v_new = -v_old * restitution).

Parameters

Reference

First 60 lines of reference.py. Full source: library_external/toxi_toxiclibs/VerletPhysics3D/reference.py

"""VerletPhysics3D — 3D mass-spring cloth simulation (Jakobsen 2001 style).

A pure-numpy port of Karsten Schmidt's toxiclibs.physics.VerletPhysics3D
(Karsten Schmidt, 2008). The algorithm uses position-based dynamics:

  * Each particle stores (x, x_prev) instead of (x, v).
  * Integration: x_new = x + (x - x_prev) * damping + a * dt²
  * Constraint relaxation: for each spring, project both ends toward
    the rest length by ±(k/2) * (current - rest) * direction.

Bending constraints span 3-particle chains and preserve the rest
angle; the floor is a hard plane (y >= 0) with restitution.

The 2D version is `toxi_openframeworks/ParticleSpring2D`. This is the
3D extension.
"""
from __future__ import annotations

import numpy as np


class VerletPhysics3D:
    """3D Verlet mass-spring system with cloth-grid initializer."""

    def __init__(self, n_grid=14, size=4.0, gravity=0.3, damping=0.98,
                 iterations=12, stiffness=1.0, wind=0.0, seed=0):
        self.n = n_grid
        self.size = size
        self.gravity = gravity
        self.damping = damping
        self.iterations = iterations
        self.stiffness = stiffness
        self.wind = wind
        self.rng = np.random.default_rng(seed)

        # Build a square cloth grid in the XZ plane at y = size/2.
        spacing = size / n_grid
        half = size / 2
        coords = np.array(
            [[(i - n_grid/2) * spacing, half, (j - n_grid/2) * spacing]
             for j in range(n_grid) for i in range(n_grid)],
            dtype=np.float32,
        )
        n_particles = n_grid * n_grid
        self.n_particles = n_particles
        self.pos = coords.copy()
        self.pos_prev = coords.copy()
        # Pin the four corners (so the cloth hangs)
        self.locked = np.zeros(n_particles, dtype=bool)
        for (i, j) in [(0, 0), (n_grid-1, 0), (0, n_grid-1), (n_grid-1, n_grid-1)]:
            self.locked[j * n_grid + i] = True
        # Slight initial perturbation so the relaxation has something to do
        self.pos += self.rng.normal(0, 0.01, size=self.pos.shape).astype(np.float32)

        # Build springs: structural (4-neighbors) + shear (diagonals)
        # + bending (skip-one neighbors). Each is a (a, b) pair.
        springs = []
        for j in range(n_grid):
            for i in range(n_grid):
                idx = j * n_grid + i