Simulation

Simulation#

This page describes the forward physics simulation in torch-diffsim, including the finite element method formulation, time integration scheme, and material model.

Finite Element Method (FEM)#

torch-diffsim uses the finite element method with linear tetrahedral elements to discretize the continuous elastic body.

Mesh Discretization#

The domain \(\Omega \subset \mathbb{R}^3\) is discretized into \(M\) tetrahedral elements:

\[\Omega \approx \bigcup_{e=1}^{M} \Omega_e\]

Each tetrahedron \(\Omega_e\) is defined by 4 vertices with positions \(\mathbf{x}_0, \mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3 \in \mathbb{R}^3\).

Deformation Gradient#

The deformation gradient \(\mathbf{F}\) maps from the reference (rest) configuration to the current (deformed) configuration. For each element, it is computed as:

\[\mathbf{F} = \mathbf{D}_s \mathbf{D}_m^{-1}\]

where:

  • \(\mathbf{D}_m = [\mathbf{X}_1 - \mathbf{X}_0, \mathbf{X}_2 - \mathbf{X}_0, \mathbf{X}_3 - \mathbf{X}_0] \in \mathbb{R}^{3 \times 3}\) contains edge vectors in the rest configuration

  • \(\mathbf{D}_s = [\mathbf{x}_1 - \mathbf{x}_0, \mathbf{x}_2 - \mathbf{x}_0, \mathbf{x}_3 - \mathbf{x}_0] \in \mathbb{R}^{3 \times 3}\) contains edge vectors in the current configuration

For linear tetrahedral elements, \(\mathbf{F}\) is constant within each element (constant strain elements).

Stable Neo-Hookean Material#

The material model defines the relationship between deformation and stress. We use the Stable Neo-Hookean hyperelastic model from Smith et al. (2018).

Energy Density#

The strain energy density \(\Psi: \mathbb{R}^{3 \times 3} \to \mathbb{R}\) is:

\[\Psi(\mathbf{F}) = \frac{\bar\mu}{2}\left[(I_C-3)-\log\left(1+\frac{I_C-3}{4}\right)\right] + \frac{\bar\lambda}{2}(J-1)^2 - \frac{3\bar\mu}{4}(J-1)\]

where:

  • \(I_C = \text{tr}(\mathbf{F}^T \mathbf{F}) = \|\mathbf{F}\|_F^2\) is the first invariant (measures stretch)

  • \(J = \det(\mathbf{F})\) is the Jacobian determinant (measures volume change)

  • \(\mu = \frac{E}{2(1+\nu)}\) is the shear modulus

  • \(\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}\) is Lamé’s first parameter

  • \(\bar\mu = \frac{4}{3}\mu\) and \(\bar\lambda = \lambda + \frac{5}{6}\mu\) are the stabilized parameters

Why “Stable”? This formulation avoids \(\log J\), so its energy and stress remain finite at singular and inverted elements. This does not make explicit integration unconditionally stable; the time step must still resolve the material stiffness.

Total Elastic Energy#

The total elastic energy is obtained by integrating over all elements:

\[E_{\text{elastic}} = \sum_{e=1}^{M} \Psi(\mathbf{F}_e) V_e^0\]

where \(V_e^0 = \frac{1}{6}|\det(\mathbf{D}_m^e)|\) is the rest volume of element \(e\).

Elastic Forces#

The elastic force on each vertex is the negative gradient of the total energy:

\[\mathbf{f}_{\text{elastic}} = -\frac{\partial E_{\text{elastic}}}{\partial \mathbf{x}}\]

This is computed using the first Piola-Kirchhoff stress tensor \(\mathbf{P}\):

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} = \bar\mu\left(1-\frac{1}{I_C+1}\right)\mathbf{F} + \left[\bar\lambda(J-1)-\frac{3\bar\mu}{4}\right] \operatorname{cof}(\mathbf{F})\]

The cofactor form remains defined when \(\mathbf{F}\) is singular.

The force on vertex \(i\) of element \(e\) is:

\[\mathbf{f}_i^e = -V_e^0 \mathbf{P}_e \frac{\partial \mathbf{F}_e}{\partial \mathbf{x}_i}\]

Semi-Implicit (Symplectic Euler) Time Integration#

Time integration advances the simulation forward in time. We use semi-implicit Euler, also known as symplectic Euler.

Integration Scheme#

Given state \((\mathbf{x}^n, \mathbf{v}^n)\) at time \(t_n\), compute the next state:

\[ \begin{align}\begin{aligned}\mathbf{v}^{n+1} &= \mathbf{v}^n + \Delta t \, \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}^n)\\\mathbf{x}^{n+1} &= \mathbf{x}^n + \Delta t \, \mathbf{v}^{n+1}\end{aligned}\end{align} \]

Key property: Velocities are updated using forces at the current position \(\mathbf{x}^n\), then positions are updated using the new velocities \(\mathbf{v}^{n+1}\).

Why Semi-Implicit?#

The scheme is:

  1. Symplectic when unmodified: The damping, velocity limit, constraints, and contact response used by thesolver make its map non-symplectic

  2. First-order accurate: Error is \(O(\Delta t)\)

  3. Conditionally stable: More stable than explicit Euler for stiff systems

  4. Simple: No iterative solver needed (unlike fully implicit methods)

The scheme is called “semi-implicit” because velocities use the current position (explicit) but positions use the new velocity (implicit dependency).

Total Force Computation#

The total force includes multiple components:

\[\mathbf{f}(\mathbf{x},\mathbf{v}) = \mathbf{f}_{\text{elastic}}(\mathbf{x}) + \mathbf{f}_{\text{gravity}} + \mathbf{f}_{\text{contact}}(\mathbf{x},\mathbf{v})\]

Elastic forces: Computed from strain energy as described above

Gravity: \(\mathbf{f}_{\text{gravity}} = \mathbf{M} \mathbf{g}\) where \(\mathbf{g} = [0, -9.8, 0]^T\)

Contact forces: A finite smooth penalty (differentiable solver), projection with friction-like damping (standard solver), or a coupled implicit IPC potential (IPCImplicitEulerSolver)

Damping: \(\mathbf{v} \leftarrow \alpha \mathbf{v}\) with \(\alpha \approx 0.99\)

Substepping#

For stability, each timestep \(\Delta t\) is subdivided into \(n_{\text{sub}}\) substeps:

\[h = \frac{\Delta t}{n_{\text{sub}}}\]

The integration scheme is applied \(n_{\text{sub}}\) times with step size \(h\). Typical values: \(\Delta t = 0.01\), \(n_{\text{sub}} = 4\).

Contact Handling#

We use three approaches depending on the solver:

  • Differentiable solver: A smooth finite penalty to maintain differentiability

  • Standard solver: Ground-plane projection with restitution and ad hoc friction-like tangential damping

  • IPC solver: Frictionless IPC barrier inside a backward-Euler minimization with continuous collision detection on every line-search segment

Barrier Potential#

For ground contact at \(y = 0\), the barrier potential for vertex \(i\) is:

\[\begin{split}b(d_i) = \begin{cases} \frac{\kappa}{3}(\hat{d}-d_i)^3 & \text{if } d_i < \hat{d} \\ 0 & \text{otherwise} \end{cases}\end{split}\]

where:

  • \(d_i = y_i\) is the distance to the ground

  • \(\hat{d}\) is the barrier activation distance (e.g., 0.01 m)

  • \(\kappa\) is the barrier stiffness (e.g., \(10^4\))

Contact Force (Differentiable)#

The normal contact force is the negative gradient of the penalty potential:

\[\mathbf{f}_{\text{contact}} = -\nabla_{\mathbf{x}} \sum_{i=1}^{N} b(d_i)\]

The potential is \(C^2\) at activation and its force is \(C^1\). We also add velocity-dependent tangential damping proportional to the normal force. This is a finite penalty and therefore does not guarantee non-penetration.

Incremental Potential Contact#

For a substep of size \(h\), the frictionless IPC solver computes \(\mathbf{x}^{n+1}\) by minimizing

\[\Pi(\mathbf{x}) = \frac{1}{2}(\mathbf{x}-\hat{\mathbf{x}})^T \mathbf{M}(\mathbf{x}-\hat{\mathbf{x}}) + h^2 E_{\mathrm{elastic}}(\mathbf{x}) + h^2 \kappa B_{\mathrm{IPC}}(\mathbf{x}),\]
\[\hat{\mathbf{x}} = \mathbf{x}^n + h\alpha\mathbf{v}^n + h^2\mathbf{M}^{-1}\mathbf{f}_{\mathrm{ext}}.\]

Here \(\alpha\) is the configured velocity damping factor and \(\mathbf{f}_{\mathrm{ext}}\) includes gravity and user/controller forces.

See the original IPC paper and the official IPC Toolkit simulation guide for the complete formulation.

Ground Projection (Standard)#

For the non-differentiable solver, vertices below the ground are projected back to \(y=0\), normal velocity is adjusted by restitution, and tangential velocity is damped to emulate friction. This preview path does not filter tetrahedral inversion, and its tangential damping is applied once per substep; use IPC when contact feasibility is required. See SemiImplicitSolver._handle_ground_collision in the code.

Mass Matrix#

The mass matrix \(\mathbf{M}\) is diagonal (lumped mass):

\[M_{ii} = \sum_{e \in \text{adj}(i)} \frac{\rho V_e^0}{4}\]

where \(\rho\) is the material density and the sum is over all elements adjacent to vertex \(i\). Each element’s mass is distributed equally to its 4 vertices.

Boundary Conditions#

Fixed vertices: Vertices can be constrained by restoring their prescribed position and setting their velocity to zero:

\[\begin{split}\mathbf{x}_i^{n+1} &= \mathbf{x}_i^{\mathrm{fixed}}, \\ \mathbf{v}_i^{n+1} &= \mathbf{0} \qquad \text{for } i \in \mathcal{F},\end{split}\]

where \(\mathcal{F}\) is the set of fixed vertex indices.

For example, the suspended-bunny demo selects the upper four percent of rest heights before simulation:

support_y = torch.quantile(mesh.vertices[:, 1], 0.96)
fixed = torch.where(mesh.vertices[:, 1] >= support_y)[0]
simulator.set_fixed_vertices(fixed)

The standard and differentiable explicit solvers restore constrained state after each substep. The IPC solver instead removes fixed degrees of freedom from its Newton system. In the differentiable explicit path, the prescribed anchor positions remain differentiable inputs even though dynamic updates cannot move them.

Implementation Notes#

In the code:

  1. \(\mathbf{D}_m^{-1}\) is precomputed and cached for efficiency

  2. Forces are accumulated using index_add_ for parallelism

  3. The explicit solvers use PyTorch tensors and support GPU acceleration

  4. The Stable Neo-Hookean implementation evaluates determinants and cofactors without clamping \(J\)

  5. The IPC solver uses the CPU IPC Toolkit for contact geometry and CCD