Data-driven subspace simulation

140 nodes · 234 triangles · 280 positional DOFs

Building the full-space solver.

It starts with one triangle. Rest edges form D_m; current edges form D_s. Their ratio F = D_sD_m⁻¹ measures local stretch, shear, and rotation. J = det(F) is the signed area ratio and must remain positive.

Differentiating ψ with respect to F gives the first Piola–Kirchhoff stress P = μ(F − F^−T) + λ ln(J)F^−T. P is the derivative ∂ψ/∂F; it converts an element's deformation into forces measured against its rest shape. Multiplying P by the rest-shape gradients yields three vertex forces whose sum is zero before external forces are added.

let F = Ds * inverse(Dm);
let J = determinant(F);
let inverseTranspose = transpose(inverse(F));
let P = mu * (F - inverseTranspose)
      + lambda * log(J) * inverseTranspose;

// Each rest-shape gradient maps P to one nodal force.
triangleForce[i] = -restArea * P * gradN[i];

01B · Many triangles become one coupled system.

A triangle only knows about its three vertices. The body appears when all 234 triangles share vertices and their energy contributions are added.

Gravity, contact, and damping join that elastic force, giving Mẍ = f(x, ẋ). The mesh now has 280 coupled positional DOFs. This is the full-space problem we want to reduce.

Choose how to build the reduced space

Let’s collect the data first.

Reduced models can be data-driven or data-free. Sharp et al., for example, learn reduced-order kinematics directly from the governing energy, without a snapshot dataset.^[8] This article follows the data-driven route: record full-order meshes, then use them to construct the PCA basis and the conditioned neural decoder.

Four area-preserving launch deformations are run at μ = 7 and μ = 15. We save the initial mesh and eleven later snapshots from each trajectory. The resulting matrix X ∈ ℝ^96×280contains the full positions the reduced models should represent.

copyBufferToBuffer(position, dataset[firstSample]);
for (let step = 0; step < 440; step++) {
  advanceFullOrder(stiffness);
  if ((step + 1) % 40 === 0) {
    copyBufferToBuffer(position, dataset[nextSample]);
  }
}

96 recorded meshes · 280 position values each

How many ways did the body really move?

Overlay the 96 meshes we just collected. Every snapshot contains 280 position values, yet those values do not wander independently: when the body bends, many vertices move together. If those coordinated displacements repeat, we can store the patterns once and describe each mesh with only a few weights.

First average the recorded meshes to get x̄. Subtract that average from each snapshot and what remains is a displacement field: one arrow per vertex showing how that shape differs from x̄. PCA searches all 96 fields for the directions that account for the most variation. Each direction u_i is a whole-body displacement mode—not one vertex moving, but all 140 moving together in a learned pattern.

For this experiment, each position is weighted by the square root of its lumped nodal mass before PCA, so the fit respects how mass is distributed across the body. Keeping the first k modes forms U_k = [u₁ … u_k]. A mesh is then rebuilt as x(q) = x̄ + U_kq, where q contains the k signed mode weights.

Recorded motion is one way to find these directions. Classical modal analysis derives small-vibration modes around an equilibrium; modal derivatives extend them toward larger deformations. We use mass-weighted PCA here because it extracts the recurring motions from the data we just collected.^[12]

var x = mean[node];
for (var k = 0u; k < latentK; k++) {
  x += basis[node * MAX_K + k] * q[k];
}

Build a full deformation by hand.

The coral outline below is one recorded mesh. Each slider changes one q_i, scaling an entire displacement field. Four numbers produce all 280 reconstructed position values.

We can rebuild a deformation from q. How do those weights become a simulation?

The decoded shape x(q) = x̄ + Uq describes every mesh shape available to the reduced body. We can then evaluate elasticity, gravity, damping, and contact on that full mesh, producing the 280-value force vector f.

What should this force do to the few weights q? Take one weight q_k and increase it slightly. Every vertex moves along mode u_k. The dot product u_k^Tf measures how much the full force favors that motion, giving the force acting on q_k. A positive value accelerates the body toward larger q_k; a negative value accelerates it toward smaller q_k. For conservative forces, the same number is −∂E/∂q_k: the downhill slope of the full energy along that mode.

fq = transpose(U) * fullForce;
// Uᵀ M U = I after mass-weighted PCA.
qAcceleration = fq;
qVelocity += dt * qAcceleration;
q += dt * qVelocity;

The linear subspace has one fixed set of directions.

U is the same for every shape and for both material conditions. A deformation family that bends away from this plane requires additional modes. Material, pose, boundary conditions, and actuation must also fit inside the same fixed generalized-coordinate system. This is where we replace the plane with a learned deformation map.

From a fixed basis to a learned manifold

Going neural.

Keep q as the compact coordinates and replace x̄ + Uq with a neural decoder x(q, c). Its local motion directions can change across the deformation family and respond to a known condition such as material stiffness c. During simulation, these local directions take the role that the PCA modes played in force projection and reduced dynamics.^[17]

Before the implementation, let's recap the architecture.

This experiment uses an auto-decoder: a latent-variable model containing a decoder network and a table of trainable latent codes. Training instance i owns an unknown code q_i. Gradient descent optimizes every q_i together with the shared decoder weights θ so D_θ(q_i, c_i) reconstructs the observed target x_i. An autoencoder obtains its code from a learned encoder q_i = E_φ(x_i), whose weights are φ.^[5]

The shared weights describe the learned family; each code selects one point on it. Several observations can share condition c while their optimized codes give them distinct latent locations. A new measured observation can receive a q by optimizing it with θ fixed or through an encoder. During simulation, q becomes the dynamic configuration advanced by the physics solve.

04A · Use the latent code as this body's generalized coordinates.

Here x is a 280-value mesh configuration, q ∈ ℝ^k is its vector of learned generalized coordinates, and c is a known material condition. The MLP g_θ maps (q, c) to a 280-value displacement. The rest configuration x_rest anchors the chart:

04B · Let material stiffness reshape the learned subspace.

A condition is an extra known input that lets the decoder change its family of shapes with a property we choose. Here that property is material stiffness. During simulation, stiffness does two jobs. The Neo-Hookean model uses it to calculate forces, and the decoder uses it to shape the space in which those forces are solved. A soft and stiff snapshot recorded from the same launch and time share one latent code q_j. The code identifies the motion sample; the stiffness condition tells the decoder which material version of that state to produce.

The dataset therefore has 96 snapshots but only 48 codes. Values between the two training stiffnesses ask the decoder to interpolate. Sharing codes assumes the paired trajectories remain comparable over time; the diagnostic in Section 02 shows when they begin to separate.

PCA initializes the latent codes and a near-linear set of MLP weights so browser training begins from a well-scaled chart. Adam then updates the codes and every decoder weight. The trained runtime map is the MLP, and its learned launch code initializes q at the start of a rollout. PCA supplies the initialization. This low-to-full neural kinematic map is the subspace representation used by direct neural reduced models.^[6][8]

let base = tanh(bias[h] + conditionWeight[h] * c);
let hidden = tanh(bias[h] + conditionWeight[h] * c
                + dot(latentWeight[h], qCode[j]));

scaledPrediction[d] += outputWeight[d, h] * (hidden - base);
scaledTarget[d] = 4.0 * (pairedSnapshot[j][materialIndex][d] - rest[d]);
error[d] = scaledPrediction[d] - scaledTarget[d];

adamUpdate(qCode[j], latentGradient);
adamUpdate(weights, decoderGradient);

The factor four is a numerical training scale used by the shader; runtime decoding multiplies the MLP output by ¼ and restores the displacement to its physical scale. Across 3,000 Adam steps, the weight learning rate decays from 0.008 toward 0.0008 and the code learning rate from 0.00025 toward 0.000025.

04C · Turn the learned map into a reduced solver.

At the current q, each column of the decoder Jacobian J_D(q, c) = ∂x/∂q is one 280-value motion produced by changing a single neural coordinate. These state-dependent directions take the place of the fixed PCA modes. Their dot products with the assembled force give the generalized forces, and their mass-weighted dot products give the reduced inertia.^[6][7][17]

G = transpose(W2) * M * W2;                 // cache after training
JD = W2 * diag(tanhPrime(z)) * Wq / 4;        // exact 280 × k Jacobian
A = diag(tanhPrime(z)) * Wq;                  // 24 × k local tangent
b = tanhSecond(z) * square(Wq * qVelocity);   // hidden curvature

hiddenForce = transpose(W2) * fullForce - G * b / 4;
fq = transpose(A) * hiddenForce / 4;
Mq = transpose(A) * G * A / 16;

The live solver uses two common safeguards. Light regularization keeps the reduced solve stable when the decoder's local motion directions become nearly redundant. A trust region keeps q near the part of the learned space represented in training. The validation panel also repeats a new motion without the main regularizer, showing how strongly that safeguard affects the rollout.

Both rollouts are compared against the full simulation over time. Training is accepted only if the meshes remain finite, preserve triangle orientation, and stay physically bounded.

Reconstructing recorded meshes is only the first test. A useful neural subspace must also remain smooth enough for forces to project through it and stable enough to roll forward. That is why we judge the decoder by complete simulations as well as reconstruction loss.^[9]

One simulation step, side by side

Stage breakdown.

A reduced simulation is a hybrid pipeline. Geometry, material forces, contact, and rendering still use the complete 140-vertex mesh. The compact coordinates take over for the acceleration solve and integration, then produce the next complete mesh.

The main difference is the route through one step. The full-order method carries x and v straight through the solve. PCA and the auto-decoder first decode q, assemble the same mesh forces, project them back to q, and advance q before decoding again. The table follows those routes stage by stage.

Here x and v are full position and velocity vectors; q and q̇ are reduced position and velocity; a and q̈ are their accelerations; subscript 0 marks initialization; c is the fixed material condition; J_D = ∂x/∂q is the decoder Jacobian; and h is one integration substep.

In this experiment, the reduced models can use less time resolution. The full mesh contains fast, element-scale deformation modes, so one 1 / 90 s display frame is divided into eight 1 / 720 s updates. PCA and the neural subspace omit those high-frequency directions and remain stable with one 1 / 90 s update. That removes seven force, projection, and integration passes per frame in addition to shrinking the state from 280 positional DOFs to k generalized coordinates. On larger, finer meshes, the full-space stability limit and the cost of each solve make both savings increasingly important.

The live comparison begins in matched-contact mode, where all three methods use the same floor-and-wall penalty forces. A stabilized-reference option adds node-wise position and velocity projection to the full solver after integration. RMS values in that mode include the resulting contact-model difference. Robust reduced contact is a separate model-reduction problem.^[14]