{
"cells": [
{
"cell_type": "markdown",
"id": "feaeef18",
"metadata": {},
"source": [
"# 22 · Fast Mass Springs — a constant-Hessian implicit solver\n",
"\n",
"A mass-spring system stepped with **implicit (backward) Euler** solves, every step,\n",
"\n",
"$$\n",
"x_{n+1} = \\arg\\min_x \\; \\tfrac{1}{2h^2}\\,\\lVert x - y\\rVert_M^2 \\;+\\; E_\\text{spring}(x),\n",
"\\qquad y = x_n + h v_n + h^2 g .\n",
"$$\n",
"\n",
"The spring energy $E_\\text{spring}=\\sum_s \\tfrac{k_s}{2}\\,(\\lVert x_i-x_j\\rVert - r_s)^2$\n",
"is **non-convex** (the rest-length term), so a straight Newton solve must rebuild and\n",
"**re-factorize** the Hessian every iteration. Liu et al. show a better way: introduce\n",
"one auxiliary direction $d_s$ per spring and write\n",
"\n",
"$$\n",
"E_\\text{spring}(x) = \\min_{\\lVert d_s\\rVert=r_s}\\; \\sum_s \\tfrac{k_s}{2}\\,\\lVert (x_i-x_j) - d_s\\rVert^2 .\n",
"$$\n",
"\n",
"Now the step alternates a **local** projection (set each $d_s$ to the rest length\n",
"along the current spring direction) and a **global** solve. With the $d_s$ fixed the\n",
"objective is quadratic and its system matrix\n",
"\n",
"$$\n",
"\\big(\\tfrac{1}{h^2}M + G^\\top W G\\big)\n",
"$$\n",
"\n",
"is built only from rest data — **constant for the whole simulation**. Factorize it\n",
"**once**; every global step is then a back-substitution. That is the \"fast\" in fast\n",
"mass springs.\n",
"\n",
"**Reference.** T. Liu, A. W. Bargteil, J. F. O'Brien, L. Kavan, *\"Fast Simulation of\n",
"Mass-Spring Systems\"*, ACM TOG (SIGGRAPH Asia) 2013. The same local/global +\n",
"prefactored-global structure underlies projective dynamics (Bouaziz et al. 2014)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "f6f9367a",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:38:01.526430Z",
"iopub.status.busy": "2026-06-09T05:38:01.526266Z",
"iopub.status.idle": "2026-06-09T05:38:02.040992Z",
"shell.execute_reply": "2026-06-09T05:38:02.040720Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import time\n",
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"import simkit\n",
"import simkit.energies as energies\n",
"from simkit.edges import edges as mesh_edges\n",
"from simkit.solvers import block_coord\n",
"import utils"
]
},
{
"cell_type": "markdown",
"id": "9f16cf38",
"metadata": {},
"source": [
"## A cloth and its constant global matrix\n",
"\n",
"The cloth is a triangulated grid; its springs are the triangle edges (structural +\n",
"shear). We pin the two **top corners** and let it drape under gravity. The simulator\n",
"class below builds, once in `__init__`, the edge-difference operator $G$ (so $Gx$\n",
"stacks the spring vectors $x_j-x_i$), the spring weights $W=\\mathrm{diag}(k_s)\\otimes I$,\n",
"the lumped mass $M$, and the **constant** system matrix $H_c = M/h^2 + G^\\top W G$,\n",
"which it factorizes once with `sp.sparse.linalg.factorized`.\n",
"\n",
"The implicit step is block coordinate descent on the split energy: a `local_step`\n",
"that projects each spring onto its rest-length sphere\n",
"($d = r\\,(x_j-x_i)/\\lVert x_j-x_i\\rVert$) and a `global_step` that is the prefactored\n",
"solve $H_c\\,x = M y / h^2 + G^\\top W d$. We drive the two with simkit's `block_coord`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "bcc14573",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:38:02.042294Z",
"iopub.status.busy": "2026-06-09T05:38:02.042206Z",
"iopub.status.idle": "2026-06-09T05:38:02.046613Z",
"shell.execute_reply": "2026-06-09T05:38:02.046422Z"
}
},
"outputs": [],
"source": [
"class FastMassSpring:\n",
" \"\"\"Hanging cloth stepped with implicit Euler via local/global (Liu et al. 2013).\n",
"\n",
" The global system matrix Hc = M/h^2 + G^T W G is built from rest data only,\n",
" so it is constant for the whole run and factorized exactly once. Each implicit\n",
" step alternates a local projection (spring -> rest length) and a prefactored\n",
" global solve, both driven by simkit.solvers.block_coord.\n",
" \"\"\"\n",
"\n",
" def __init__(self, X, T, ym=80.0, h=1.0 / 30.0, gravity=9.8, corners_idx=None):\n",
" self.X, self.T = X, T\n",
" self.n, self.dim = X.shape\n",
" self.h = h\n",
"\n",
" # springs = unique triangle edges\n",
" self.Edge = mesh_edges(T)\n",
" self.ne = self.Edge.shape[0]\n",
"\n",
" # rest lengths, per-spring stiffness k = ym / l0^2\n",
" self.l0 = np.linalg.norm(X[self.Edge[:, 1]] - X[self.Edge[:, 0]], axis=1)[:, None]\n",
" self.ym = np.full((self.ne, 1), ym)\n",
" self.k = (self.ym / self.l0 ** 2).flatten()\n",
"\n",
" # edge-difference operator G : flattened x -> stacked spring vectors d_e = x_j - x_i\n",
" rows = np.repeat(np.arange(self.ne), 2)\n",
" cols = self.Edge.reshape(-1)\n",
" vals = np.tile([-1.0, 1.0], self.ne)\n",
" Sinc = sp.sparse.csc_matrix((vals, (rows, cols)), shape=(self.ne, self.n))\n",
" self.G = sp.sparse.kron(Sinc, sp.sparse.eye(self.dim)).tocsc()\n",
" self.W = sp.sparse.kron(sp.sparse.diags(self.k), sp.sparse.eye(self.dim)).tocsc()\n",
"\n",
" # lumped mass and gravity acceleration\n",
" self.M = sp.sparse.eye(self.n * self.dim).tocsc()\n",
" self.gvec = np.tile([0.0, -gravity], self.n).reshape(-1, 1)\n",
"\n",
" # pin the two top corners (default) and keep their DOFs fixed at rest\n",
" if corners_idx is None:\n",
" ytop = X[:, 1].max()\n",
" corners_idx = np.where((X[:, 1] >= ytop - 1e-9) &\n",
" ((X[:, 0] <= X[:, 0].min() + 1e-9) | (X[:, 0] >= X[:, 0].max() - 1e-9)))[0]\n",
" self.corners = corners_idx\n",
" self.pin_dof = (corners_idx[:, None] * self.dim + np.arange(self.dim)).ravel()\n",
" self.free_dof = np.setdiff1d(np.arange(self.n * self.dim), self.pin_dof)\n",
" self.xp = X.flatten().reshape(-1, 1)[self.pin_dof] # corners stay fixed at rest\n",
"\n",
" # THE CONSTANT global matrix Hc = M/h^2 + G^T W G (independent of x) -> factor once\n",
" self.Hc = (self.M / h ** 2 + self.G.T @ self.W @ self.G).tocsc()\n",
" self.H_ff = self.Hc[self.free_dof][:, self.free_dof].tocsc()\n",
" self.H_fp = self.Hc[self.free_dof][:, self.pin_dof].tocsc()\n",
" self.solve = sp.sparse.linalg.factorized(self.H_ff) # the ONLY factorization\n",
"\n",
" # ---- energy of the implicit step: inertia term + spring term ----\n",
" def energy(self, x, y):\n",
" d = (self.G @ x).reshape(-1, self.dim)\n",
" l = np.linalg.norm(d, axis=1, keepdims=True)\n",
" E_inertia = 0.5 / self.h ** 2 * float((x - y).T @ (self.M @ (x - y)))\n",
" E_spring = float(np.sum(0.5 * self.k[:, None] * (l - self.l0) ** 2))\n",
" return E_inertia + E_spring\n",
"\n",
" def _make_full(self, xf):\n",
" x = np.empty((self.n * self.dim, 1))\n",
" x[self.free_dof] = xf\n",
" x[self.pin_dof] = self.xp\n",
" return x\n",
"\n",
" # ---- local step: project each spring onto its rest-length sphere ----\n",
" def local_step(self, xf):\n",
" d = (self.G @ self._make_full(xf)).reshape(-1, self.dim)\n",
" d = self.l0 * d / np.linalg.norm(d, axis=1, keepdims=True)\n",
" return d.reshape(-1, 1)\n",
"\n",
" # ---- global step: the prefactored constant solve (back-substitution) ----\n",
" def global_step(self, xf, d, rhs_inertia):\n",
" rhs = (rhs_inertia + self.G.T @ (self.W @ d))[self.free_dof] - self.H_fp @ self.xp\n",
" return self.solve(rhs).reshape(-1, 1)\n",
"\n",
" def step(self, x, v, n_inner=8):\n",
" # one backward-Euler step via local/global (fast mass springs)\n",
" y = x + self.h * v + self.h ** 2 * self.gvec # momentum target\n",
" rhs_inertia = (self.M / self.h ** 2) @ y\n",
" global_step = lambda xf, d: self.global_step(xf, d, rhs_inertia)\n",
" xf = block_coord(x[self.free_dof].copy(), global_step, self.local_step,\n",
" max_iter=n_inner, tolerance=0.0)\n",
" x_new = self._make_full(xf)\n",
" v_new = (x_new - x) / self.h\n",
" return x_new, v_new\n",
"\n",
" def simulate(self, n_steps=240, n_inner=8):\n",
" x = self.X.flatten().reshape(-1, 1).copy()\n",
" v = np.zeros_like(x)\n",
" states = [x.reshape(self.n, self.dim).copy()]\n",
" for _ in range(n_steps):\n",
" x, v = self.step(x, v, n_inner=n_inner)\n",
" states.append(x.reshape(self.n, self.dim).copy())\n",
" return states\n",
"\n",
" # ---- a Newton implicit Hessian, for the comparison below ----\n",
" def newton_hessian(self, x):\n",
" d = (self.G @ x).reshape(-1, self.dim)\n",
" He = energies.mass_springs_hessian_element_d(d, self.ym, self.l0) # per-spring, depends on x\n",
" Hd = sp.sparse.block_diag(He).tocsc()\n",
" return (self.M / self.h ** 2 + self.G.T @ Hd @ self.G).tocsc()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "1b5834d8",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:38:02.047486Z",
"iopub.status.busy": "2026-06-09T05:38:02.047433Z",
"iopub.status.idle": "2026-06-09T05:38:02.052705Z",
"shell.execute_reply": "2026-06-09T05:38:02.052523Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"625 nodes, 1776 springs\n"
]
}
],
"source": [
"X, T = utils.triangulated_grid(nx=25, ny=25, width=2.0, height=2.0)\n",
"sim = FastMassSpring(X, T, ym=80.0, h=1.0 / 30.0, gravity=9.8)\n",
"print(f\"{sim.n} nodes, {sim.ne} springs\")"
]
},
{
"cell_type": "markdown",
"id": "00183198",
"metadata": {},
"source": [
"## \"Constant Hessian\": local/global vs. Newton\n",
"\n",
"The local/global matrix $H_c$ is built only from rest data, so it never changes\n",
"— one factorization serves the whole run. A Newton implicit solve instead uses\n",
"$M/h^2 + G^\\top H_\\text{spring}(x)\\,G$, whose spring block depends on the current\n",
"positions and must be reassembled and re-factorized at every iteration. Let's verify."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8b1de8ba",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:38:02.053658Z",
"iopub.status.busy": "2026-06-09T05:38:02.053599Z",
"iopub.status.idle": "2026-06-09T05:38:02.063199Z",
"shell.execute_reply": "2026-06-09T05:38:02.062998Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Newton implicit Hessian, ||H(x_a) - H(x_b)|| = 9.610e+05 (nonzero -> rebuild + refactorize every iter)\n",
"Local/global matrix Hc is state-independent = 0.000e+00 (always 0 -> factor once)\n"
]
}
],
"source": [
"rng = np.random.default_rng(0)\n",
"x_a = X.flatten().reshape(-1, 1)\n",
"x_b = (X + 0.3 * rng.standard_normal(X.shape)).flatten().reshape(-1, 1)\n",
"dHn = sp.sparse.linalg.norm(sim.newton_hessian(x_a) - sim.newton_hessian(x_b))\n",
"print(f\"Newton implicit Hessian, ||H(x_a) - H(x_b)|| = {dHn:.3e} (nonzero -> rebuild + refactorize every iter)\")\n",
"print(f\"Local/global matrix Hc is state-independent = {sp.sparse.linalg.norm(sim.Hc - sim.Hc):.3e} (always 0 -> factor once)\")"
]
},
{
"cell_type": "markdown",
"id": "684fe69a",
"metadata": {},
"source": [
"## Drop the cloth\n",
"\n",
"Released flat from rest, held at the two top corners, settling under gravity. The\n",
"solver is unconditionally stable (implicit), so a large $h=1/30$ and only 8 inner\n",
"local/global iterations per step are plenty."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "53db2e81",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:38:02.064247Z",
"iopub.status.busy": "2026-06-09T05:38:02.064200Z",
"iopub.status.idle": "2026-06-09T05:38:02.227676Z",
"shell.execute_reply": "2026-06-09T05:38:02.227466Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"240 implicit steps in 0.16 s (one factorization reused throughout)\n"
]
}
],
"source": [
"n_steps = 240\n",
"t0 = time.perf_counter()\n",
"states = sim.simulate(n_steps=n_steps, n_inner=8)\n",
"print(f\"{n_steps} implicit steps in {time.perf_counter()-t0:.2f} s (one factorization reused throughout)\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "c05e2a44",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:38:02.228559Z",
"iopub.status.busy": "2026-06-09T05:38:02.228507Z",
"iopub.status.idle": "2026-06-09T05:38:04.519542Z",
"shell.execute_reply": "2026-06-09T05:38:04.519280Z"
}
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def thin(states, target=120):\n",
" k = max(1, len(states) // target); idx = list(range(0, len(states), k))\n",
" if idx[-1] != len(states) - 1: idx.append(len(states) - 1)\n",
" return [states[i] for i in idx]\n",
"\n",
"lims = ((-1.6, 1.6), (-2.7, 1.3))\n",
"fig, anim = utils.animate_mesh(thin(states), T, lims=lims, fps=30,\n",
" title=\"Hanging cloth, fast mass springs (implicit, prefactored)\", pin_pts=X[sim.corners])\n",
"utils.show_video(fig, anim, \"media/22_cloth.mp4\", fps=30)"
]
},
{
"cell_type": "markdown",
"id": "d7846964",
"metadata": {},
"source": [
"## What the constant matrix buys per step\n",
"\n",
"Per implicit step, the local/global solver does a handful of back-substitutions\n",
"through one cached factorization; a Newton step instead reassembles and\n",
"re-factorizes its state-dependent Hessian at every inner iteration. Same problem,\n",
"very different per-step cost."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ee4417cb",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:38:04.521782Z",
"iopub.status.busy": "2026-06-09T05:38:04.521699Z",
"iopub.status.idle": "2026-06-09T05:38:04.568531Z",
"shell.execute_reply": "2026-06-09T05:38:04.568329Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"per step (8 inner iters): fast local/global 0.9 ms vs refactorizing 43.5 ms\n"
]
}
],
"source": [
"xt = states[5].flatten().reshape(-1, 1)\n",
"vt = (states[5] - states[4]).flatten().reshape(-1, 1) / sim.h\n",
"y = xt + sim.h * vt + sim.h ** 2 * sim.gvec\n",
"rhs_inertia = (sim.M / sim.h ** 2) @ y\n",
"\n",
"# fast: 8 prefactored back-substitutions (reuse the one factorization)\n",
"t0 = time.perf_counter()\n",
"for _ in range(8):\n",
" d = (sim.G @ xt).reshape(-1, sim.dim); d = sim.l0 * d / np.linalg.norm(d, axis=1, keepdims=True)\n",
" _ = sim.solve((rhs_inertia + sim.G.T @ (sim.W @ d.reshape(-1, 1)))[sim.free_dof] - sim.H_fp @ sim.xp)\n",
"t_fast = time.perf_counter() - t0\n",
"\n",
"# newton-style: assemble + factorize the state-dependent Hessian each iteration\n",
"t0 = time.perf_counter()\n",
"for _ in range(8):\n",
" Hn = sim.newton_hessian(xt)[sim.free_dof][:, sim.free_dof].tocsc()\n",
" _ = sp.sparse.linalg.spsolve(Hn, np.ones((len(sim.free_dof), 1)))\n",
"t_newton = time.perf_counter() - t0\n",
"\n",
"print(f\"per step (8 inner iters): fast local/global {1e3*t_fast:6.1f} ms vs refactorizing {1e3*t_newton:6.1f} ms\")"
]
},
{
"cell_type": "markdown",
"id": "00d8e86f",
"metadata": {},
"source": [
"## Takeaways\n",
"\n",
"* **The reformulation makes the system matrix constant.** Splitting each spring's\n",
" energy with an auxiliary direction turns the global step into a quadratic solve\n",
" with a fixed matrix $M/h^2 + G^\\top W G$. Factor it once at setup; every step and\n",
" every inner iteration thereafter is a cheap back-substitution.\n",
"* **This is the same machinery as tutorial 21.** Local step = closed-form per-element\n",
" projection (here, a spring onto its rest length; for ARAP, a triangle onto its\n",
" closest rotation). Global step = one prefactored linear solve. Block coordinate\n",
" descent (`simkit.solvers.block_coord`) drives both.\n",
"* **Implicit + prefactored = fast and stable.** A big timestep stays stable, and the\n",
" per-step cost is dominated by sparse back-substitutions rather than repeated\n",
" factorizations — exactly the win reported by Liu et al. (2013)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "simkit",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}