"""
Material models for finite element simulation
This module implements hyperelastic material models for FEM simulation.
The primary model is the Stable Neo-Hookean formulation, which provides
numerical stability even under large deformations and element inversion.
The strain energy density is derived from the deformation gradient :math:`\\mathbf{F}`,
and forces are computed as the negative gradient of the total elastic energy.
"""
import math
import torch
[docs]
class StableNeoHookean:
"""
Stable Neo-Hookean hyperelastic material model
This class implements the stable Neo-Hookean formulation from Smith et al.
(2018), "Stable Neo-Hookean Flesh Simulation". Up to an additive constant, its
energy density is:
.. math::
\\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:
- :math:`\\mathbf{F}` is the deformation gradient tensor (3×3)
- :math:`I_C = \\text{tr}(\\mathbf{F}^T \\mathbf{F}) = \\|\\mathbf{F}\\|_F^2` is the first invariant
- :math:`J = \\det(\\mathbf{F})` is the Jacobian determinant (volume ratio)
- :math:`\\mu = \\frac{E}{2(1+\\nu)}` is the shear modulus
- :math:`\\lambda = \\frac{E\\nu}{(1+\\nu)(1-2\\nu)}` is Lamé's first parameter
Here :math:`\\bar\\mu=4\\mu/3` and
:math:`\\bar\\lambda=\\lambda+5\\mu/6`. This formulation remains finite for
singular and inverted elements; explicit time integration still requires a
sufficiently small time step.
Parameters:
youngs_modulus (float): Young's modulus :math:`E` in Pascals (default: 1e6)
poissons_ratio (float): Poisson's ratio :math:`\\nu` (default: 0.45)
Attributes:
E (float): Young's modulus
nu (float): Poisson's ratio
mu (float): Shear modulus (Lamé's second parameter)
lam (float): Lamé's first parameter
Reference:
Smith, B., De Goes, F., & Kim, T. (2018). Stable neo-hookean flesh simulation.
ACM Transactions on Graphics (TOG), 37(2), 1-15.
"""
[docs]
def __init__(self, youngs_modulus=1e6, poissons_ratio=0.45):
"""
Initialize material with elastic constants
Args:
youngs_modulus: Young's modulus (E)
poissons_ratio: Poisson's ratio (ν)
"""
if not math.isfinite(float(youngs_modulus)) or youngs_modulus < 0:
raise ValueError("youngs_modulus must be finite and nonnegative")
if not math.isfinite(float(poissons_ratio)) or not (
-1.0 < poissons_ratio < 0.5
):
raise ValueError("poissons_ratio must satisfy -1 < nu < 0.5")
self.E = youngs_modulus
self.nu = poissons_ratio
# Convert to Lamé parameters
# μ = E / (2(1+ν))
# λ = E*ν / ((1+ν)(1-2ν))
self.mu = self.E / (2.0 * (1.0 + self.nu))
self.lam = self.E * self.nu / ((1.0 + self.nu) * (1.0 - 2.0 * self.nu))
[docs]
def energy_density(self, F):
"""
Compute strain energy density for deformation gradient F
Args:
F: :math:`(M, 3, 3)` deformation gradient tensor
Returns:
psi: :math:`(M,)` energy density for each element
"""
# Compute first invariant Ic = trace(F^T F) = ||F||_F^2
Ic = torch.sum(F * F, dim=(1, 2)) # (M,)
J, _ = self._determinant_and_cofactor(F)
stable_mu = 4.0 * self.mu / 3.0
stable_lam = self.lam + 5.0 * self.mu / 6.0
Ic_offset = Ic - 3.0
psi = (
stable_mu / 2.0 * (Ic_offset - torch.log1p(Ic_offset / 4.0))
+ stable_lam / 2.0 * (J - 1.0) ** 2
- 3.0 * stable_mu / 4.0 * (J - 1.0)
)
return psi
@staticmethod
def _determinant_and_cofactor(F):
"""Compute determinant and cofactor without a matrix inverse."""
a, b, c = F.unbind(dim=2)
cofactor = torch.stack(
(
torch.cross(b, c, dim=1),
torch.cross(c, a, dim=1),
torch.cross(a, b, dim=1),
),
dim=2,
)
J = torch.sum(a * cofactor[:, :, 0], dim=1)
return J, cofactor
[docs]
def first_piola_kirchhoff_stress(self, F):
"""
Compute first Piola-Kirchhoff stress tensor P = ∂Ψ/∂F
Args:
F: :math:`(M, 3, 3)` deformation gradient tensor
Returns:
P: :math:`(M, 3, 3)` first Piola-Kirchhoff stress tensor
"""
return self.first_piola_kirchhoff_stress_analytic(F)
[docs]
def first_piola_kirchhoff_stress_analytic(self, F):
"""
Compute the exact first Piola-Kirchhoff stress tensor analytically.
Args:
F: :math:`(M, 3, 3)` deformation gradient tensor
Returns:
P: :math:`(M, 3, 3)` first Piola-Kirchhoff stress tensor
"""
Ic = torch.sum(F * F, dim=(1, 2))
J, cofactor = self._determinant_and_cofactor(F)
stable_mu = 4.0 * self.mu / 3.0
stable_lam = self.lam + 5.0 * self.mu / 6.0
return (
stable_mu * (1.0 - 1.0 / (Ic + 1.0)).view(-1, 1, 1) * F
+ (stable_lam * (J - 1.0) - 3.0 * stable_mu / 4.0).view(-1, 1, 1) * cofactor
)
[docs]
def compute_elastic_forces(self, F, Dm_inv, volume):
"""
Compute elastic forces from stress tensor
Args:
F: :math:`(M, 3, 3)` deformation gradient
Dm_inv: :math:`(M, 3, 3)` inverse rest shape matrix
volume: :math:`(M,)` rest volume of each element
Returns:
forces: :math:`(M, 4, 3)` forces on vertices of each element
"""
# Compute stress tensor
P = self.first_piola_kirchhoff_stress_analytic(F) # (M, 3, 3)
# H = -volume * P * Dm_inv^T (force matrix)
H = -volume.unsqueeze(-1).unsqueeze(-1) * torch.bmm(
P, Dm_inv.transpose(1, 2)
) # (M, 3, 3)
# Extract forces for each vertex
# f1, f2, f3 are columns of H
# f0 = -(f1 + f2 + f3) (force balance)
f1 = H[:, :, 0] # (M, 3)
f2 = H[:, :, 1] # (M, 3)
f3 = H[:, :, 2] # (M, 3)
f0 = -(f1 + f2 + f3) # (M, 3)
forces = torch.stack([f0, f1, f2, f3], dim=1) # (M, 4, 3)
return forces