explicit: hyperelastic materials on the shared finite-strain material point

Neo-Hookean, Mooney–Rivlin and Yeoh materials in explicit dynamics, in
total form from the deformation gradient (objective by construction):
the one-point hex8 takes its mean F over the undeformed element and the
stress τ/J; tet4 takes the isochoric stress of its own F with the
average nodal pressure of the decoupled volumetric energy (Bonet–Burton
at finite strain). Each element's frequency bound follows the tangent's
largest normal stiffness. Small-displacement steps and thermal strain of
hyperelastic elements are refused by name.

Verified by uniaxial strain to λ = 1.5 of all three laws on hex8 and
tet4 against the material point (to rounding), and the oracle case
explicit-neo-hookean-implicit: a quasi-static explicit pull against the
implicit nonlinear solver (lateral contraction to 0.09%, stress to 0.03%).

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/explicit/src/elements.rs

@@ -54,6 +54,7 @@
 //! reaches its material's erosion plastic strain is removed: it carries no stress from then on,
 //! its nodes keep their mass and momentum, and its internal energy at that moment is reported.
 
+use fieldcraft_fem::constitutive::MaterialModel;
 use fieldcraft_fem::material::stress;
 use fieldcraft_fem::plasticity::{Increment, PointState};
 use fieldcraft_fem::rotation::{self, IDENTITY, Mat3, Polar};
@@ -76,6 +77,10 @@ pub struct Workspace {
     /// volume.
     theta: Vec<f64>,
     tet_volume: Vec<f64>,
+    /// Per node, for hyperelastic tetrahedra: the current and the rest volume of its
+    /// tetrahedra, then their ratio.
+    tet_ratio: Vec<f64>,
+    tet_rest: Vec<f64>,
     /// Per tetrahedron (in [`System::tets`] order): its increment between the two passes.
     tet_increment: Vec<Option<TetIncrement>>,
 }
@@ -84,6 +89,7 @@ impl Workspace {
     pub fn new(sys: &System) -> Workspace {
         let n = sys.labels.len();
         let dofs = sys.dofs();
+        let hyper_tets = sys.tets.iter().any(|t| sys.materials.get(t.material).is_some_and(|m| m.hyperelastic.is_some()));
         Workspace {
             stress_force: vec![0.0; dofs],
             hourglass_force: vec![0.0; dofs],
@@ -91,6 +97,8 @@ impl Workspace {
             contact_force: if sys.contacts.is_empty() { Vec::new() } else { vec![0.0; dofs] },
             theta: vec![0.0; n],
             tet_volume: vec![0.0; n],
+            tet_ratio: if hyper_tets { vec![1.0; n] } else { Vec::new() },
+            tet_rest: if hyper_tets { vec![0.0; n] } else { Vec::new() },
             tet_increment: vec![None; sys.tets.len()],
         }
     }
@@ -314,6 +322,52 @@ fn update(mat: &Material, es: &mut ElementState, q: &Mat3, d: &[f64; 6], vol: f6
     }
 }
 
+/// The hyperelastic update of an element at deformation gradient `f` (the hexahedron's mean,
+/// or the tetrahedron's with its averaged volume ratio): the material point's Kirchhoff stress τ
+/// gives the stress of the nodal forces τ/`j` (`j` the element's own volume ratio, so the forces
+/// are ∫ τ ∇ₓN dV₀), the logarithmic strain is kept for output, and the stress work over the
+/// increment `d` (on the midpoint volume `volume`) is added to the element's energy. Returns the
+/// ratio of the tangent's largest normal stiffness (with the stress) to the material's
+/// small-strain one, at least 1, which scales the element's frequency bound.
+#[allow(clippy::too_many_arguments)]
+fn hyperelastic(
+    model: &MaterialModel,
+    mat: &Material,
+    es: &mut ElementState,
+    f: &Mat3,
+    j: f64,
+    pressure: Option<f64>,
+    d: &[f64; 6],
+    volume: f64,
+) -> Result<f64> {
+    let unstable = |e: fieldcraft_fem::Error| Error::Unstable(e.to_string());
+    if !(j > 0.0 && j.is_finite()) {
+        return Err(Error::Unstable(format!("the material is inverted (volume ratio {j:.3e})")));
+    }
+    let r = model.response(f, None, &[], &mut []).map_err(unstable)?;
+    let before = es.point.stress;
+    es.point.stress = match pressure {
+        // The isochoric part of the element's own F with the averaged pressure.
+        Some(p) => {
+            let iso = model.isochoric(f, None, &[], &mut []).map_err(unstable)?;
+            let mut s = iso.stress.map(|t| t / j);
+            for v in s.iter_mut().take(3) {
+                *v += p;
+            }
+            s
+        }
+        None => r.stress.map(|t| t / j),
+    };
+    let work: f64 = (0..6).map(|k| 0.5 * (before[k] + es.point.stress[k]) * d[k]).sum();
+    es.energy += work * volume;
+    if let Some(e) = fieldcraft_fem::constitutive::logarithmic_strain(f) {
+        es.strain = e;
+    }
+    let normal = (0..3).map(|k| r.tangent[k][k]).fold(0.0, f64::max) / j;
+    let stress = es.point.stress.iter().fold(0.0f64, |m, s| m.max(s.abs()));
+    Ok(((normal + stress) / mat.stiffest).max(1.0))
+}
+
 /// The mechanical part of a strain increment: less the thermal strain increment `thermal` in
 /// every normal direction (isotropic expansion, the same in any frame).
 fn mechanical(mut d: [f64; 6], thermal: f64) -> [f64; 6] {
@@ -493,7 +547,20 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element
     let (d, _) = rotation::split(&g);
     let d = mechanical(d, h.thermal_rate * m.dt);
     let (rho, c) = current(mat, h.kinematics.volume, now.volume);
-    update(mat, es, &q, &d, d[0] + d[1] + d[2], m.dt, rho, mid.volume);
+    let stiffer = match &mat.hyperelastic {
+        Some(model) => {
+            let mut f = h.kinematics.gradient(&e.u);
+            for k in 0..3 {
+                f[k][k] += 1.0;
+            }
+            let j = fieldcraft_fem::dense::det3(&f);
+            hyperelastic(model, mat, es, &f, j, None, &d, mid.volume).map_err(|_| inverted(h.label, m.time))?
+        }
+        None => {
+            update(mat, es, &q, &d, d[0] + d[1] + d[2], m.dt, rho, mid.volume);
+            1.0
+        }
+    };
     if let Some(failure) = es.failed {
         return Ok(Outcome::Eroded(failure));
     }
@@ -526,7 +593,7 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element
             }
         }
     }
-    let element = Element { now: &now, mid: &mid, density: rho, wave_speed: c, length: 2.0 * c / scales.omega2.sqrt() };
+    let element = Element { now: &now, mid: &mid, density: rho, wave_speed: c, length: 2.0 * c / (stiffer * scales.omega2).sqrt() };
     let mut forces = element.finish(&es.point.stress, &e.v, m.bulk);
     forces.hourglass = Some(f);
     Ok(Outcome::Forces(forces))
@@ -554,6 +621,7 @@ fn tetrahedron(
     t: &crate::system::Tet,
     inc: Option<TetIncrement>,
     theta: &[f64],
+    ratio: &[f64],
     m: &Motion,
     es: &mut ElementState,
 ) -> Result<Outcome<4>> {
@@ -564,11 +632,34 @@ fn tetrahedron(
     let thermal = t.thermal_rate * m.dt;
     let d = mechanical(d, thermal);
     let (rho, c) = current(mat, t.kinematics.volume, inc.now.volume);
-    update(mat, es, &q, &d, averaged - 3.0 * thermal, m.dt, rho, inc.mid.volume);
+    let stiffer = match &mat.hyperelastic {
+        // The average nodal pressure of a decoupled material (Bonet & Burton 1998, at finite
+        // strain): each node's pressure U′(J_a) from its tetrahedra's volume ratio, the element's
+        // the mean of its nodes', with the isochoric stress of its own F.
+        Some(model) => {
+            let u = gather(m.u, &t.nodes);
+            let mut f = t.kinematics.gradient(&u);
+            for k in 0..3 {
+                f[k][k] += 1.0;
+            }
+            let j = fieldcraft_fem::dense::det3(&f);
+            let pressure = if model.is_decoupled() {
+                let p: Option<f64> = t.nodes.iter().map(|n| model.volumetric(ratio.get(*n).copied().unwrap_or(1.0), false).map(|(p, _)| p)).sum();
+                p.map(|p| p / 4.0)
+            } else {
+                None
+            };
+            hyperelastic(model, mat, es, &f, j, pressure, &d, inc.mid.volume).map_err(|_| inverted(t.label, m.time))?
+        }
+        None => {
+            update(mat, es, &q, &d, averaged - 3.0 * thermal, m.dt, rho, inc.mid.volume);
+            1.0
+        }
+    };
     if let Some(failure) = es.failed {
         return Ok(Outcome::Eroded(failure));
     }
-    let omega2 = if sys.nonlinear { frequency_bound(&inc.now, inc.now.gradient_norm(), mat.modulus, &t.mass) } else { t.omega2 };
+    let omega2 = stiffer * if sys.nonlinear { frequency_bound(&inc.now, inc.now.gradient_norm(), mat.modulus, &t.mass) } else { t.omega2 };
     let element = Element { now: &inc.now, mid: &inc.mid, density: rho, wave_speed: c, length: 2.0 * c / omega2.sqrt() };
     Ok(Outcome::Forces(element.finish(&es.point.stress, &inc.v, m.bulk)))
 }
@@ -632,16 +723,32 @@ pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Wor
                 *x /= w;
             }
         }
+        if !ws.tet_ratio.is_empty() {
+            ws.tet_ratio.fill(0.0);
+            ws.tet_rest.fill(0.0);
+            for (t, inc) in tets.iter().zip(&ws.tet_increment) {
+                let Some(inc) = inc else { continue };
+                for node in t.nodes {
+                    if let (Some(x), Some(w)) = (ws.tet_ratio.get_mut(node), ws.tet_rest.get_mut(node)) {
+                        *x += inc.now.volume;
+                        *w += t.kinematics.volume;
+                    }
+                }
+            }
+            for (x, w) in ws.tet_ratio.iter_mut().zip(&ws.tet_rest) {
+                *x = if *w > 0.0 { *x / w } else { 1.0 };
+            }
+        }
         // Pass 2: each tetrahedron's update.
         let State { elements, eroded } = state;
-        let (theta, increments) = (&ws.theta, &ws.tet_increment);
+        let (theta, ratio, increments) = (&ws.theta, &ws.tet_ratio, &ws.tet_increment);
         let mut forces: Vec<(usize, Forces<4>)> = Vec::new();
         run(
             elements,
             tets.len(),
             |k| tets.get(k).map_or(usize::MAX, |t| t.index),
             |k, es| match tets.get(k) {
-                Some(t) => tetrahedron(sys, t, increments.get(k).copied().flatten(), theta, m, es),
+                Some(t) => tetrahedron(sys, t, increments.get(k).copied().flatten(), theta, ratio, m, es),
                 None => Ok(Outcome::Skipped),
             },
             |k, outcome| {

crates/explicit/src/system.rs

@@ -3,6 +3,7 @@
 //! ([`crate::discrete`]), prescribed velocities, rigid walls, initial velocities and the load
 //! vector of an explicit dynamic step.
 
+use fieldcraft_fem::constitutive::MaterialModel;
 use fieldcraft_fem::material::Elasticity;
 use fieldcraft_fem::plasticity::J2;
 use fieldcraft_fem::uniform::{Hex8, Tet4, modulus_bound};
@@ -38,6 +39,9 @@ pub struct Material {
     /// The hypoelastic(–plastic) update of an isotropic material; `None` for orthotropic
     /// elasticity, updated with [`Material::elasticity`].
     pub j2: Option<J2>,
+    /// The finite-strain material point of a hyperelastic material (its stress from the
+    /// deformation gradient; `None` for the rate-form materials).
+    pub hyperelastic: Option<MaterialModel>,
     /// The equivalent plastic strain at which elements erode.
     pub erosion: Option<f64>,
     /// The (isotropic) expansion coefficient and reference temperature, when the material
@@ -143,9 +147,6 @@ fn invalid(m: impl Into<String>) -> Error {
 fn material(model: &Model, name: &str) -> Result<Material> {
     let m = model.material(name).ok_or_else(|| invalid(format!("material `{name}` is missing")))?;
     let elastic = m.elastic.as_ref().ok_or_else(|| invalid(format!("material `{name}` has no elastic constants")))?;
-    if let Some(h) = m.hyperelastic {
-        return Err(Error::NotSupported(format!("{} hyperelastic material `{name}` in explicit dynamics", h.name())));
-    }
     let density = m
         .density
         .filter(|d| *d > 0.0 && d.is_finite())
@@ -174,6 +175,12 @@ fn material(model: &Model, name: &str) -> Result<Material> {
         stiffest,
         j2,
         erosion: m.erosion.map(|e| e.plastic_strain),
+        hyperelastic: m
+            .hyperelastic
+            .as_ref()
+            .map(|h| MaterialModel::from_model(elastic, Some(h)))
+            .transpose()
+            .map_err(|e| invalid(format!("material `{name}`: {e}")))?,
     })
 }
 
@@ -368,9 +375,33 @@ impl System {
         sys.load(model, resolved, step)?;
         sys.temperatures(model, resolved, step, settings.end_time)?;
         sys.conditions(resolved, settings)?;
+        sys.check_hyperelastic()?;
         Ok(sys)
     }
 
+    /// Refuse what hyperelastic elements cannot take yet: small-displacement steps (the
+    /// material is a finite-strain law) and thermal strain.
+    fn check_hyperelastic(&self) -> Result<()> {
+        let solids =
+            self.hexes.iter().map(|h| (h.label, h.material, h.thermal_rate)).chain(self.tets.iter().map(|t| (t.label, t.material, t.thermal_rate)));
+        for (label, material, thermal) in solids {
+            let Some(m) = self.materials.get(material).filter(|m| m.hyperelastic.is_some()) else { continue };
+            if !self.nonlinear {
+                return Err(Error::NotSupported(format!(
+                    "hyperelastic material `{}` (element {label}) without nonlinear geometry; hyperelastic laws are finite-strain, set the step's nonlinear_geometry to true",
+                    m.name
+                )));
+            }
+            if thermal != 0.0 {
+                return Err(Error::NotSupported(format!(
+                    "thermal strain in hyperelastic material `{}` (element {label}) in explicit dynamics",
+                    m.name
+                )));
+            }
+        }
+        Ok(())
+    }
+
     /// Add an element's row-sum lumped masses.
     fn lump(&mut self, nodes: &[usize], mass: &[f64]) {
         for (node, v) in nodes.iter().zip(mass) {

crates/explicit/tests/hyperelastic.rs

@@ -0,0 +1,148 @@
+//! Hyperelastic materials in explicit dynamics: the stress of the shared finite-strain material
+//! point (neo-Hookean, Mooney–Rivlin, Yeoh) from the deformation gradient, on one-point hex8 and
+//! F-bar tet4, against the closed form of uniaxial strain reached quasi-statically, and the
+//! refusals.
+//!
+//! In uniaxial strain F = diag(λ, 1, 1): the hyperelastic Cauchy stress is the material point's
+//! Kirchhoff stress over J = λ (fieldcraft_fem::constitutive, whose laws are checked there
+//! against their strain energies), so the explicit elements must reproduce it once the motion has
+//! settled; the wave dynamics of a slow pull leave a fraction of a percent.
+
+use fieldcraft_explicit::{ExplicitResult, Limits, solve_explicit};
+use fieldcraft_fem::constitutive::MaterialModel;
+use fieldcraft_model::*;
+
+const N: usize = 2;
+
+/// A unit cube of N³ hex8 (or six tet4 per cell), held on its faces x = 0, y = 0, y = 1, z = 0,
+/// z = 1 in their normal directions and pulled on x = 1 to a stretch λ over `end`.
+fn block(tet: bool, hyper: Hyperelastic, stretch: f64, end: f64, nonlinear: bool) -> Model {
+    let mut m = Model::default();
+    let h = 1.0 / N as f64;
+    let label = |i: usize, j: usize, k: usize| (1 + i + (N + 1) * (j + (N + 1) * k)) as u32;
+    for k in 0..=N {
+        for j in 0..=N {
+            for i in 0..=N {
+                m.set_node(label(i, j, k), [i as f64 * h, j as f64 * h, k as f64 * h]).unwrap();
+            }
+        }
+    }
+    let corners = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]];
+    let mut next = 1;
+    for k in 0..N {
+        for j in 0..N {
+            for i in 0..N {
+                let at = |c: [usize; 3]| label(i + c[0], j + c[1], k + c[2]);
+                if !tet {
+                    m.set_element(next, ElementKind::Hex8, corners.iter().map(|c| at(*c)).collect()).unwrap();
+                    next += 1;
+                    continue;
+                }
+                for perm in [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]] {
+                    let mut v = [[0usize; 3]; 4];
+                    for s in 0..3 {
+                        v[s + 1] = v[s];
+                        v[s + 1][perm[s]] = 1;
+                    }
+                    let e = |a: usize| [0, 1, 2].map(|d| v[a][d] as f64 - v[0][d] as f64);
+                    let (a, b, c) = (e(1), e(2), e(3));
+                    let vol = a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]);
+                    if vol < 0.0 {
+                        v.swap(1, 2);
+                    }
+                    m.set_element(next, ElementKind::Tet4, v.iter().map(|c| at(*c)).collect()).unwrap();
+                    next += 1;
+                }
+            }
+        }
+    }
+    let all: Vec<u32> = m.elements().keys().copied().collect();
+    m.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
+    let nodes = |f: &dyn Fn([f64; 3]) -> bool| -> Vec<u32> { m.nodes().iter().filter(|(_, n)| f(n.xyz)).map(|(l, _)| *l).collect() };
+    let sets = [
+        ("x0", nodes(&|x| x[0] == 0.0)),
+        ("x1", nodes(&|x| x[0] == 1.0)),
+        ("y", nodes(&|x| x[1] == 0.0 || x[1] == 1.0)),
+        ("z", nodes(&|x| x[2] == 0.0 || x[2] == 1.0)),
+    ];
+    for (name, list) in sets {
+        m.add_set(name, Set { nodes: list, ..Set::default() }).unwrap();
+    }
+    let elastic = Elastic::Isotropic { young: 1.0, poisson: 0.45 };
+    m.add_material("Rubber", Material { elastic: Some(elastic), density: Some(1.0), hyperelastic: Some(hyper), ..Material::default() }).unwrap();
+    let integration = if tet { Integration::Full } else { Integration::Reduced };
+    let formulation = Formulation { integration, ..Formulation::default() };
+    m.add_section("all", Section { set: "all".into(), material: "Rubber".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
+    let held = |set: &str, k: usize, v: f64| {
+        let mut values = [None; 6];
+        values[k] = Some(v);
+        BoundaryCondition::Displacement { target: Target::Set(set.into()), values }
+    };
+    let mut step = ExplicitDynamic::new(end);
+    step.nonlinear_geometry = nonlinear;
+    step.output_interval = Some(end / 2.0);
+    let mut s = Step::new("Pull", Procedure::ExplicitDynamic(step));
+    for (name, bc) in [("x0", held("x0", 0, 0.0)), ("x1", held("x1", 0, stretch - 1.0)), ("y", held("y", 1, 0.0)), ("z", held("z", 2, 0.0))] {
+        m.add_boundary_condition(name, bc).unwrap();
+        s.boundary_conditions.insert(name.into());
+    }
+    m.add_step(s, None).unwrap();
+    m
+}
+
+/// σ of uniaxial strain λ from the material point.
+fn closed_form(hyper: Hyperelastic, stretch: f64) -> [f64; 6] {
+    let model = MaterialModel::from_model(&Elastic::Isotropic { young: 1.0, poisson: 0.45 }, Some(&hyper)).unwrap();
+    let f = [[stretch, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
+    model.response(&f, None, &[], &mut []).unwrap().stress.map(|t| t / stretch)
+}
+
+fn mean_stress(r: &ExplicitResult) -> [f64; 6] {
+    let f = r.frames.last().unwrap();
+    let mut s = [0.0; 6];
+    for e in &f.stress {
+        for k in 0..6 {
+            s[k] += e[k] / f.stress.len() as f64;
+        }
+    }
+    s
+}
+
+#[test]
+fn hyperelastic_uniaxial_strain_matches_the_material_point() {
+    let stretch = 1.5;
+    // Dilatational wave speed √((κ + 4μ/3)/ρ) ≈ 2.0: a pull over 400 transits of the cube.
+    let end = 200.0;
+    let materials = [Hyperelastic::NeoHookean {}, Hyperelastic::MooneyRivlin { c01: 0.08 }, Hyperelastic::Yeoh { c20: -0.01, c30: 0.002 }];
+    for hyper in materials {
+        let reference = closed_form(hyper, stretch);
+        for tet in [false, true] {
+            let r = solve_explicit(&block(tet, hyper, stretch, end, true), None, &Limits::default()).unwrap();
+            assert!(r.completed);
+            let s = mean_stress(&r);
+            let spread = r.frames.last().unwrap().stress.iter().map(|e| (e[0] - s[0]).abs()).fold(0.0, f64::max);
+            eprintln!(
+                "{} {}: σ_xx {:.5} vs {:.5}, σ_yy {:.5} vs {:.5}; spread {spread:.2e}; balance {:.2e}; {} increments",
+                hyper.name(),
+                if tet { "tet4" } else { "hex8" },
+                s[0],
+                reference[0],
+                s[1],
+                reference[1],
+                r.energy.relative_error(r.initial_kinetic),
+                r.increments
+            );
+            for k in 0..3 {
+                assert!((s[k] - reference[k]).abs() < 0.01 * reference[0].abs(), "{} tet {tet}: {s:?} vs {reference:?}", hyper.name());
+            }
+            assert!(spread < 0.01 * reference[0].abs(), "uniform state: {spread}");
+            assert!(r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy);
+        }
+    }
+}
+
+#[test]
+fn hyperelastic_materials_need_nonlinear_geometry() {
+    let e = solve_explicit(&block(false, Hyperelastic::NeoHookean {}, 1.1, 1.0, false), None, &Limits::default()).unwrap_err().to_string();
+    assert!(e.contains("nonlinear geometry") && e.contains("Rubber"), "{e}");
+}