diff --git a/ROADMAP.md b/ROADMAP.md
index cc908f2..4cb66c1 100644
--- a/ROADMAP.md
+++ b/ROADMAP.md
@@ -28,7 +28,7 @@ hex meshes today and it gives the best elements where stresses are read.
 | M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, springs, point masses; amplitudes, multiple steps; nonlinear static steps in decks | planned |
 | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; heat transfer | planned |
 | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned |
-| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `modelcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (objective rates, J2, Johnson–Cook, damage, erosion); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started |
+| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `modelcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started |
 | M11 Equations of state and shells | Polynomial, ideal-gas, Mie–Grüneisen, JWL with programmed burn, tabulated; energy coupling; shock verification (Rankine–Hugoniot, Sod, Noh, Sedov); Belytschko–Lin–Tsay shells with self-contact on quad meshes | planned |
 | M12 CFD: lattice Boltzmann | D3Q19/D3Q27 on grids voxelised from the face-tagged surface; TRT, regularised and cumulant collision; interpolated walls, inlets and outlets; LES; octree refinement shared with the hex mesher; wgpu backend (deterministic per device, checked against the CPU) (CPU core: done — `modelcraft-lbm` with D3Q19 TRT and Guo forcing, halfway bounce-back and moving walls, slip and periodic faces, Zou–He velocity and pressure faces, exact voxelisation of obstacle bodies, momentum-exchange forces per face; fluid materials, fluid domains and fluid steps in the model; `model.fluid_domain.*`, `solve.fluid`; verified by Poiseuille channel and pipe, Couette, Taylor–Green decay, sphere-array Stokes drag and mass conservation; Schäfer–Turek 2D-1 as a benchmark) | started |
 | M13 Fluid–structure interaction | Immersed-boundary coupling of the lattice Boltzmann fluid with the explicit or implicit structure; interface force transfer; sub-cycling; Aitken then interface quasi-Newton coupling; free-surface flow | planned |
diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs
index 88a7252..b1c10fd 100644
--- a/crates/decks/src/abaqus.rs
+++ b/crates/decks/src/abaqus.rs
@@ -82,6 +82,21 @@
 //! materials outside solid sections (shell and beam material axes need orientations), and
 //! distributing couplings that leave a translation free.
 //!
+//! # Materials
+//!
+//! Isotropic elasticity is `*ELASTIC, TYPE=ISOTROPIC`, engineering constants
+//! `TYPE=ENGINEERING CONSTANTS`, density `*DENSITY`. Von Mises plasticity is `*PLASTIC`: a
+//! hardening table as its (yield stress, plastic strain) lines; bilinear isotropic hardening as
+//! two lines, the second at a plastic strain of [`BILINEAR_REACH`] (Abaqus holds the last yield
+//! stress beyond a table's end, so the line is written out further than any analysis strains);
+//! bilinear kinematic hardening as `HARDENING=KINEMATIC` with its yield stress at plastic strains
+//! 0 and 1 (Abaqus's linear kinematic model extends that line); Johnson–Cook as
+//! `HARDENING=JOHNSON COOK` (A, B, n, m, melting and transition temperature, the transition
+//! temperature being the model's room temperature) with `*RATE DEPENDENT, TYPE=JOHNSON COOK`
+//! (C, ε̇₀) when C is not zero, and adiabatic heating as `*SPECIFIC HEAT` and
+//! `*INELASTIC HEAT FRACTION`. Damage and erosion are failure models of Abaqus/Explicit, and this
+//! writer writes Abaqus/Standard procedures: they are reported, not written.
+//!
 //! Node order: the model stores elements in Abaqus order, so connectivity passes through; the
 //! model's face numbers are Abaqus's S1–S6 for solids, and SPOS/SNEG for shell faces 1/2.
 
@@ -89,14 +104,17 @@ use std::collections::{BTreeMap, BTreeSet};
 use std::fmt::Write as _;
 
 use modelcraft_model::{
-    BeamProfile, BoundaryCondition, CouplingKind, Elastic, Element, ElementKind, Integration, Load, Model, OutputVariable, Procedure, Resolved,
-    Section, SectionKind, Step, Symmetry, Target,
+    BeamProfile, BoundaryCondition, CouplingKind, Elastic, Element, ElementKind, Hardening, Integration, Load, Model, OutputVariable, Plastic,
+    Procedure, Resolved, Section, SectionKind, Step, Symmetry, Target,
 };
 
 use crate::{Error, Result, unsupported};
 
 const FORMAT: &str = "Abaqus";
 
+/// The plastic strain a bilinear isotropic hardening line is written out to (1000%).
+pub const BILINEAR_REACH: f64 = 10.0;
+
 /// Longest line Abaqus reads.
 pub const MAX_LINE: usize = 256;
 
@@ -251,6 +269,13 @@ fn check_expressible(model: &Model) -> Result<()> {
     for step in model.steps().iter().filter(|s| matches!(s.procedure, Procedure::ExplicitDynamic(_))) {
         bad.insert(format!("step `{}`: explicit dynamic steps (not written to Abaqus input yet)", step.name));
     }
+    for (name, m) in model.materials() {
+        if m.damage.is_some() || m.erosion.is_some() {
+            bad.insert(format!(
+                "material `{name}`: damage and erosion (Abaqus/Explicit failure models; this writer writes Abaqus/Standard procedures so far)"
+            ));
+        }
+    }
     for (name, c) in model.couplings() {
         if c.kind == CouplingKind::Distributing && !(1..=3).all(|d| c.dofs.contains(d)) {
             bad.insert(format!("coupling `{name}`: an Abaqus distributing coupling always couples translations 1-3"));
@@ -621,6 +646,50 @@ impl Deck {
                 self.line("*DENSITY");
                 self.numbers_line("", &[rho])?;
             }
+            if let Some(p) = &m.plastic {
+                self.plastic(p)?;
+            }
+        }
+        Ok(())
+    }
+
+    /// Von Mises plasticity (see the module documentation for the mapping).
+    fn plastic(&mut self, p: &Plastic) -> Result<()> {
+        match p {
+            Plastic::Bilinear { yield_stress, hardening_modulus, hardening } => {
+                if *hardening_modulus == 0.0 {
+                    self.line("*PLASTIC");
+                    self.numbers_line("", &[*yield_stress, 0.0])?;
+                } else if *hardening == Hardening::Kinematic {
+                    self.line("*PLASTIC, HARDENING=KINEMATIC");
+                    self.numbers_line("", &[*yield_stress, 0.0])?;
+                    self.numbers_line("", &[yield_stress + hardening_modulus, 1.0])?;
+                } else {
+                    self.line("*PLASTIC");
+                    self.numbers_line("", &[*yield_stress, 0.0])?;
+                    self.numbers_line("", &[yield_stress + hardening_modulus * BILINEAR_REACH, BILINEAR_REACH])?;
+                }
+            }
+            Plastic::Table { points } => {
+                self.line("*PLASTIC");
+                for point in points {
+                    self.numbers_line("", &[point.stress, point.strain])?;
+                }
+            }
+            Plastic::JohnsonCook { a, b, n, c, m, reference_strain_rate, room_temperature, melt_temperature, specific_heat, taylor_quinney } => {
+                self.line("*PLASTIC, HARDENING=JOHNSON COOK");
+                self.numbers_line("", &[*a, *b, *n, *m, *melt_temperature, *room_temperature])?;
+                if *c != 0.0 {
+                    self.line("*RATE DEPENDENT, TYPE=JOHNSON COOK");
+                    self.numbers_line("", &[*c, *reference_strain_rate])?;
+                }
+                if let Some(cp) = specific_heat {
+                    self.line("*SPECIFIC HEAT");
+                    self.numbers_line("", &[*cp])?;
+                    self.line("*INELASTIC HEAT FRACTION");
+                    self.numbers_line("", &[*taylor_quinney])?;
+                }
+            }
         }
         Ok(())
     }
diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs
index 8721b34..8e6cd6a 100644
--- a/crates/decks/src/lsdyna/mod.rs
+++ b/crates/decks/src/lsdyna/mod.rs
@@ -8,7 +8,13 @@
 //!   `*SECTION_BEAM` (rods as trusses) and the element formulation mapped from the section's
 //!   ([`elements`]); `*HOURGLASS` where a reduced-integration formulation needs it;
 //! - `*MAT_ELASTIC` for isotropic and `*MAT_ORTHOTROPIC_ELASTIC` (global material axes) for
-//!   orthotropic materials;
+//!   orthotropic materials; von Mises plasticity as `*MAT_PLASTIC_KINEMATIC` (bilinear: the
+//!   tangent modulus E·H/(E + H) of the uniaxial curve, β = 1 isotropic or 0 kinematic, the
+//!   erosion plastic strain as its failure strain) or `*MAT_PIECEWISE_LINEAR_PLASTICITY` (a
+//!   hardening table as a `*DEFINE_CURVE` of yield stress against plastic strain, the erosion
+//!   plastic strain as its failure strain). Johnson–Cook plasticity and damage are reported: their
+//!   LS-DYNA material needs an equation of state and belongs to explicit analyses, which this
+//!   writer does not write yet;
 //! - `*ELEMENT_SOLID` (tetrahedra and wedges as degenerate hexahedra, 10-node tetrahedra),
 //!   `*ELEMENT_SOLID_H20`, `*ELEMENT_SHELL` and `*ELEMENT_BEAM`;
 //! - `*SET_NODE_LIST` and `*SET_SEGMENT` for the sets boundary conditions, loads and couplings
@@ -47,8 +53,8 @@ mod elements;
 use std::collections::{BTreeMap, BTreeSet};
 
 use modelcraft_model::{
-    BoundaryCondition, CouplingKind, Diagnostic, Dofs, Elastic, Element, ElementKind, Load, Model, OutputVariable, Procedure, Resolved, Section,
-    SectionKind, Step, Target,
+    BoundaryCondition, CouplingKind, Diagnostic, Dofs, Elastic, Element, ElementKind, Hardening, Load, Material, Model, OutputVariable, Plastic,
+    Procedure, Resolved, Section, SectionKind, Step, Target,
 };
 
 use crate::{Error, Result};
@@ -117,6 +123,8 @@ struct Plan<'a> {
     uniform: BTreeMap<&'a str, ([f64; 3], u32)>,
     /// The part set gravity acts on, when it acts on a set.
     gravity_parts: Option<Vec<u32>>,
+    /// Material name → the curve id of its hardening table.
+    curves: BTreeMap<&'a str, u32>,
 }
 
 /// Write the model as an LS-DYNA keyword deck.
@@ -304,13 +312,25 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result<Plan<'a>> {
             None
         }
     };
-    for m in materials.keys() {
-        if model.material(m).is_some_and(|m| m.density.is_none()) {
-            refused.insert(format!("material `{m}` has no density, and LS-DYNA requires a mass density on every material"));
+    let mut curves = BTreeMap::new();
+    for name in materials.keys() {
+        let Some(m) = model.material(name) else { continue };
+        if m.density.is_none() {
+            refused.insert(format!("material `{name}` has no density, and LS-DYNA requires a mass density on every material"));
+        }
+        if matches!(m.plastic, Some(Plastic::JohnsonCook { .. })) || m.damage.is_some() {
+            refused.insert(format!(
+                "material `{name}`: Johnson–Cook plasticity and damage (LS-DYNA's Johnson–Cook material needs an equation of state and is for explicit analyses, which this writer does not write yet)"
+            ));
+        }
+        if matches!(m.plastic, Some(Plastic::Table { .. })) {
+            // Curve ids after the steps' load ramps (ids 1 to the number of steps).
+            let lcid = u32::try_from(model.steps().len() + 1 + curves.len()).map_err(|_| Error::Limit("too many curves".into()))?;
+            curves.insert(*name, lcid);
         }
     }
     crate::unsupported(FORMAT, refused)?;
-    Ok(Plan { parts, materials, hourglasses, sets, node_sets, segment_sets, uniform, gravity_parts })
+    Ok(Plan { parts, materials, hourglasses, sets, node_sets, segment_sets, uniform, gravity_parts, curves })
 }
 
 /// A 1-based id from a 0-based index.
@@ -374,6 +394,7 @@ fn materials(d: &mut Deck, model: &Model, plan: &Plan) -> Result<()> {
         let m = model.material(name).ok_or_else(|| broken(format!("material `{name}`")))?;
         let rho = m.density.unwrap_or_default();
         match m.elastic {
+            Some(Elastic::Isotropic { young, poisson }) if m.plastic.is_some() => plastic_material(d, name, mid, rho, young, poisson, m, plan)?,
             Some(Elastic::Isotropic { young, poisson }) => {
                 d.line("*MAT_ELASTIC_TITLE");
                 d.line(name);
@@ -401,6 +422,56 @@ fn materials(d: &mut Deck, model: &Model, plan: &Plan) -> Result<()> {
     Ok(())
 }
 
+/// The plastic materials: `*MAT_PLASTIC_KINEMATIC` or `*MAT_PIECEWISE_LINEAR_PLASTICITY` with its
+/// curve, in place of `*MAT_ELASTIC`.
+fn plastic_material(d: &mut Deck, name: &str, mid: u32, rho: f64, young: f64, poisson: f64, m: &Material, plan: &Plan) -> Result<()> {
+    let fail = m.erosion.map_or(Field::Blank, |e| Field::from(e.plastic_strain));
+    match &m.plastic {
+        Some(Plastic::Bilinear { yield_stress, hardening_modulus, hardening }) => {
+            d.line("*MAT_PLASTIC_KINEMATIC_TITLE");
+            d.line(name);
+            d.header(&["mid", "ro", "e", "pr", "sigy", "etan", "beta"]);
+            let tangent = young * hardening_modulus / (young + hardening_modulus);
+            let beta = if *hardening == Hardening::Kinematic { 0.0 } else { 1.0 };
+            d.card([
+                Field::from(mid),
+                Field::from(rho),
+                Field::from(young),
+                Field::from(poisson),
+                Field::from(*yield_stress),
+                Field::from(tangent),
+                Field::from(beta),
+            ]);
+            d.header(&["src", "srp", "fs"]);
+            d.card([Field::Blank, Field::Blank, fail]);
+        }
+        Some(Plastic::Table { points }) => {
+            let lcid = plan.curves.get(name).copied().ok_or_else(|| broken(format!("material `{name}` has no curve")))?;
+            d.line("*MAT_PIECEWISE_LINEAR_PLASTICITY_TITLE");
+            d.line(name);
+            d.header(&["mid", "ro", "e", "pr", "sigy", "etan", "fail"]);
+            let first = m.plastic.as_ref().map_or(0.0, Plastic::initial_yield);
+            d.card([Field::from(mid), Field::from(rho), Field::from(young), Field::from(poisson), Field::from(first), Field::Blank, fail]);
+            d.header(&["c", "p", "lcss"]);
+            d.card([Field::Blank, Field::Blank, Field::from(lcid)]);
+            d.header(&["eps1"]);
+            d.card([Field::Blank]);
+            d.header(&["es1"]);
+            d.card([Field::Blank]);
+            d.line("*DEFINE_CURVE_TITLE");
+            d.line(&format!("{name} hardening"));
+            d.header(&["lcid"]);
+            d.card([Field::from(lcid)]);
+            d.header(&["a1", "o1"]);
+            for p in points {
+                d.card([Field::from(p.strain), Field::from(p.stress)]);
+            }
+        }
+        _ => return Err(broken(format!("material `{name}`: a plasticity LS-DYNA cannot take"))),
+    }
+    Ok(())
+}
+
 fn parts(d: &mut Deck, plan: &Plan) -> Result<()> {
     for p in &plan.parts {
         match p.section.kind {
diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs
index be07584..297f2a7 100644
--- a/crates/decks/src/nastran/mod.rs
+++ b/crates/decks/src/nastran/mod.rs
@@ -14,9 +14,12 @@
 //! - the element formulation: `PSOLID`'s integration network and scheme (see [`solid_rule`]);
 //!   formulations Nastran's linear elements have no counterpart for are reported, not
 //!   approximated;
-//! - isotropic materials: `MAT1`; orthotropic engineering constants: `MAT9` (the stiffness matrix
-//!   in the basic axes) for solids and `MAT8` for shells lying in planes of constant Z, with the
-//!   material X axis tied to the basic X axis through a `CORD2R`;
+//! - isotropic materials: `MAT1`, with `MATS1` for von Mises plasticity (bilinear by its
+//!   hardening slope, a hardening table as a `TABLES1` stress–strain curve; Johnson–Cook, damage
+//!   and erosion are reported, not written); a linear solution reads the elastic entry alone;
+//! - orthotropic engineering constants: `MAT9` (the stiffness matrix in the basic axes) for
+//!   solids and `MAT8` for shells lying in planes of constant Z, with the material X axis tied to
+//!   the basic X axis through a `CORD2R`;
 //! - boundary conditions: per step, one constraint set of `SPC1` (held degrees of freedom) and
 //!   `SPC` (prescribed values); in a modes step every constrained degree of freedom is held;
 //! - forces, moments and tractions: `FORCE` and `MOMENT` at the nodal values the model resolves,
@@ -44,8 +47,8 @@ use std::collections::{BTreeMap, BTreeSet};
 use std::fmt::Write as _;
 
 use modelcraft_model::{
-    BeamProfile, Coupling, CouplingKind, Dofs, Elastic, ElementKind, Formulation, Integration, Load, MAX_LABEL, Model, OutputVariable, Procedure,
-    Resolved, SectionKind, Step,
+    BeamProfile, Coupling, CouplingKind, Dofs, Elastic, ElementKind, Formulation, Hardening, Integration, Load, MAX_LABEL, Model, OutputVariable,
+    Plastic, Procedure, Resolved, SectionKind, Step,
 };
 
 use crate::{Error, Result};
@@ -184,6 +187,8 @@ struct MaterialIds {
     solid: Option<u32>,
     /// `MAT8`, for orthotropic materials of shells.
     shell: Option<u32>,
+    /// `TABLES1` of a tabulated plasticity's stress–strain curve.
+    table: Option<u32>,
 }
 
 /// Which element families use a material.
@@ -307,7 +312,7 @@ impl<'m> Deck<'m> {
     }
 
     fn plan(&mut self) -> Result<()> {
-        self.plan_materials();
+        self.plan_materials()?;
         self.plan_sections()?;
         self.plan_couplings()?;
         self.plan_constraints()?;
@@ -321,7 +326,7 @@ impl<'m> Deck<'m> {
     }
 
     /// Material ids, for the forms the sections need.
-    fn plan_materials(&mut self) {
+    fn plan_materials(&mut self) -> Result<()> {
         let model = self.model;
         let mut uses: BTreeMap<&str, MaterialUse> = BTreeMap::new();
         for s in model.sections().values() {
@@ -353,8 +358,24 @@ impl<'m> Deck<'m> {
                 }
                 None => MaterialIds::default(),
             };
+            let mut ids = ids;
+            if let Some(m) = model.material(name) {
+                match &m.plastic {
+                    Some(Plastic::JohnsonCook { .. }) => {
+                        self.problems
+                            .insert(format!("material `{name}`: Johnson–Cook plasticity (MATS1 has no rate- and temperature-dependent flow stress)"));
+                    }
+                    Some(Plastic::Table { .. }) => ids.table = Some(next_id(&mut self.next_sid)?),
+                    _ => {}
+                }
+                if m.damage.is_some() || m.erosion.is_some() {
+                    self.problems
+                        .insert(format!("material `{name}`: damage and erosion (Nastran's linear and implicit solutions have no element failure)"));
+                }
+            }
             self.materials.insert(name, ids);
         }
+        Ok(())
     }
 
     /// A new property for `elements`.
@@ -706,6 +727,9 @@ impl<'m> Deck<'m> {
             }
             if let (Some(card), Some(mid)) = (card, ids.isotropic) {
                 write_named(out, &format!("material {name}"), mid, card)?;
+                if let (Some(plastic), Some(Elastic::Isotropic { young, .. })) = (&m.plastic, m.elastic) {
+                    plasticity(out, name, mid, young, plastic, ids.table)?;
+                }
             }
         }
 
@@ -907,6 +931,29 @@ fn beam_dimensions(p: BeamProfile) -> (&'static str, Vec<f64>) {
     }
 }
 
+/// `MATS1` (von Mises yield, `YF` 1) of an isotropic material's plasticity: a bilinear law by its
+/// work-hardening slope `H` (stress against plastic strain) and initial yield stress `LIMIT1`,
+/// isotropic (`HR` 1) or kinematic (`HR` 2); a hardening table as a `TABLES1` stress–strain
+/// curve through the origin, with the total strain σ/E + ε_p of each point.
+fn plasticity(out: &mut String, name: &str, mid: u32, young: f64, plastic: &Plastic, table: Option<u32>) -> Result<()> {
+    let what = format!("plasticity of {name}");
+    match (plastic, table) {
+        (Plastic::Bilinear { yield_stress, hardening_modulus, hardening }, _) => {
+            let hr = if *hardening == Hardening::Kinematic { 2 } else { 1 };
+            let card = Card::new("MATS1").int(mid).blank().text("PLASTIC").real(*hardening_modulus).int(1).int(hr).real(*yield_stress);
+            write_named(out, &what, mid, card)
+        }
+        (Plastic::Table { points }, Some(tid)) => {
+            write_named(out, &what, mid, Card::new("MATS1").int(mid).int(tid).text("PLASTIC").blank().int(1).int(1))?;
+            let curve = std::iter::once([0.0, 0.0]).chain(points.iter().map(|p| [p.stress / young + p.strain, p.stress])).flatten();
+            write_named(out, &format!("stress-strain curve of {name}"), tid, Card::new("TABLES1").int(tid).end_line().reals(curve).text("ENDT"))
+        }
+        (Plastic::Table { .. }, None) | (Plastic::JohnsonCook { .. }, _) => {
+            Err(Error::Model(modelcraft_model::Error::Invalid(format!("material `{name}`: its plasticity was not planned for MATS1"))))
+        }
+    }
+}
+
 /// An entry preceded by a comment naming what it holds (`material Steel: MAT1 1`).
 fn write_named(out: &mut String, what: &str, id: u32, card: Card) -> Result<()> {
     comment(out, &format!("{what}: {} {id}", card.name()));
diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs
index 3e4a380..2cc6ad0 100644
--- a/crates/decks/tests/abaqus.rs
+++ b/crates/decks/tests/abaqus.rs
@@ -15,7 +15,7 @@ use modelcraft_model::{
 // Models
 
 fn steel() -> Material {
-    Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None }
+    Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() }
 }
 
 fn faces(list: &[(u32, u8)]) -> Set {
@@ -91,7 +91,7 @@ fn plate() -> Model {
     m.add_set("top", faces(&[(1, 1), (2, 1), (3, 1), (4, 1)])).unwrap();
     m.add_set("edge", Set { nodes: vec![1, 4, 7], ..Set::default() }).unwrap();
     m.add_set("mid", Set { nodes: vec![2, 5, 8], ..Set::default() }).unwrap();
-    let aluminium = Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9), fluid: None };
+    let aluminium = Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9), ..Material::default() };
     m.add_material("Aluminium", aluminium).unwrap();
     let kind = SectionKind::Shell { thickness: 2.0, offset: 0.5, points: 5 };
     let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() };
@@ -523,8 +523,11 @@ fn largest_labels_longest_numbers_and_names_fit_a_line() {
     m.set_element(modelcraft_model::MAX_LABEL, ElementKind::Hex20, (0..20).map(|i| base + i).collect()).unwrap();
     let long = format!("A{}", "x".repeat(modelcraft_model::MAX_NAME_LEN - 1));
     m.add_set(&long, Set { elements: vec![modelcraft_model::MAX_LABEL], nodes: (0..20).map(|i| base + i).collect(), ..Set::default() }).unwrap();
-    m.add_material(&long, Material { elastic: Some(Elastic::Isotropic { young: ugly.abs(), poisson: 0.3 }), density: Some(1.0), fluid: None })
-        .unwrap();
+    m.add_material(
+        &long,
+        Material { elastic: Some(Elastic::Isotropic { young: ugly.abs(), poisson: 0.3 }), density: Some(1.0), ..Material::default() },
+    )
+    .unwrap();
     m.add_section(&long, Section { set: long.clone(), material: long.clone(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
         .unwrap();
     m.add_reference_point(&long, ReferencePoint { xyz: [ugly; 3], label: 0 }).unwrap();
@@ -640,7 +643,7 @@ fn what_abaqus_cannot_express_is_named() {
 
     let mut ortho = plate();
     let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] };
-    ortho.replace_material("Aluminium", Material { elastic: Some(carbon), density: Some(1.6e-9), fluid: None }).unwrap();
+    ortho.replace_material("Aluminium", Material { elastic: Some(carbon), density: Some(1.6e-9), ..Material::default() }).unwrap();
     assert!(unsupported(&ortho).contains("orthotropic material `Aluminium` in a shell section"));
 
     let mut loose = frequency();
@@ -658,7 +661,7 @@ fn what_abaqus_cannot_express_is_named() {
 fn orthotropic_solids_beams_and_rods_are_written() {
     let mut m = cantilever();
     let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] };
-    m.replace_material("Steel", Material { elastic: Some(carbon), density: None, fluid: None }).unwrap();
+    m.replace_material("Steel", Material { elastic: Some(carbon), density: None, ..Material::default() }).unwrap();
     let deck = abaqus::write(&m).unwrap();
     assert!(deck.contains("*ELASTIC, TYPE=ENGINEERING CONSTANTS\n140000.0, 10000.0, 10000.0, 0.3, 0.3, 0.4, 5000.0, 5000.0\n3800.0\n"), "{deck}");
 
@@ -699,3 +702,54 @@ fn a_model_that_is_not_ready_is_refused() {
     assert!(matches!(abaqus::write(&m), Err(Error::Model(_))));
     assert!(matches!(abaqus::write(&Model::default()), Err(Error::Model(_))));
 }
+
+/// The cantilever with its steel given each plasticity in turn.
+fn plastic_cantilever(p: modelcraft_model::Plastic) -> Model {
+    let mut m = cantilever();
+    m.replace_material("Steel", Material { plastic: Some(p), ..steel() }).unwrap();
+    m
+}
+
+/// Plasticity is written as `*PLASTIC` (and Johnson–Cook's rate and heat keywords); damage and
+/// erosion, failure models of Abaqus/Explicit, are reported.
+#[test]
+fn plasticity_is_written_and_failure_reported() {
+    use modelcraft_model::{Damage, Erosion, Hardening, HardeningPoint, Plastic};
+    let section = |deck: &str| deck.split("*MATERIAL, NAME=Steel").nth(1).and_then(|s| s.split("*SOLID SECTION").next()).unwrap().to_string();
+    let deck =
+        abaqus::write(&plastic_cantilever(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Isotropic }))
+            .unwrap();
+    assert!(section(&deck).contains("*PLASTIC\n250.0, 0.0\n10250.0, 10.0\n"), "{}", section(&deck));
+    let deck =
+        abaqus::write(&plastic_cantilever(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Kinematic }))
+            .unwrap();
+    assert!(section(&deck).contains("*PLASTIC, HARDENING=KINEMATIC\n250.0, 0.0\n1250.0, 1.0\n"), "{}", section(&deck));
+    let table = Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 200.0 }, HardeningPoint { strain: 0.1, stress: 300.0 }] };
+    let deck = abaqus::write(&plastic_cantilever(table)).unwrap();
+    assert!(section(&deck).contains("*PLASTIC\n200.0, 0.0\n300.0, 0.1\n"), "{}", section(&deck));
+    let jc = Plastic::JohnsonCook {
+        a: 90.0,
+        b: 292.0,
+        n: 0.31,
+        c: 0.025,
+        m: 1.09,
+        reference_strain_rate: 1.0,
+        room_temperature: 294.0,
+        melt_temperature: 1356.0,
+        specific_heat: Some(3.83e8),
+        taylor_quinney: 0.9,
+    };
+    let deck = abaqus::write(&plastic_cantilever(jc.clone())).unwrap();
+    let s = section(&deck);
+    assert!(
+        s.contains("*PLASTIC, HARDENING=JOHNSON COOK\n90.0, 292.0, 0.31, 1.09, 1356.0, 294.0\n*RATE DEPENDENT, TYPE=JOHNSON COOK\n0.025, 1.0\n*SPECIFIC HEAT\n383000000.0\n*INELASTIC HEAT FRACTION\n0.9\n"),
+        "{s}"
+    );
+    let mut failing = plastic_cantilever(jc);
+    let mut steel = failing.material("Steel").unwrap().clone();
+    steel.damage = Some(Damage::JohnsonCook { d1: 0.54, d2: 4.89, d3: -3.03, d4: 0.014, d5: 1.12 });
+    steel.erosion = Some(Erosion { plastic_strain: 1.0 });
+    failing.replace_material("Steel", steel).unwrap();
+    let Err(Error::Unsupported { what, .. }) = abaqus::write(&failing) else { panic!("damage is not written") };
+    assert!(what.contains("material `Steel`: damage and erosion"), "{what}");
+}
diff --git a/crates/decks/tests/lsdyna.rs b/crates/decks/tests/lsdyna.rs
index 62dc6d9..eb4eae0 100644
--- a/crates/decks/tests/lsdyna.rs
+++ b/crates/decks/tests/lsdyna.rs
@@ -13,7 +13,7 @@ use modelcraft_model::*;
 // Our own models.
 
 fn steel() -> Material {
-    Material { elastic: Some(Elastic::Isotropic { young: 210000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None }
+    Material { elastic: Some(Elastic::Isotropic { young: 210000.0, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() }
 }
 
 fn step(name: &str, procedure: Procedure, bcs: &[&str], loads: &[&str]) -> Step {
@@ -72,7 +72,7 @@ fn cantilever_modes() -> Model {
     let mut m = Model::default();
     cantilever_mesh(&mut m);
     let ortho = Elastic::EngineeringConstants { young: [140000.0, 10000.0, 10000.0], poisson: [0.3, 0.3, 0.45], shear: [5000.0, 5000.0, 3500.0] };
-    m.add_material("Ply", Material { elastic: Some(ortho), density: Some(1.6e-9), fluid: None }).unwrap();
+    m.add_material("Ply", Material { elastic: Some(ortho), density: Some(1.6e-9), ..Material::default() }).unwrap();
     let f = Formulation { incompatible_modes: true, ..Formulation::default() };
     m.add_section("beam", Section { set: "all".into(), material: "Ply".into(), kind: SectionKind::Solid {}, formulation: f }).unwrap();
     m.add_reference_point("RP", ReferencePoint { xyz: [5.0, 0.5, 0.5], label: 0 }).unwrap();
@@ -155,8 +155,11 @@ fn solids() -> Model {
     let h: Vec<u32> = h.iter().map(|p| node(&mut m, *p)).collect();
     m.set_element(4, ElementKind::Hex20, h).unwrap();
 
-    m.add_material("Aluminium", Material { elastic: Some(Elastic::Isotropic { young: 70e9, poisson: 0.33 }), density: Some(2700.0), fluid: None })
-        .unwrap();
+    m.add_material(
+        "Aluminium",
+        Material { elastic: Some(Elastic::Isotropic { young: 70e9, poisson: 0.33 }), density: Some(2700.0), ..Material::default() },
+    )
+    .unwrap();
     for (e, name) in [(1, "tet4"), (2, "tet10"), (3, "wedge"), (4, "hex20")] {
         m.add_set(name, Set { elements: vec![e], ..Set::default() }).unwrap();
         let f = Formulation { hybrid: e == 1, ..Formulation::default() };
@@ -658,7 +661,11 @@ fn hostile_numbers_stay_within_their_fields() {
     m.set_node(2, [f64::MIN_POSITIVE, -f64::MIN_POSITIVE, 5e-324]).unwrap();
     m.replace_material(
         "Steel",
-        Material { elastic: Some(Elastic::Isotropic { young: 1e300, poisson: 0.499_999_999_999_999_9 }), density: Some(1e-300), fluid: None },
+        Material {
+            elastic: Some(Elastic::Isotropic { young: 1e300, poisson: 0.499_999_999_999_999_9 }),
+            density: Some(1e-300),
+            ..Material::default()
+        },
     )
     .unwrap();
     m.replace_load(
@@ -687,3 +694,46 @@ fn explicit_steps_are_refused_by_every_format() {
         assert!(e.to_string().contains("step `Impact`: explicit dynamic"), "{}: {e}", f.name);
     }
 }
+
+/// Plasticity: `*MAT_PLASTIC_KINEMATIC` (bilinear) and `*MAT_PIECEWISE_LINEAR_PLASTICITY` with a
+/// curve numbered after the steps' ramps; the erosion strain is the failure strain; Johnson–Cook
+/// and damage are reported.
+#[test]
+fn plasticity_is_written_as_plastic_materials() {
+    let with = |p: Plastic, erosion: Option<f64>| {
+        let mut m = cantilever();
+        m.replace_material("Steel", Material { plastic: Some(p), erosion: erosion.map(|e| Erosion { plastic_strain: e }), ..steel() }).unwrap();
+        m
+    };
+    let w = lsdyna::write(&with(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Isotropic }, Some(0.5)))
+        .unwrap();
+    let mat = w.text.split("*MAT_PLASTIC_KINEMATIC_TITLE\nSteel\n").nth(1).unwrap();
+    let cards: Vec<Vec<&str>> = mat.lines().filter(|l| !l.starts_with('$')).take(2).map(|l| l.split_whitespace().collect()).collect();
+    // E_t = E H / (E + H) = 210000 · 1000 / 211000.
+    let etan: f64 = cards[0][5].parse().unwrap();
+    assert!((etan - 210000.0 * 1000.0 / 211000.0).abs() < 1e-9, "{cards:?}");
+    assert_eq!((cards[0][4], cards[0][6], cards[1][0]), ("250.0", "1.0", "0.5"), "{cards:?}");
+    assert!(!w.text.contains("*MAT_ELASTIC"));
+    let table = Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 200.0 }, HardeningPoint { strain: 0.1, stress: 300.0 }] };
+    let w = lsdyna::write(&with(table, None)).unwrap();
+    let mat = w.text.split("*MAT_PIECEWISE_LINEAR_PLASTICITY_TITLE\nSteel\n").nth(1).unwrap();
+    let cards: Vec<&str> = mat.lines().filter(|l| !l.starts_with('$')).take(2).collect();
+    assert_eq!(cards[1].split_whitespace().collect::<Vec<_>>(), ["2"], "the curve after the one step's ramp: {cards:?}");
+    let curve = w.text.split("*DEFINE_CURVE_TITLE\nSteel hardening\n").nth(1).unwrap();
+    let rows: Vec<Vec<&str>> = curve.lines().filter(|l| !l.starts_with('$')).take(3).map(|l| l.split_whitespace().collect()).collect();
+    assert_eq!(rows, [vec!["2"], vec!["0.0", "200.0"], vec!["0.1", "300.0"]]);
+    let jc = Plastic::JohnsonCook {
+        a: 90.0,
+        b: 292.0,
+        n: 0.31,
+        c: 0.0,
+        m: 1.0,
+        reference_strain_rate: 1.0,
+        room_temperature: 294.0,
+        melt_temperature: 1356.0,
+        specific_heat: None,
+        taylor_quinney: 0.9,
+    };
+    let Err(Error::Unsupported { what, .. }) = lsdyna::write(&with(jc, None)) else { panic!() };
+    assert!(what.contains("Johnson–Cook"), "{what}");
+}
diff --git a/crates/decks/tests/nastran.rs b/crates/decks/tests/nastran.rs
index b2856ff..ca79a66 100644
--- a/crates/decks/tests/nastran.rs
+++ b/crates/decks/tests/nastran.rs
@@ -131,7 +131,7 @@ impl Read {
 // Our own models.
 
 fn steel() -> Material {
-    Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None }
+    Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() }
 }
 
 /// A bar of `n` unit hex8 along X; node (i, y, z) = 1 + 4i + 2y + z.
@@ -190,7 +190,7 @@ fn plate() -> Model {
     m.add_set("top", Set { faces: (1..=4).map(|element| ElementFace { element, face: 1 }).collect(), ..Set::default() }).unwrap();
     m.add_material(
         "Aluminium",
-        Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9), fluid: None },
+        Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9), ..Material::default() },
     )
     .unwrap();
     let shell = SectionKind::Shell { thickness: 2.0, offset: 0.0, points: 5 };
@@ -521,7 +521,7 @@ fn beams_rods_orthotropy_couplings_and_prescribed_motion() {
     let ortho = Material {
         elastic: Some(Elastic::EngineeringConstants { young: [140e3, 10e3, 12e3], poisson: [0.3, 0.28, 0.4], shear: [5e3, 4.5e3, 3.8e3] }),
         density: Some(1.6e-9),
-        fluid: None,
+        ..Material::default()
     };
     m.add_material("Carbon", ortho.clone()).unwrap();
     m.replace_section(
@@ -676,7 +676,7 @@ fn what_nastran_cannot_express_is_named() {
     // Orthotropic beams, tilted orthotropic shells, one-point shells.
     let mut m = plate();
     let ortho = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] };
-    m.replace_material("Aluminium", Material { elastic: Some(ortho), density: Some(1.6e-9), fluid: None }).unwrap();
+    m.replace_material("Aluminium", Material { elastic: Some(ortho), density: Some(1.6e-9), ..Material::default() }).unwrap();
     m.set_node(9, [10.0, 10.0, 0.5]).unwrap();
     assert!(unsupported(&m).contains("planes of constant Z"));
     m.set_node(9, [10.0, 10.0, 0.0]).unwrap();
@@ -695,7 +695,7 @@ fn what_nastran_cannot_express_is_named() {
     let mut m = plate();
     m.set_element(9, ElementKind::Rod2, vec![1, 9]).unwrap();
     m.add_set("bars", Set { elements: vec![9], ..Set::default() }).unwrap();
-    m.add_material("Fibre", Material { elastic: Some(ortho), density: Some(1.6e-9), fluid: None }).unwrap();
+    m.add_material("Fibre", Material { elastic: Some(ortho), density: Some(1.6e-9), ..Material::default() }).unwrap();
     m.add_section(
         "bars",
         Section { set: "bars".into(), material: "Fibre".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() },
@@ -735,3 +735,50 @@ fn conflicting_values_and_ids_at_the_top_of_the_range() {
     m.add_coupling("pad", Coupling { reference_point: "RP".into(), set: "end".into(), kind: CouplingKind::Kinematic, dofs: Dofs::all() }).unwrap();
     assert!(matches!(nastran::write(&m), Err(Error::Limit(_))));
 }
+
+/// Plasticity: `MATS1` beside `MAT1` (bilinear by slope, a table through `TABLES1`); Johnson–Cook,
+/// damage and erosion are reported.
+#[test]
+fn plasticity_is_written_as_mats1() {
+    use modelcraft_model::{Erosion, Hardening, HardeningPoint, Plastic};
+    let with = |p: Plastic| {
+        let mut m = cantilever();
+        m.replace_material("Steel", Material { plastic: Some(p), ..steel() }).unwrap();
+        m
+    };
+    let deck = nastran::write(&with(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Kinematic })).unwrap();
+    let r = read(&deck);
+    let (_, f) = r.cards.iter().find(|(n, _)| n == "MATS1").unwrap();
+    assert_eq!((f[0].as_str(), f[1].as_str(), f[2].as_str(), int(&f[4]), int(&f[5])), ("1", "", "PLASTIC", 1, 2), "{f:?}");
+    assert_eq!((real(&f[3]), real(&f[6])), (1000.0, 250.0), "{f:?}");
+    let table = Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 210.0 }, HardeningPoint { strain: 0.1, stress: 315.0 }] };
+    let deck = nastran::write(&with(table)).unwrap();
+    let r = read(&deck);
+    let (_, mats1) = r.cards.iter().find(|(n, _)| n == "MATS1").unwrap();
+    let (_, table) = r.cards.iter().find(|(n, _)| n == "TABLES1").unwrap();
+    assert_eq!(mats1[1], table[0], "MATS1 names its table");
+    // The curve after the header line: the origin, then total strain σ/E + ε_p against σ.
+    let data: Vec<&String> = table[8..].iter().filter(|f| !f.is_empty()).collect();
+    assert_eq!(data.last().map(|s| s.as_str()), Some("ENDT"));
+    let values: Vec<f64> = data[..data.len() - 1].iter().map(|s| real(s)).collect();
+    let want = [0.0, 0.0, 0.001, 210.0, 0.1015, 315.0];
+    assert!(values.len() == 6 && values.iter().zip(want).all(|(a, b)| (a - b).abs() < 1e-12), "{values:?}");
+    let jc = Plastic::JohnsonCook {
+        a: 90.0,
+        b: 292.0,
+        n: 0.31,
+        c: 0.0,
+        m: 1.0,
+        reference_strain_rate: 1.0,
+        room_temperature: 294.0,
+        melt_temperature: 1356.0,
+        specific_heat: None,
+        taylor_quinney: 0.9,
+    };
+    let mut m = with(jc);
+    let mut steel = m.material("Steel").unwrap().clone();
+    steel.erosion = Some(Erosion { plastic_strain: 1.0 });
+    m.replace_material("Steel", steel).unwrap();
+    let Err(Error::Unsupported { what, .. }) = nastran::write(&m) else { panic!() };
+    assert!(what.contains("Johnson–Cook plasticity") && what.contains("damage and erosion"), "{what}");
+}
diff --git a/crates/decks/tests/planar.rs b/crates/decks/tests/planar.rs
index 01bd09f..3da995c 100644
--- a/crates/decks/tests/planar.rs
+++ b/crates/decks/tests/planar.rs
@@ -32,8 +32,11 @@ fn plate(kind: SectionKind, element: ElementKind, formulation: Formulation) -> M
     m.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
     m.add_set("left", Set { nodes: vec![1, 4], ..Set::default() }).unwrap();
     m.add_set("right", Set { edges: vec![right], ..Set::default() }).unwrap();
-    m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None })
-        .unwrap();
+    m.add_material(
+        "Steel",
+        Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() },
+    )
+    .unwrap();
     m.add_section("plate", Section { set: "all".into(), material: "Steel".into(), kind, formulation }).unwrap();
     m.add_boundary_condition("hold", BoundaryCondition::Fixed { target: Target::Set("left".into()) }).unwrap();
     m.add_load("press", Load::Pressure { set: "right".into(), magnitude: 5.0 }).unwrap();
diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs
index b71ede9..caac4a0 100644
--- a/crates/engine/src/catalogue.rs
+++ b/crates/engine/src/catalogue.rs
@@ -203,9 +203,9 @@ pub static CATALOGUE: &[Capability] = &[
         "Elastic (isotropic, orthotropic) material with density",
         &["model.material.edit", "model.material.delete"],
     ),
-    soon("model.material.plasticity", "Plasticity", Properties, "Materials", "M8", "Plastic hardening data"),
+    ready("model.material.plasticity", "Plasticity", Properties, "Materials", "Von Mises plasticity: bilinear, hardening table or Johnson–Cook"),
     soon("model.material.hyperelastic", "Hyperelastic", Properties, "Materials", "M8", "Hyperelastic material models"),
-    soon("model.material.damage", "Damage", Properties, "Materials", "M8", "Damage initiation and evolution"),
+    ready("model.material.damage", "Damage", Properties, "Materials", "Johnson–Cook damage and element erosion"),
     soon("model.material.thermal", "Thermal", Properties, "Materials", "M8", "Conductivity, specific heat and expansion"),
     soon(
         "model.material.eos",
diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs
index 4c6d09c..c95217a 100644
--- a/crates/engine/src/cmd/define.rs
+++ b/crates/engine/src/cmd/define.rs
@@ -80,7 +80,7 @@ definition_commands!(
     Material,
     "material",
     "Material",
-    r#"{"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3} | {"type": "engineering_constants", "young": [E1, E2, E3], "poisson": [ν12, ν13, ν23], "shear": [G12, G13, G23]}, "density": 7.85e-9, "fluid": {"type": "newtonian", "viscosity": 1.0e-9} (dynamic viscosity; a fluid needs a density)} — in the model's unit system (model.units), never converted"#,
+    r#"{"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3} | {"type": "engineering_constants", "young": [E1, E2, E3], "poisson": [ν12, ν13, ν23], "shear": [G12, G13, G23]}, "density": 7.85e-9, "fluid": {"type": "newtonian", "viscosity": 1.0e-9} (dynamic viscosity; a fluid needs a density), "plastic": {…} (see model.material.plasticity), "damage": {…}, "erosion": {…} (see model.material.damage)} — in the model's unit system (model.units), never converted"#,
     add_material,
     replace_material,
     remove_material,
@@ -89,6 +89,51 @@ definition_commands!(
     MATERIAL_DELETE
 );
 
+/// A copy of material `name` with `change` applied, put back through the model's validation.
+fn change_material(s: &mut Session, p: &Value, change: impl FnOnce(&mut Material, &Value) -> Result<()>) -> Result<Value> {
+    let n = name(p)?;
+    let mut m = s.model.material(n).cloned().ok_or_else(|| EngineError::Params(format!("material `{n}` does not exist")))?;
+    change(&mut m, p)?;
+    s.model.replace_material(n, m)?;
+    Ok(json!({"name": n}))
+}
+
+pub const MATERIAL_PLASTICITY: CommandSpec = CommandSpec {
+    id: "model.material.plasticity",
+    label: "Plasticity",
+    params: r#"{"name": "Steel", "plastic": {"type": "bilinear", "yield_stress": 250, "hardening_modulus": 1000, "hardening": "isotropic|kinematic"} | {"type": "table", "points": [{"strain": 0, "stress": 250}, {"strain": 0.1, "stress": 300}]} | {"type": "johnson_cook", "a": 792, "b": 510, "n": 0.26, "c": 0.014, "m": 1.03, "reference_strain_rate": 1, "room_temperature": 294, "melt_temperature": 1793, "specific_heat": 4.77e8, "taylor_quinney": 0.9} | null} — von Mises plasticity of an isotropic material (null removes it); explicit dynamics uses it, the linear solver the elastic constants alone"#,
+    kind: CommandKind::Model,
+    run: |s: &mut Session, p: &Value| {
+        change_material(s, p, |m, p| {
+            m.plastic = req(p, "plastic")?;
+            Ok(())
+        })
+    },
+};
+
+pub const MATERIAL_DAMAGE: CommandSpec = CommandSpec {
+    id: "model.material.damage",
+    label: "Damage",
+    params: r#"{"name": "Steel", "damage": {"type": "johnson_cook", "d1": 0.05, "d2": 3.44, "d3": -2.12, "d4": 0.002, "d5": 0.61} | null, "erosion": {"plastic_strain": 1.0} | null} — damage to failure and element erosion of a plastic material (each key optional; null removes it); a failed element is eroded in explicit dynamics"#,
+    kind: CommandKind::Model,
+    run: |s: &mut Session, p: &Value| {
+        change_material(s, p, |m, p| {
+            // Absent leaves a field as it is; null removes it.
+            let (damage, erosion) = (opt(p, "damage")?, opt(p, "erosion")?);
+            if damage.is_none() && erosion.is_none() {
+                return Err(EngineError::Params("give `damage`, `erosion` or both".into()));
+            }
+            if let Some(d) = damage {
+                m.damage = d;
+            }
+            if let Some(e) = erosion {
+                m.erosion = e;
+            }
+            Ok(())
+        })
+    },
+};
+
 definition_commands!(
     "section",
     Section,
@@ -423,4 +468,35 @@ pub(crate) mod tests {
         let e = s.execute("model.set.resolve", &json!({"name": "root"})).unwrap_err().to_string();
         assert!(e.contains("face left of body Bar is not on the mesh (the mesh carries no geometry tags)"), "{e}");
     }
+
+    #[test]
+    fn plasticity_and_damage_edit_a_material() {
+        let mut s = Session::default();
+        defined_bar(&mut s);
+        let bilinear = json!({"type": "bilinear", "yield_stress": 250, "hardening_modulus": 1000});
+        s.execute("model.material.plasticity", &json!({"name": "Steel", "plastic": bilinear})).unwrap();
+        assert!(s.model.material("Steel").unwrap().plastic.is_some());
+        // Damage needs plasticity; D4 and D5 need Johnson–Cook: refused, the model unchanged.
+        let before = s.model.clone();
+        let jc_damage = json!({"type": "johnson_cook", "d1": 0.05, "d2": 3.44, "d3": -2.12, "d4": 0.002, "d5": 0.61});
+        assert!(s.execute("model.material.damage", &json!({"name": "Steel", "damage": jc_damage})).is_err());
+        assert_eq!(s.model, before);
+        s.execute(
+            "model.material.damage",
+            &json!({"name": "Steel", "damage": {"type": "johnson_cook", "d1": 0.3, "d2": 0, "d3": 0}, "erosion": {"plastic_strain": 0.5}}),
+        )
+        .unwrap();
+        let m = s.model.material("Steel").unwrap();
+        assert!(m.damage.is_some() && m.erosion.is_some());
+        s.execute("model.material.damage", &json!({"name": "Steel", "erosion": null})).unwrap();
+        assert!(s.model.material("Steel").unwrap().erosion.is_none() && s.model.material("Steel").unwrap().damage.is_some());
+        // Removing plasticity under damage is refused.
+        assert!(s.execute("model.material.plasticity", &json!({"name": "Steel", "plastic": null})).is_err());
+        for bad in
+            [json!({"name": "Steel"}), json!({"name": "Nope", "plastic": null}), json!({"name": "Steel", "plastic": {"type": "table", "points": []}})]
+        {
+            assert!(s.execute("model.material.plasticity", &bad).is_err(), "{bad}");
+        }
+        assert!(s.execute("model.material.damage", &json!({"name": "Steel"})).is_err());
+    }
 }
diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs
index 7cf562d..f23d4eb 100644
--- a/crates/engine/src/cmd/mod.rs
+++ b/crates/engine/src/cmd/mod.rs
@@ -43,6 +43,8 @@ pub static ALL: &[CommandSpec] = &[
     define::MATERIAL_CREATE,
     define::MATERIAL_EDIT,
     define::MATERIAL_DELETE,
+    define::MATERIAL_PLASTICITY,
+    define::MATERIAL_DAMAGE,
     define::SECTION_CREATE,
     define::SECTION_EDIT,
     define::SECTION_DELETE,
diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs
index 5e9990f..3a0aa79 100644
--- a/crates/engine/src/cmd/solve.rs
+++ b/crates/engine/src/cmd/solve.rs
@@ -54,7 +54,7 @@ pub const FLUID: CommandSpec = CommandSpec {
 pub const PROBE: CommandSpec = CommandSpec {
     id: "solve.probe",
     label: "Probe Results",
-    params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, reaction, stress and von Mises (static), the shape of mode `mode` (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default); a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#,
+    params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, reaction, stress and von Mises (static), the shape of mode `mode` (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default) at `nodes` and, with `elements` (explicit), each element's stress, von Mises stress, strain, equivalent plastic strain, damage, temperature and whether it has eroded; a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#,
     kind: CommandKind::Model,
     run: probe,
 };
@@ -175,6 +175,7 @@ fn explicit_summary(r: &ExplicitResult) -> Value {
         "energy_balance": r.energy.relative_error(r.initial_kinetic),
         "frames": r.frames.len(),
         "history_samples": r.history.time.len(),
+        "eroded": r.eroded,
         "warnings": r.warnings,
     })
 }
@@ -344,15 +345,15 @@ fn probe(s: &mut Session, p: &Value) -> Result<Value> {
             format!("step `{step}` has no result; solve it first")
         })
     })?;
-    if let StepResult::Fluid(r) = result.as_ref() {
-        return probe_fluid(r, p);
-    }
-    let nodes = labels(p, "nodes", MAX_PROBE)?;
     let result = match result.as_ref() {
-        StepResult::Explicit(r) => return probe_frame(r, &step, &nodes, p),
+        StepResult::Explicit(r) => {
+            let list = |key| if p.get(key).is_some() { labels(p, key, MAX_PROBE) } else { Ok(Vec::new()) };
+            return probe_frame(r, &step, &list("nodes")?, &list("elements")?, p);
+        }
         StepResult::Fluid(r) => return probe_fluid(r, p),
         StepResult::Implicit(r) => r,
     };
+    let nodes = labels(p, "nodes", MAX_PROBE)?;
     let locate = |label: u32| {
         result
             .locate(label)
@@ -399,11 +400,38 @@ fn probe(s: &mut Session, p: &Value) -> Result<Value> {
     Ok(json!({"step": step, "nodes": out}))
 }
 
-/// Displacement and velocity at nodes in one field frame of an explicit result.
-fn probe_frame(r: &ExplicitResult, step: &str, nodes: &[u32], p: &Value) -> Result<Value> {
+/// Displacement and velocity at nodes, and element results, in one field frame of an explicit
+/// result.
+fn probe_frame(r: &ExplicitResult, step: &str, nodes: &[u32], elements: &[u32], p: &Value) -> Result<Value> {
     let last = r.frames.len().saturating_sub(1);
     let k = opt::<usize>(p, "frame")?.unwrap_or(last);
     let f = r.frames.get(k).ok_or_else(|| EngineError::Params(format!("`frame` must be from 0 to {last}")))?;
+    if nodes.is_empty() && elements.is_empty() {
+        return Err(EngineError::Params("give `nodes`, `elements` or both".into()));
+    }
+    if !elements.is_empty() && f.stress.is_empty() {
+        return Err(EngineError::Params(format!("step `{step}` does not output stresses")));
+    }
+    let element_values: Vec<Value> = elements
+        .iter()
+        .map(|&l| {
+            let i =
+                r.elements.binary_search(&l).map_err(|_| EngineError::Params(format!("element {l} has no result (it is not a solved element)")))?;
+            Ok(json!({
+                "element": l,
+                "stress": f.stress.get(i),
+                "von_mises": f.stress.get(i).map(modelcraft_fem::von_mises),
+                "strain": f.strain.get(i),
+                "plastic_strain": f.plastic_strain.get(i),
+                "damage": f.damage.get(i),
+                "temperature": f.temperature.get(i),
+                "eroded": f.failed.get(i),
+            }))
+        })
+        .collect::<Result<_>>()?;
+    if nodes.is_empty() {
+        return Ok(json!({"step": step, "frame": k, "time": f.time, "elements": element_values}));
+    }
     if f.displacement.is_empty() {
         return Err(EngineError::Params(format!("step `{step}` does not output displacements")));
     }
@@ -415,7 +443,7 @@ fn probe_frame(r: &ExplicitResult, step: &str, nodes: &[u32], p: &Value) -> Resu
             Ok(json!({"node": l, "displacement": f.displacement.get(i), "velocity": f.velocity.get(i)}))
         })
         .collect::<Result<_>>()?;
-    Ok(json!({"step": step, "frame": k, "time": f.time, "nodes": out}))
+    Ok(json!({"step": step, "frame": k, "time": f.time, "nodes": out, "elements": element_values}))
 }
 
 #[cfg(test)]
@@ -850,6 +878,11 @@ pub(crate) mod tests {
         assert_eq!(probe["nodes"][1]["velocity"][0].as_f64(), Some(-0.1), "{probe}");
         let last = s.execute("solve.probe", &json!({"step": "Impact", "nodes": [9]})).unwrap();
         assert_eq!(last["frame"], 4, "{last}");
+        let el = s.execute("solve.probe", &json!({"step": "Impact", "elements": [1, 8]})).unwrap();
+        let e = &el["elements"][0];
+        assert!(e["stress"].as_array().unwrap().len() == 6 && e["von_mises"].as_f64().unwrap() >= 0.0 && e["eroded"] == false, "{el}");
+        assert!(e["plastic_strain"].is_null(), "an elastic model has no plastic strain: {el}");
+        assert_eq!(r["eroded"].as_array().map(Vec::len), Some(0), "{r}");
         let h = s.execute("solve.history", &json!({"step": "Impact", "every": 2})).unwrap();
         assert_eq!(h["time"].as_array().unwrap().len(), h["energies"].as_array().unwrap().len());
         assert!(h["nodes"][0]["velocity"].as_array().unwrap().len() > 5, "{h}");
@@ -875,6 +908,11 @@ pub(crate) mod tests {
             json!({"step": "Impact", "nodes": [1], "frame": 99}),
             json!({"step": "Impact", "nodes": [1], "frame": -1}),
             json!({"step": "Impact", "nodes": [0, 4294967296u64]}),
+            json!({"step": "Impact"}),
+            json!({"step": "Impact", "elements": [999]}),
+            json!({"step": "Impact", "elements": [0, -1, 1.5]}),
+            json!({"step": "Impact", "elements": vec![1; 20_000]}),
+            json!({"step": "Impact", "elements": "all", "nodes": [1]}),
         ];
         for id in ["solve.explicit", "solve.history", "solve.probe"] {
             for p in &hostile {
diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs
index d47392d..dcce9d3 100644
--- a/crates/engine/src/lib.rs
+++ b/crates/engine/src/lib.rs
@@ -349,6 +349,14 @@ mod tests {
             json!({"name": "A", "elastic": {"type": "engineering_constants", "young": [1e-300, 1e300, 1], "poisson": [1e300, 0, 0], "shear": [0, 0, 0]}}),
             json!({"name": "A", "elastic": {"type": "isotropic", "young": "stiff"}, "density": -0.0}),
             json!({"name": "A", "set": "A", "material": "A", "kind": {"type": "solid"}}),
+            json!({"name": "Steel", "plastic": {"type": "bilinear", "yield_stress": -1e308, "hardening_modulus": f64::MAX}}),
+            json!({"name": "Steel", "plastic": {"type": "table", "points": vec![json!({"strain": 0, "stress": 1}); 2_000]}}),
+            json!({"name": "Steel", "plastic": {"type": "table", "points": [{"strain": 0, "stress": 1}, {"strain": 0, "stress": 2}]}}),
+            json!({"name": "Steel", "plastic": {"type": "johnson_cook", "a": 1, "b": 1, "n": -1, "room_temperature": 1e308, "melt_temperature": -1e308}}),
+            json!({"name": "Steel", "plastic": {"type": "johnson_cook", "a": 1, "b": 1e308, "n": 1e308, "c": 1e308, "m": 1e-308, "room_temperature": 0, "melt_temperature": 1, "taylor_quinney": 2}}),
+            json!({"name": "Steel", "plastic": "hard"}),
+            json!({"name": "Steel", "damage": {"type": "johnson_cook", "d1": 1e308, "d2": -1e308, "d3": 1e308}, "erosion": {"plastic_strain": -1}}),
+            json!({"name": "Steel", "damage": null, "erosion": {"plastic_strain": 1e-308}}),
             json!({"name": "A", "set": "all", "material": "Steel", "kind": {"type": "beam", "profile": {"type": "pipe", "radius": 1, "thickness": 2}, "orientation": [0, 0, 0]}}),
             json!({"name": "A", "set": "all", "material": "Steel", "kind": {"type": "shell", "thickness": 1e308, "offset": -0.5, "points": 4294967295u64}}),
             json!({"name": "A", "xyz": [1e308, 0, 0]}),
diff --git a/crates/explicit/src/elements.rs b/crates/explicit/src/elements.rs
index f718f40..211eb2e 100644
--- a/crates/explicit/src/elements.rs
+++ b/crates/explicit/src/elements.rs
@@ -1,24 +1,41 @@
-//! Element internal forces and critical time increments.
+//! Element internal forces, the material update and critical time increments.
+//!
+//! **Kinematics.** With nonlinear geometry the elements are updated Lagrangian (T. Belytschko,
+//! W. K. Liu, B. Moran, K. Elkhodary, *Nonlinear Finite Elements for Continua and Structures*,
+//! 2nd ed., Wiley 2014, ch. 4 and 6): the strain and spin increments of an increment come from
+//! the displacement increment Δu and the gradients of the **midpoint** configuration
+//! x_{n−½} = X + u_n − ½Δu, ΔG = (1/V) Σ_I Δu_I ⊗ b_I(x_{n−½}), which makes the update
+//! incrementally objective (T. J. R. Hughes, J. Winget, IJNME 15(12):1862–1867, 1980); the nodal
+//! forces of the Cauchy stress use the gradients of the current configuration,
+//! f_iI = σ_ij b_jI(x_n). Without it, both are the gradients of the undeformed mesh and the
+//! update is small-strain elasticity (and plasticity).
+//!
+//! **Stress update** in rate form: the material state (stress, back stress, the accumulated
+//! strain and the stiffness-form hourglass forces) is rotated with the increment's rotation —
+//! the Cayley transform of the spin increment (Jaumann rate, Hughes & Winget 1980) or the
+//! increment of the polar rotation (Green–Naghdi rate, D. P. Flanagan, L. M. Taylor, CMAME
+//! 62(3):305–320, 1987), see [`modelcraft_fem::rotation`] — then the constitutive increment is
+//! added: hypoelastic, or J2 radial return with Johnson–Cook damage
+//! ([`modelcraft_fem::plasticity`]). Orthotropic elasticity is updated in the polar frame (its
+//! axes turn with the material) whichever rate the step names.
 //!
 //! **One-point hexahedron** (D. P. Flanagan, T. Belytschko, IJNME 17(5):679–706, 1981): the
-//! stress of the mean strain gives the forces f_iI = σ_ij b_jI ([`modelcraft_fem::uniform`]).
-//! The hourglass modes are resisted through the shape vectors γ_α: with the hourglass rates
-//! q̇_iα = Σ_I v_iI γ_αI (or displacements q_iα = Σ_I u_iI γ_αI for the stiffness form), the
-//! generalized forces Q_iα = C q̇_iα (viscous) or k q_iα (stiffness) give the nodal forces
-//! f_iI = Σ_α Q_iα γ_αI. Because the γ_α are orthogonal to linear fields, these forces do no work
-//! in a constant strain state. The coefficients are scaled from the dimensionless hourglass
-//! coefficient κ with the element's size h = V / √(2 λ_max(Σ b bᵀ)) (the edge, for a cube),
-//! density ρ and wave speed c: C = κ ρ c h² / 16 damps a cube's hourglass modes at a fraction κ
-//! of critical at the frequency 2c/h, and k = κ ρ c² h / 16 puts their frequency at √κ · 2c/h.
-//! (These normalizations are ours; the paper leaves the scale of the coefficient to the user.)
+//! hourglass modes are resisted through the shape vectors γ_α of the current configuration:
+//! with the hourglass rates q̇_iα = Σ_I v_iI γ_αI, the generalized forces are Q_iα = C q̇_iα
+//! (viscous) or accumulate as Q_iα += k Σ_I Δu_iI γ_αI (stiffness, rotated with the element),
+//! and give the nodal forces f_iI = Σ_α Q_iα γ_αI. Because the γ_α are orthogonal to linear
+//! fields, these forces do no work in a constant strain state. The coefficients are scaled from
+//! the dimensionless hourglass coefficient κ with the element's current size
+//! h = V / √(2 λ_max(Σ b bᵀ)) (the edge, for a cube), density ρ and wave speed c:
+//! C = κ ρ c h² / 16 damps a cube's hourglass modes at a fraction κ of critical at the frequency
+//! 2c/h, and k = κ ρ c² h / 16 puts their frequency at √κ · 2c/h. (These normalizations are ours;
+//! the paper leaves the scale of the coefficient to the user.)
 //!
-//! **Average nodal pressure tetrahedron** (J. Bonet, A. J. Burton, *A simple average nodal
-//! pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit
-//! applications*, Commun. Numer. Methods Eng. 14(5):437–449, 1998): the volumetric strain of each
-//! node is the volume-weighted mean of its tetrahedra's, and each tetrahedron takes the mean of
-//! its nodes' as its own: σ = 2μ dev ε + K θ̄ I. The forces are those of the strain energy
-//! Σ_e V_e μ ‖dev ε_e‖² + ½ K Σ_a V_a θ_a², so the element conserves energy, and it does not lock
-//! when the material is nearly incompressible.
+//! **Average nodal pressure tetrahedron** (J. Bonet, A. J. Burton, Commun. Numer. Methods Eng.
+//! 14(5):437–449, 1998), in rate form: the volumetric strain increment of each node is the
+//! volume-weighted mean of its tetrahedra's, and each tetrahedron takes the mean of its nodes' as
+//! its own, so its pressure increment is K times that mean while the deviator is its own. It does
+//! not lock when the material is nearly incompressible or flows plastically.
 //!
 //! **Bulk viscosity** (J. von Neumann, R. D. Richtmyer, J. Appl. Phys. 21(3):232–237, 1950; M. L.
 //! Wilkins, J. Comput. Phys. 36(3):281–303, 1980): in compression (volumetric strain rate
@@ -27,60 +44,118 @@
 //! Δt = L / (Q + √(Q² + c²)), Q = C_l c + C_q L |ε̇ᵥ| (Wilkins 1980).
 //!
 //! **Critical increment.** Each element's stable length L = 2c/ω_max comes from the bound on its
-//! highest frequency ([`modelcraft_fem::uniform::UniformStrain::frequency_bound`], plus the
-//! hourglass stiffness when it is used), so without bulk viscosity its critical increment is
-//! 2/ω_max (T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, 2014, §6.3; the element eigenvalue
-//! inequality makes the smallest of them stable for the mesh). For the nodal-pressure
-//! tetrahedron the averaging only lowers the volumetric energy (the mean of squares bounds the
-//! square of the mean), so the plain tetrahedron's bound holds.
+//! highest frequency in its current shape ([`modelcraft_fem::uniform`], with the element's lumped
+//! masses, plus the hourglass stiffness when it is used), so without bulk viscosity its critical
+//! increment is 2/ω_max (Belytschko et al. 2014, §6.3). The elastic moduli bound the tangent
+//! moduli of plastic flow, so the bound holds for the inelastic materials too.
+//!
+//! **Erosion** (G. R. Johnson, R. A. Stryk, Int. J. Impact Eng. 5:411–421, 1987; T. Belytschko,
+//! J. I. Lin, ibid. 5:111–127, 1987): an element whose material point fails (damage 1) or
+//! 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 modelcraft_fem::material::stress;
-use modelcraft_fem::uniform::UniformStrain;
-use modelcraft_model::{BulkViscosity, HourglassKind};
+use modelcraft_fem::plasticity::{Increment, PointState};
+use modelcraft_fem::rotation::{self, IDENTITY, Mat3, Polar};
+use modelcraft_fem::uniform::{Hex8, Tet4, UniformStrain};
+use modelcraft_model::{BulkViscosity, HourglassKind, StressRate};
 
 use crate::system::{Material, NODE_DOFS, System};
+use crate::{Error, Result};
 
-/// Scratch vectors and the element state of the latest evaluation.
+/// Scratch vectors of one evaluation.
 #[derive(Debug, Clone)]
 pub struct Workspace {
     /// Internal forces per degree of freedom, by origin.
     pub stress_force: Vec<f64>,
     pub hourglass_force: Vec<f64>,
     pub viscous_force: Vec<f64>,
-    /// Per element (in [`System::elements`] order): the material stress and the strain.
-    pub stress: Vec<[f64; 6]>,
-    pub strain: Vec<[f64; 6]>,
-    /// Per node: the mean volumetric strain of its tetrahedra (a volume-weighted sum until it is
-    /// divided by `tet_volume`).
+    /// Per node: the volume-weighted volumetric strain increments of its tetrahedra, and their
+    /// volume.
     theta: Vec<f64>,
     tet_volume: Vec<f64>,
+    /// Per tetrahedron (in [`System::tets`] order): its increment between the two passes.
+    tet_increment: Vec<Option<TetIncrement>>,
 }
 
 impl Workspace {
     pub fn new(sys: &System) -> Workspace {
         let n = sys.labels.len();
-        let mut tet_volume = vec![0.0; n];
-        for t in &sys.tets {
-            for node in t.nodes {
-                if let Some(v) = tet_volume.get_mut(node) {
-                    *v += t.kinematics.volume;
-                }
-            }
-        }
         let dofs = sys.dofs();
-        let elements = sys.elements.len();
         Workspace {
             stress_force: vec![0.0; dofs],
             hourglass_force: vec![0.0; dofs],
             viscous_force: vec![0.0; dofs],
-            stress: vec![[0.0; 6]; elements],
-            strain: vec![[0.0; 6]; elements],
             theta: vec![0.0; n],
-            tet_volume,
+            tet_volume: vec![0.0; n],
+            tet_increment: vec![None; sys.tets.len()],
         }
     }
 }
 
+/// Why an element was eroded.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize)]
+#[serde(rename_all = "snake_case")]
+pub enum Failure {
+    /// Its damage reached 1.
+    Damage,
+    /// Its equivalent plastic strain reached the material's erosion strain.
+    PlasticStrain,
+}
+
+/// The state of one element (its single material point).
+#[derive(Debug, Clone)]
+pub struct ElementState {
+    pub point: PointState,
+    /// The accumulated (rotated) strain increments: engineering shears.
+    pub strain: [f64; 6],
+    /// The polar rotation and stretch (Green–Naghdi rate and orthotropic materials).
+    pub polar: Polar,
+    /// Stiffness-form hourglass generalized forces Q_α (one-point hexahedra).
+    pub hourglass: [[f64; 3]; 4],
+    /// Work of the stress on the element (its internal energy, hourglass apart).
+    pub energy: f64,
+    /// Plastic work on the element.
+    pub plastic_work: f64,
+    pub failed: Option<Failure>,
+}
+
+/// The element states and the erosions of a run.
+#[derive(Debug, Clone)]
+pub struct State {
+    /// Per element, in [`System::elements`] order.
+    pub elements: Vec<ElementState>,
+    /// Elements eroded during the latest evaluation: (index, failure).
+    pub eroded: Vec<(usize, Failure)>,
+}
+
+impl State {
+    pub fn new(sys: &System) -> State {
+        let mut material = vec![usize::MAX; sys.elements.len()];
+        for (index, m) in sys.hexes.iter().map(|h| (h.index, h.material)).chain(sys.tets.iter().map(|t| (t.index, t.material))) {
+            if let Some(slot) = material.get_mut(index) {
+                *slot = m;
+            }
+        }
+        let elements = material
+            .into_iter()
+            .map(|m| {
+                let temperature = sys.materials.get(m).and_then(|m| m.j2.as_ref()).map_or(0.0, |j| j.initial_temperature());
+                ElementState {
+                    point: PointState::new(temperature),
+                    strain: [0.0; 6],
+                    polar: Polar::default(),
+                    hourglass: [[0.0; 3]; 4],
+                    energy: 0.0,
+                    plastic_work: 0.0,
+                    failed: None,
+                }
+            })
+            .collect();
+        State { elements, eroded: Vec::new() }
+    }
+}
+
 pub(crate) fn gather<const N: usize>(x: &[f64], nodes: &[usize; N]) -> [[f64; 3]; N] {
     nodes.map(|n| [0, 1, 2].map(|k| x.get(NODE_DOFS * n + k).copied().unwrap_or(0.0)))
 }
@@ -112,102 +187,328 @@ pub fn critical_increment(rate: f64, length: f64, wave_speed: f64, b: &BulkVisco
     length / (q + (q * q + wave_speed * wave_speed).sqrt())
 }
 
-/// What the hexahedron and the tetrahedron share once their stress is known: the stress and bulk
-/// viscosity forces, the element state and the critical increment.
-struct Finish<'a, const N: usize> {
-    kinematics: &'a UniformStrain<N>,
-    nodes: &'a [usize; N],
-    index: usize,
-    length: f64,
+/// The largest absolute row sum of the hourglass vectors' Gram matrix Σ_I γ_αI γ_βI / m_I (an
+/// upper bound of its largest eigenvalue, exact for parallelepipeds where the γ_α are
+/// orthogonal).
+pub fn hourglass_gram_bound(gamma: &[[f64; 8]; 4], mass: &[f64; 8]) -> f64 {
+    let mut worst: f64 = 0.0;
+    for a in gamma {
+        let row: f64 = gamma.iter().map(|b| a.iter().zip(b).zip(mass).map(|((x, y), m)| x * y / m).sum::<f64>().abs()).sum();
+        worst = worst.max(row);
+    }
+    worst
 }
 
-impl<const N: usize> Finish<'_, N> {
-    fn apply(&self, ws: &mut Workspace, mat: &Material, v: &[f64], sigma: [f64; 6], strain: [f64; 6], bulk: &BulkViscosity) -> f64 {
-        let k = self.kinematics;
-        let rate = k.divergence(&gather(v, self.nodes));
-        let q = bulk_pressure(rate, mat.density, self.length, mat.wave_speed, bulk);
-        scatter(&mut ws.stress_force, self.nodes, &k.forces(&sigma));
-        if q != 0.0 {
-            scatter(&mut ws.viscous_force, self.nodes, &k.pressure_forces(q));
-        }
-        if let Some(s) = ws.stress.get_mut(self.index) {
-            *s = sigma;
-        }
-        if let Some(s) = ws.strain.get_mut(self.index) {
-            *s = strain;
-        }
-        critical_increment(rate, self.length, mat.wave_speed, bulk)
+/// What one evaluation of the internal forces is given.
+pub struct Motion<'a> {
+    /// Displacements u^n.
+    pub u: &'a [f64],
+    /// The displacement increment that led to `u`, Δt^{n−½} v^{n−½}.
+    pub du: &'a [f64],
+    /// The velocity v^{n−½}.
+    pub v: &'a [f64],
+    /// The increment's length Δt^{n−½} (0 at the start).
+    pub dt: f64,
+    pub time: f64,
+    pub bulk: &'a BulkViscosity,
+}
+
+/// The element's nodal values: rest positions, displacements, displacement increments and
+/// velocities, gathered once per evaluation.
+struct Nodal<const N: usize> {
+    x: [[f64; 3]; N],
+    u: [[f64; 3]; N],
+    du: [[f64; 3]; N],
+    v: [[f64; 3]; N],
+}
+
+impl<const N: usize> Nodal<N> {
+    fn gather(sys: &System, nodes: &[usize; N], m: &Motion) -> Nodal<N> {
+        Nodal { x: nodes.map(|n| sys.xyz.get(n).copied().unwrap_or_default()), u: gather(m.u, nodes), du: gather(m.du, nodes), v: gather(m.v, nodes) }
+    }
+
+    /// The positions X + u − s·Δu.
+    fn positions(&self, s: f64) -> [[f64; 3]; N] {
+        std::array::from_fn(|i| [0, 1, 2].map(|k| self.x[i][k] + self.u[i][k] - s * self.du[i][k]))
     }
 }
 
-/// Evaluate the internal forces at displacements `u` and velocities `v` into `ws` and return the
-/// smallest critical increment of the elements.
-pub fn internal_forces(sys: &System, u: &[f64], v: &[f64], bulk: &BulkViscosity, ws: &mut Workspace) -> f64 {
+fn inverted(label: u32, time: f64) -> Error {
+    Error::Unstable(format!("element {label} turned inside out at time {time:.6e} (its volume is no longer positive)"))
+}
+
+/// The rotation of an increment for an element's state, advancing its polar decomposition when
+/// it is tracked.
+fn rotation_of(sys: &System, mat: &Material, es: &mut ElementState, g: &Mat3, label: u32, time: f64) -> Result<Mat3> {
+    if !sys.nonlinear {
+        return Ok(IDENTITY);
+    }
+    if sys.stress_rate == StressRate::GreenNaghdi || mat.j2.is_none() {
+        es.polar.advance(g).ok_or_else(|| inverted(label, time))
+    } else {
+        let (_, w) = rotation::split(g);
+        rotation::cayley(&w).ok_or_else(|| inverted(label, time))
+    }
+}
+
+/// Rotate an element's state with `q`, then add the constitutive increment for the strain
+/// increment `d` (engineering shears) with volumetric part `vol` (the tetrahedron's averaged
+/// one; the trace of `d` for the hexahedron), track its energy at the midpoint volume `volume`,
+/// and check for erosion.
+#[allow(clippy::too_many_arguments)]
+fn update(mat: &Material, es: &mut ElementState, q: &Mat3, d: &[f64; 6], vol: f64, dt: f64, density: f64, volume: f64) {
+    let p = &mut es.point;
+    if *q != IDENTITY {
+        // Only the parts of the state that are in use (a back stress with kinematic hardening,
+        // hourglass forces with stiffness-form control).
+        p.stress = rotation::rotate(q, &p.stress);
+        if p.back != [0.0; 6] {
+            p.back = rotation::rotate(q, &p.back);
+        }
+        es.strain = rotation::rotate_strain(q, &es.strain);
+        if es.hourglass != [[0.0; 3]; 4] {
+            for h in &mut es.hourglass {
+                *h = rotation::apply(q, h);
+            }
+        }
+    }
+    let before = p.stress;
+    let inc = match &mat.j2 {
+        Some(j2) => j2.update(p, d, vol, dt, density),
+        None => {
+            // Orthotropic: in the polar frame, σ̂ += D d̂.
+            let r = es.polar.rotation;
+            if r == IDENTITY {
+                for (s, x) in p.stress.iter_mut().zip(stress(&mat.elasticity, d)) {
+                    *s += x;
+                }
+            } else {
+                let local = rotation::rotate_strain(&rotation::transpose(&r), d);
+                for (s, x) in p.stress.iter_mut().zip(rotation::rotate(&r, &stress(&mat.elasticity, &local))) {
+                    *s += x;
+                }
+            }
+            Increment::default()
+        }
+    };
+    // Stress work over the increment, ½(σ_n + σ_{n+1}) : Δε, with the volumetric strain the
+    // stress took.
+    let shift = (d[0] + d[1] + d[2] - vol) / 3.0;
+    let mut work = 0.0;
+    for k in 0..6 {
+        let dk = if k < 3 { d[k] - shift } else { d[k] };
+        work += 0.5 * (before[k] + p.stress[k]) * dk;
+        es.strain[k] += dk;
+    }
+    es.energy += work * volume;
+    es.plastic_work += inc.plastic_work * volume;
+    if es.failed.is_none() {
+        if p.damage >= 1.0 {
+            es.failed = Some(Failure::Damage);
+        } else if mat.erosion.is_some_and(|e| p.plastic_strain >= e) {
+            es.failed = Some(Failure::PlasticStrain);
+        }
+    }
+}
+
+/// The density and dilatational wave speed of an element of rest volume `rest`, now of volume
+/// `volume`.
+fn current(mat: &Material, rest: f64, volume: f64) -> (f64, f64) {
+    let rho = mat.density * rest / volume;
+    (rho, (mat.stiffest / rho).sqrt())
+}
+
+/// The bound ω² ≤ k λ_max(Σ b bᵀ) / (V m_min) on an element's highest frequency
+/// ([`modelcraft_fem::uniform`]), with the lumped masses of the undeformed element and
+/// `gradient_norm` = λ_max(Σ b bᵀ).
+pub fn frequency_bound<const N: usize>(k: &UniformStrain<N>, gradient_norm: f64, modulus: f64, mass: &[f64; N]) -> f64 {
+    let lightest = mass.iter().copied().fold(f64::INFINITY, f64::min);
+    modulus * gradient_norm / (k.volume * lightest)
+}
+
+/// A hexahedron's size h, hourglass coefficient (damping or stiffness) and frequency bound ω²
+/// (with the hourglass stiffness), at density `density`.
+#[derive(Debug, Clone, Copy)]
+pub struct HexScales {
+    pub size: f64,
+    pub coefficient: f64,
+    pub omega2: f64,
+}
+
+/// The [`HexScales`] of a hexahedron in configuration `k` with shape vectors `gamma`.
+pub fn hex_scales(
+    k: &Hex8,
+    gamma: &[[f64; 8]; 4],
+    mass: &[f64; 8],
+    mat: &Material,
+    density: f64,
+    hourglass: modelcraft_model::Hourglass,
+) -> HexScales {
+    let norm = k.gradient_norm();
+    // h: the edge of the cube with the element's volume and gradient norm.
+    let size = k.volume / (2.0 * norm).sqrt();
+    let wave2 = mat.stiffest / density;
+    let mut omega2 = frequency_bound(k, norm, mat.modulus, mass);
+    let kappa = hourglass.coefficient;
+    let coefficient = match hourglass.kind {
+        HourglassKind::Viscous => kappa * density * wave2.sqrt() * size * size / 16.0,
+        HourglassKind::Stiffness => {
+            let stiffness = kappa * density * wave2 * size / 16.0;
+            omega2 += stiffness * hourglass_gram_bound(gamma, mass);
+            stiffness
+        }
+    };
+    HexScales { size, coefficient, omega2 }
+}
+
+/// Evaluate the internal forces at the motion `m` into `ws`, updating the element states, and
+/// return the smallest critical increment of the elements still in the analysis.
+pub fn internal_forces(sys: &System, m: &Motion, state: &mut State, ws: &mut Workspace) -> Result<f64> {
     for f in [&mut ws.stress_force, &mut ws.hourglass_force, &mut ws.viscous_force] {
         f.fill(0.0);
     }
+    state.eroded.clear();
     let mut critical = f64::INFINITY;
+    let hourglass = modelcraft_model::Hourglass { kind: sys.hourglass, coefficient: sys.hourglass_coefficient };
     for h in &sys.hexes {
-        let Some(mat) = sys.materials.get(h.material) else { continue };
-        let ue = gather(u, &h.nodes);
-        let strain = h.kinematics.strain(&ue);
-        let sigma = stress(&mat.elasticity, &strain);
-        // Hourglass forces on the displacement (stiffness) or velocity (viscous) modes.
-        let modes = match sys.hourglass {
-            HourglassKind::Viscous => gather(v, &h.nodes),
-            HourglassKind::Stiffness => ue,
+        let (Some(mat), Some(es)) = (sys.materials.get(h.material), state.elements.get_mut(h.index)) else { continue };
+        if es.failed.is_some() {
+            continue;
+        }
+        let e = Nodal::gather(sys, &h.nodes, m);
+        let x_now = e.positions(0.0);
+        let (mid, now) = if sys.nonlinear {
+            let at = |x: &[[f64; 3]; 8]| Hex8::hex8_deformed(x).map_err(|_| inverted(h.label, m.time));
+            (at(&e.positions(0.5))?, at(&x_now)?)
+        } else {
+            (h.kinematics, h.kinematics)
+        };
+        let g = mid.gradient(&e.du);
+        let q = rotation_of(sys, mat, es, &g, h.label, m.time)?;
+        let (d, _) = rotation::split(&g);
+        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);
+        if let Some(failure) = es.failed {
+            state.eroded.push((h.index, failure));
+            continue;
+        }
+        // Hourglass forces with the current shape vectors and size (fixed without nonlinear
+        // geometry).
+        let (gamma, scales) = if sys.nonlinear {
+            let gamma = now.hourglass_vectors(&x_now);
+            let scales = hex_scales(&now, &gamma, &h.mass, mat, rho, hourglass);
+            (gamma, scales)
+        } else {
+            (h.gamma, h.rest)
         };
         let mut f = [[0.0; 3]; 8];
-        for g in &h.gamma {
-            let mut generalized = [0.0; 3];
-            for (x, gk) in modes.iter().zip(g) {
-                for i in 0..3 {
-                    generalized[i] += x[i] * gk;
+        let modes = if sys.hourglass == HourglassKind::Viscous { &e.v } else { &e.du };
+        for (k, g) in gamma.iter().enumerate() {
+            let rate: [f64; 3] = [0, 1, 2].map(|i| modes.iter().zip(g).map(|(x, gk)| x[i] * gk).sum());
+            // Viscous: Q = C q̇; stiffness: Q accumulates k Δq.
+            let generalized = match (sys.hourglass, es.hourglass.get_mut(k)) {
+                (HourglassKind::Stiffness, Some(qa)) => {
+                    for i in 0..3 {
+                        qa[i] += scales.coefficient * rate[i];
+                    }
+                    *qa
                 }
-            }
+                _ => rate.map(|r| scales.coefficient * r),
+            };
             for (fk, gk) in f.iter_mut().zip(g) {
                 for i in 0..3 {
-                    fk[i] += h.hourglass * generalized[i] * gk;
+                    fk[i] += generalized[i] * gk;
                 }
             }
         }
         scatter(&mut ws.hourglass_force, &h.nodes, &f);
-        let finish = Finish { kinematics: &h.kinematics, nodes: &h.nodes, index: h.index, length: h.length };
-        critical = critical.min(finish.apply(ws, mat, v, sigma, strain, bulk));
+        let element = Element { now: &now, mid: &mid, nodes: &h.nodes, density: rho, wave_speed: c, length: 2.0 * c / scales.omega2.sqrt() };
+        critical = critical.min(element.finish(ws, &es.point.stress, &e.v, m.bulk));
     }
     if !sys.tets.is_empty() {
-        // Pass 1: each tetrahedron's strain, and each node's mean volumetric strain.
+        // Pass 1: each tetrahedron's increment, and each node's mean volumetric increment.
         ws.theta.fill(0.0);
-        for t in &sys.tets {
-            let strain = t.kinematics.strain(&gather(u, &t.nodes));
+        ws.tet_volume.fill(0.0);
+        for (k, t) in sys.tets.iter().enumerate() {
+            let Some(slot) = ws.tet_increment.get_mut(k) else { continue };
+            *slot = None;
+            if state.elements.get(t.index).is_none_or(|e| e.failed.is_some()) {
+                continue;
+            }
+            let e = Nodal::gather(sys, &t.nodes, m);
+            let (mid, now) = if sys.nonlinear {
+                let at = |x: &[[f64; 3]; 4]| Tet4::tet4(x).map_err(|_| inverted(t.label, m.time));
+                (at(&e.positions(0.5))?, at(&e.positions(0.0))?)
+            } else {
+                (t.kinematics, t.kinematics)
+            };
+            let g = mid.gradient(&e.du);
             for node in t.nodes {
-                if let Some(x) = ws.theta.get_mut(node) {
-                    *x += t.kinematics.volume * (strain[0] + strain[1] + strain[2]);
+                if let (Some(x), Some(w)) = (ws.theta.get_mut(node), ws.tet_volume.get_mut(node)) {
+                    *x += mid.volume * (g[0][0] + g[1][1] + g[2][2]);
+                    *w += mid.volume;
                 }
             }
-            if let Some(s) = ws.strain.get_mut(t.index) {
-                *s = strain;
+            *slot = Some(TetIncrement { g, mid, now, v: e.v });
+        }
+        for (x, w) in ws.theta.iter_mut().zip(&ws.tet_volume) {
+            if *w > 0.0 {
+                *x /= w;
             }
         }
-        for (x, v) in ws.theta.iter_mut().zip(&ws.tet_volume) {
-            if *v > 0.0 {
-                *x /= v;
-            }
-        }
-        // Pass 2: each tetrahedron's pressure from the mean of its nodes' volumetric strains.
-        for t in &sys.tets {
-            let Some(mat) = sys.materials.get(t.material) else { continue };
-            let strain = ws.strain.get(t.index).copied().unwrap_or_default();
+        // Pass 2: each tetrahedron's update with the mean of its nodes' volumetric increments.
+        for (k, t) in sys.tets.iter().enumerate() {
+            let (Some(mat), Some(es), Some(Some(inc))) =
+                (sys.materials.get(t.material), state.elements.get_mut(t.index), ws.tet_increment.get(k).copied())
+            else {
+                continue;
+            };
             let averaged = t.nodes.iter().map(|n| ws.theta.get(*n).copied().unwrap_or(0.0)).sum::<f64>() / 4.0;
-            let own = strain[0] + strain[1] + strain[2];
-            let mut sigma = stress(&mat.elasticity, &strain);
-            let shift = mat.bulk.unwrap_or(0.0) * (averaged - own);
-            for s in sigma.iter_mut().take(3) {
-                *s += shift;
+            let q = rotation_of(sys, mat, es, &inc.g, t.label, m.time)?;
+            let (d, _) = rotation::split(&inc.g);
+            let (rho, c) = current(mat, t.kinematics.volume, inc.now.volume);
+            update(mat, es, &q, &d, averaged, m.dt, rho, inc.mid.volume);
+            if let Some(failure) = es.failed {
+                state.eroded.push((t.index, failure));
+                continue;
             }
-            let finish = Finish { kinematics: &t.kinematics, nodes: &t.nodes, index: t.index, length: t.length };
-            critical = critical.min(finish.apply(ws, mat, v, sigma, strain, bulk));
+            let omega2 = 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, nodes: &t.nodes, density: rho, wave_speed: c, length: 2.0 * c / omega2.sqrt() };
+            critical = critical.min(element.finish(ws, &es.point.stress, &inc.v, m.bulk));
         }
     }
-    critical
+    Ok(critical)
+}
+
+/// A tetrahedron's increment between the two passes: the gradient of its displacement increment,
+/// its midpoint and current kinematics, and its nodal velocities.
+#[derive(Debug, Clone, Copy)]
+struct TetIncrement {
+    g: Mat3,
+    mid: Tet4,
+    now: Tet4,
+    v: [[f64; 3]; 4],
+}
+
+/// What the hexahedron and the tetrahedron share once their stress is known.
+struct Element<'a, const N: usize> {
+    now: &'a UniformStrain<N>,
+    mid: &'a UniformStrain<N>,
+    nodes: &'a [usize; N],
+    density: f64,
+    wave_speed: f64,
+    length: f64,
+}
+
+impl<const N: usize> Element<'_, N> {
+    /// Add the stress and bulk viscosity forces and return the critical increment.
+    fn finish(&self, ws: &mut Workspace, sigma: &[f64; 6], v: &[[f64; 3]; N], bulk: &BulkViscosity) -> f64 {
+        let rate = self.mid.divergence(v);
+        let q = bulk_pressure(rate, self.density, self.length, self.wave_speed, bulk);
+        scatter(&mut ws.stress_force, self.nodes, &self.now.forces(sigma));
+        if q != 0.0 {
+            scatter(&mut ws.viscous_force, self.nodes, &self.now.pressure_forces(q));
+        }
+        critical_increment(rate, self.length, self.wave_speed, bulk)
+    }
 }
diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs
index 77d1b57..9c78e50 100644
--- a/crates/explicit/src/integrate.rs
+++ b/crates/explicit/src/integrate.rs
@@ -22,8 +22,8 @@
 
 use modelcraft_model::{ExplicitDynamic, Model, OutputVariable, Procedure, Step};
 
-use crate::elements::{Workspace, internal_forces};
-use crate::result::{ControllingElement, Energies, ExplicitResult, Frame, History};
+use crate::elements::{Motion, State, Workspace, internal_forces};
+use crate::result::{ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History};
 use crate::system::{NODE_DOFS, System};
 use crate::{Error, Limits, Result};
 
@@ -178,13 +178,17 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
     let total_mass: f64 = sys.mass.iter().sum();
     let mut ws = Workspace::new(&sys);
     let mut previous = Workspace::new(&sys);
+    let mut state = State::new(&sys);
+    let mut eroded: Vec<ErodedElement> = Vec::new();
     let mut u = vec![0.0; n];
     let mut du = vec![0.0; n];
     let v0 = sys.initial_velocity.clone();
     let initial_kinetic = 0.5 * v0.iter().zip(&mass).map(|(v, m)| m * v * v).sum::<f64>();
+    let inelastic = sys.materials.iter().any(|m| m.j2.as_ref().is_some_and(|j| j.flow.is_some()));
 
     // Start: forces at rest with the initial velocities, then the half increment.
-    let crit = internal_forces(&sys, &u, &v0, &bulk, &mut ws);
+    let start = Motion { u: &u, du: &du, v: &v0, dt: 0.0, time: 0.0, bulk: &bulk };
+    let crit = internal_forces(&sys, &start, &mut state, &mut ws)?;
     let mut dt_half = (scale * crit).min(end);
     let initial_increment = dt_half;
     let mut v_half = vec![0.0; n];
@@ -222,17 +226,28 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
         history.displacement.push(pick(u));
         history.velocity.push(pick(v));
     };
-    let frame = |t: f64, u: &[f64], v: &[f64], ws: &Workspace| Frame {
-        time: t,
-        displacement: if output.displacement { per_node(u) } else { Vec::new() },
-        velocity: if output.displacement { per_node(v) } else { Vec::new() },
-        stress: if output.stress { ws.stress.clone() } else { Vec::new() },
-        strain: if output.strain { ws.strain.clone() } else { Vec::new() },
+    let frame = |t: f64, u: &[f64], v: &[f64], state: &State| {
+        let each = |on: bool, f: &dyn Fn(&crate::elements::ElementState) -> f64| -> Vec<f64> {
+            if on { state.elements.iter().map(f).collect() } else { Vec::new() }
+        };
+        let inelastic = output.stress && inelastic;
+        Frame {
+            time: t,
+            displacement: if output.displacement { per_node(u) } else { Vec::new() },
+            velocity: if output.displacement { per_node(v) } else { Vec::new() },
+            stress: if output.stress { state.elements.iter().map(|e| e.point.stress).collect() } else { Vec::new() },
+            strain: if output.strain { state.elements.iter().map(|e| e.strain).collect() } else { Vec::new() },
+            plastic_strain: each(inelastic, &|e| e.point.plastic_strain),
+            damage: each(inelastic, &|e| e.point.damage),
+            temperature: each(inelastic, &|e| e.point.temperature),
+            failed: if output.stress { state.elements.iter().map(|e| e.failed.is_some()).collect() } else { Vec::new() },
+        }
     };
+    let plastic = |state: &State| state.elements.iter().map(|e| e.plastic_work).sum::<f64>();
     let mut histories = Schedule { interval: settings.history_interval(), next: settings.history_interval() };
     let mut frames_due = Schedule { interval: settings.output_interval(), next: settings.output_interval() };
     sample(&mut history, 0.0, dt_half, energy, &u, &v0, &reaction);
-    out_frames.push(frame(0.0, &u, &v0, &ws));
+    out_frames.push(frame(0.0, &u, &v0, &state));
 
     let mut t = 0.0;
     let mut increments = 0u64;
@@ -253,7 +268,14 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
         t += dt_half;
         increments += 1;
         std::mem::swap(&mut ws, &mut previous);
-        let crit = internal_forces(&sys, &u, &v_half, &bulk, &mut ws);
+        let motion = Motion { u: &u, du: &du, v: &v_half, dt: dt_half, time: t, bulk: &bulk };
+        let crit = internal_forces(&sys, &motion, &mut state, &mut ws)?;
+        for &(index, failure) in &state.eroded {
+            let element = state.elements.get(index);
+            let internal = element.map_or(0.0, |e| e.energy);
+            energy.eroded += internal;
+            eroded.push(ErodedElement { label: sys.elements.get(index).copied().unwrap_or(0), time: t, failure, internal_energy: internal });
+        }
         let remaining = end - t;
         let last = remaining <= 1e-12 * end;
         let dt_next = if last { 0.0 } else { (scale * crit).min(remaining) };
@@ -269,8 +291,8 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
         let [mut internal, mut hourglass, mut viscous, mut external, mut constraint, mut kinetic] = [0.0; 6];
         let forces = ws.stress_force.iter().zip(&ws.hourglass_force).zip(&ws.viscous_force);
         let before = previous.stress_force.iter().zip(&previous.hourglass_force).zip(&previous.viscous_force);
-        let state = v_now.iter_mut().zip(&v_half).zip(&v_new).zip(&mass).zip(&du).zip(&sys.external).zip(&reaction);
-        for ((((((((now, vh), vn), m), dx), e), r), ((s, h), q)), ((s0, h0), q0)) in state.zip(forces).zip(before) {
+        let per_dof = v_now.iter_mut().zip(&v_half).zip(&v_new).zip(&mass).zip(&du).zip(&sys.external).zip(&reaction);
+        for ((((((((now, vh), vn), m), dx), e), r), ((s, h), q)), ((s0, h0), q0)) in per_dof.zip(forces).zip(before) {
             *now = vh + weight * (vn - vh);
             internal += dx * (s0 + s);
             hourglass += dx * (h0 + h);
@@ -287,14 +309,16 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
         energy.kinetic = 0.5 * kinetic;
         peak_internal = peak_internal.max(energy.internal.abs());
         energy.balance = energy.kinetic + energy.internal + energy.hourglass + energy.viscous - energy.external - energy.constraint - initial_kinetic;
-        if !(energy.balance.is_finite() && crit.is_finite()) {
+        // Every element eroded leaves an infinite critical increment: the nodes fly freely.
+        if !(energy.balance.is_finite() && !crit.is_nan()) {
             return Err(Error::Unstable(format!("the solution is not finite at time {t:.6e} (increment {increments})")));
         }
         if histories.due(t) || last {
+            energy.plastic = plastic(&state);
             sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction);
         }
         if frames_due.due(t) || last {
-            out_frames.push(frame(t, &u, &v_now, &ws));
+            out_frames.push(frame(t, &u, &v_now, &state));
         }
         std::mem::swap(&mut v_half, &mut v_new);
         dt_half = dt_next;
@@ -307,9 +331,17 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
             sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction);
         }
         if out_frames.last().map(|f| f.time) != Some(t) {
-            out_frames.push(frame(t, &u, &v_now, &ws));
+            out_frames.push(frame(t, &u, &v_now, &state));
         }
     }
+    energy.plastic = plastic(&state);
+    if !eroded.is_empty() {
+        warnings.push(format!(
+            "{} elements eroded (their material failed); their internal energy at erosion was {:.4e}",
+            eroded.len(),
+            energy.eroded
+        ));
+    }
 
     let error = energy.relative_error(initial_kinetic);
     if error > BALANCE_TOLERANCE {
@@ -338,6 +370,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res
         total_mass,
         energy,
         initial_kinetic,
+        eroded,
         warnings,
     })
 }
diff --git a/crates/explicit/src/lib.rs b/crates/explicit/src/lib.rs
index 4dbe072..c6dbfb0 100644
--- a/crates/explicit/src/lib.rs
+++ b/crates/explicit/src/lib.rs
@@ -17,10 +17,12 @@
 //! - **Integration** ([`integrate`]): the time loop, kinematic constraints, energy accounting
 //!   with a balance check, history samples and field frames ([`result`]).
 //!
-//! This first version is geometrically linear (small displacements and strains; the elements'
-//! gradients are those of the undeformed mesh) with linear elastic materials, which is what the
-//! analytic verification problems need. Large deformation kinematics with objective stress rates
-//! and inelastic materials come next (docs/solvers/literature.md, Part F, X2).
+//! With nonlinear geometry (the default) the elements are updated Lagrangian and the material
+//! state is integrated in rate form with an objective stress rate (Jaumann by Hughes–Winget,
+//! Green–Naghdi by Flanagan–Taylor); materials are hypoelastic or J2 plastic (bilinear,
+//! tabulated, Johnson–Cook) with Johnson–Cook damage and element erosion ([`elements`];
+//! docs/solvers/literature.md, Parts F and J). Without it the gradients stay those of the
+//! undeformed mesh (small displacements).
 //!
 //! Everything is sequential, with elements in label order and no parallel reductions, so the
 //! same model gives the same numbers bit for bit; sizes, increments and output are capped by
@@ -35,8 +37,9 @@ pub mod integrate;
 pub mod result;
 pub mod system;
 
+pub use elements::Failure;
 pub use integrate::solve_explicit;
-pub use result::{ControllingElement, Energies, ExplicitResult, Frame, History};
+pub use result::{ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History};
 
 #[derive(Debug, Clone, PartialEq, thiserror::Error)]
 pub enum Error {
diff --git a/crates/explicit/src/result.rs b/crates/explicit/src/result.rs
index f4519f7..70bc4b1 100644
--- a/crates/explicit/src/result.rs
+++ b/crates/explicit/src/result.rs
@@ -3,6 +3,8 @@
 
 use serde::Serialize;
 
+use crate::elements::Failure;
+
 /// The energy account at one time (Belytschko, Liu, Moran, Elkhodary 2014, §6.2.3). The works of
 /// the internal and external forces are integrated with the trapezoidal rule over each
 /// increment's displacement; the constraints' with their impulses (see [`crate::integrate`]).
@@ -20,6 +22,10 @@ pub struct Energies {
     pub external: f64,
     /// Work of the boundary conditions and rigid walls on the model.
     pub constraint: f64,
+    /// Plastic work: the part of `internal` dissipated by plastic flow.
+    pub plastic: f64,
+    /// Internal energy of the eroded elements when they eroded (part of `internal`).
+    pub eroded: f64,
     /// kinetic + internal + hourglass + viscous − external − constraint − the initial kinetic
     /// energy: zero when energy is conserved.
     pub balance: f64,
@@ -45,7 +51,26 @@ pub struct Frame {
     /// Per element in [`ExplicitResult::elements`] order (one-point elements: one value each):
     /// xx, yy, zz, xy, xz, yz, engineering shear strains.
     pub stress: Vec<[f64; 6]>,
+    /// The accumulated strain increments (small strain: the strain; large deformation: the
+    /// integral of the rate of deformation, rotated with the material).
     pub strain: Vec<[f64; 6]>,
+    /// Per element, with stress output in a model with plasticity: equivalent plastic strain,
+    /// damage and temperature.
+    pub plastic_strain: Vec<f64>,
+    pub damage: Vec<f64>,
+    pub temperature: Vec<f64>,
+    /// Per element, with stress output: whether it has eroded.
+    pub failed: Vec<bool>,
+}
+
+/// An element removed from the analysis.
+#[derive(Debug, Clone, Copy, PartialEq, Serialize)]
+pub struct ErodedElement {
+    pub label: u32,
+    pub time: f64,
+    pub failure: Failure,
+    /// Its internal energy (stress work) when it eroded.
+    pub internal_energy: f64,
 }
 
 /// Quantities sampled through the run.
@@ -95,5 +120,7 @@ pub struct ExplicitResult {
     /// The energy account at the end.
     pub energy: Energies,
     pub initial_kinetic: f64,
+    /// Elements eroded during the run, in the order they eroded.
+    pub eroded: Vec<ErodedElement>,
     pub warnings: Vec<String>,
 }
diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs
index 51c3f70..e326e4f 100644
--- a/crates/explicit/src/system.rs
+++ b/crates/explicit/src/system.rs
@@ -3,11 +3,14 @@
 //! vector of an explicit dynamic step.
 
 use modelcraft_fem::material::Elasticity;
+use modelcraft_fem::plasticity::J2;
 use modelcraft_fem::uniform::{Hex8, Tet4, modulus_bound};
 use modelcraft_fem::{DofKind, DofLayout, DofSet, face};
-use modelcraft_model::{Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, Target};
+use modelcraft_model::{
+    Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, StressRate, Target,
+};
 
-use crate::elements::add_at;
+use crate::elements::{HexScales, add_at, frequency_bound, hex_scales};
 use crate::{Error, Result};
 
 /// Degrees of freedom per node: the three translations, node-major (see [`System`]).
@@ -27,8 +30,13 @@ pub struct Material {
     pub wave_speed: f64,
     /// The bound of [`modelcraft_fem::uniform::modulus_bound`].
     pub modulus: f64,
-    /// Bulk modulus, for the nodal-pressure tetrahedron (isotropic materials).
-    pub bulk: Option<f64>,
+    /// The largest normal stiffness M (λ + 2μ when isotropic): c = √(M/ρ).
+    pub stiffest: f64,
+    /// The hypoelastic(–plastic) update of an isotropic material; `None` for orthotropic
+    /// elasticity, updated with [`Material::elasticity`].
+    pub j2: Option<J2>,
+    /// The equivalent plastic strain at which elements erode.
+    pub erosion: Option<f64>,
 }
 
 /// A one-point hexahedron.
@@ -41,9 +49,10 @@ pub struct Hex {
     pub material: usize,
     pub kinematics: Hex8,
     pub gamma: [[f64; 8]; 4],
-    /// The hourglass coefficient: a damping (viscous) or a stiffness (stiffness form), see
-    /// [`crate::elements`].
-    pub hourglass: f64,
+    /// Row-sum lumped masses of the element's nodes.
+    pub mass: [f64; 8],
+    /// Size, hourglass coefficient and frequency bound of the undeformed element.
+    pub rest: HexScales,
     /// The element's stable length: its critical increment without bulk viscosity is
     /// `length / wave_speed`.
     pub length: f64,
@@ -57,6 +66,9 @@ pub struct Tet {
     pub nodes: [usize; 4],
     pub material: usize,
     pub kinematics: Tet4,
+    pub mass: [f64; 4],
+    /// The frequency bound ω² of the undeformed element.
+    pub omega2: f64,
     pub length: f64,
 }
 
@@ -97,6 +109,10 @@ pub struct System {
     /// Nodes whose histories are recorded.
     pub history_nodes: Vec<usize>,
     pub hourglass: HourglassKind,
+    pub hourglass_coefficient: f64,
+    /// Large deformation (updated Lagrangian) or the undeformed mesh's gradients.
+    pub nonlinear: bool,
+    pub stress_rate: StressRate,
     pub warnings: Vec<String>,
 }
 
@@ -114,8 +130,11 @@ fn material(model: &Model, name: &str) -> Result<Material> {
         .ok_or_else(|| invalid(format!("material `{name}` has no density, and explicit dynamics needs one")))?;
     let d = modelcraft_fem::elasticity(elastic).map_err(|e| invalid(format!("material `{name}`: {e}")))?;
     let stiffest = d[0][0].max(d[1][1]).max(d[2][2]);
-    let bulk = match *elastic {
-        Elastic::Isotropic { young, poisson } => Some(young / (3.0 * (1.0 - 2.0 * poisson))),
+    // The model keeps plasticity, damage and erosion to isotropic materials.
+    let j2 = match *elastic {
+        Elastic::Isotropic { .. } => {
+            Some(J2::from_model(elastic, m.plastic.as_ref(), m.damage.as_ref()).map_err(|e| invalid(format!("material `{name}`: {e}")))?)
+        }
         Elastic::EngineeringConstants { .. } => None,
     };
     Ok(Material {
@@ -124,22 +143,12 @@ fn material(model: &Model, name: &str) -> Result<Material> {
         density,
         wave_speed: (stiffest / density).sqrt(),
         modulus: modulus_bound(elastic, &d),
-        bulk,
+        stiffest,
+        j2,
+        erosion: m.erosion.map(|e| e.plastic_strain),
     })
 }
 
-/// The largest absolute row sum of the hourglass vectors' Gram matrix Σ_I γ_αI γ_βI / v_I (an
-/// upper bound of its largest eigenvalue, exact for parallelepipeds where the γ_α are
-/// orthogonal).
-fn hourglass_gram_bound(gamma: &[[f64; 8]; 4], nodal_volume: &[f64; 8]) -> f64 {
-    let mut worst: f64 = 0.0;
-    for a in gamma {
-        let row: f64 = gamma.iter().map(|b| a.iter().zip(b).zip(nodal_volume).map(|((x, y), v)| x * y / v).sum::<f64>().abs()).sum();
-        worst = worst.max(row);
-    }
-    worst
-}
-
 /// The section's elements of a kind the explicit solver has, or why not.
 fn check_formulation(kind: ElementKind, section_name: &str, f: &modelcraft_model::Formulation) -> Result<()> {
     let refuse = |why: &str| Err(Error::NotSupported(format!("section `{section_name}`: {why}")));
@@ -183,7 +192,7 @@ impl System {
             for &label in &set.elements {
                 let el = model.element(label).ok_or_else(|| invalid(format!("element {label} is missing")))?;
                 check_formulation(el.kind, name, &section.formulation)?;
-                if el.kind == ElementKind::Tet4 && materials.get(mi).is_some_and(|m| m.bulk.is_none()) {
+                if el.kind == ElementKind::Tet4 && materials.get(mi).is_some_and(|m| m.j2.is_none()) {
                     return Err(Error::NotSupported(format!(
                         "nodal-pressure tet4 elements of the orthotropic material `{}` (section `{name}`)",
                         section.material
@@ -222,6 +231,9 @@ impl System {
             external: vec![0.0; dofs],
             history_nodes: Vec::new(),
             hourglass: settings.hourglass.kind,
+            hourglass_coefficient: settings.hourglass.coefficient,
+            nonlinear: settings.nonlinear_geometry,
+            stress_rate: settings.stress_rate,
             warnings: resolved.warnings.iter().map(ToString::to_string).collect(),
             labels,
             xyz,
@@ -237,28 +249,20 @@ impl System {
                     let x = nodes.map(|i| sys.xyz.get(i).copied().unwrap_or_default());
                     let k = Hex8::hex8(&x).map_err(at)?;
                     let gamma = k.hourglass_vectors(&x);
-                    // h: the edge of the cube with the element's volume and gradient norm.
-                    let h = k.volume / (2.0 * k.gradient_norm()).sqrt();
-                    let kappa = settings.hourglass.coefficient;
-                    let mut omega2 = k.frequency_bound(mat.modulus, rho);
-                    let hourglass = match settings.hourglass.kind {
-                        HourglassKind::Viscous => kappa * rho * c * h * h / 16.0,
-                        HourglassKind::Stiffness => {
-                            let stiffness = kappa * rho * c * c * h / 16.0;
-                            omega2 += stiffness * hourglass_gram_bound(&gamma, &k.nodal_volume) / rho;
-                            stiffness
-                        }
-                    };
-                    sys.lump(&nodes, &k.nodal_volume, rho);
-                    sys.hexes.push(Hex { label, index, nodes, material: mi, kinematics: k, gamma, hourglass, length: 2.0 * c / omega2.sqrt() });
+                    let mass = k.nodal_volume.map(|v| rho * v);
+                    let rest = hex_scales(&k, &gamma, &mass, mat, rho, settings.hourglass);
+                    sys.lump(&nodes, &mass);
+                    let length = 2.0 * c / rest.omega2.sqrt();
+                    sys.hexes.push(Hex { label, index, nodes, material: mi, kinematics: k, gamma, mass, rest, length });
                 }
                 ElementKind::Tet4 => {
                     let nodes: [usize; 4] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 4 nodes")))?;
                     let x = nodes.map(|i| sys.xyz.get(i).copied().unwrap_or_default());
                     let k = Tet4::tet4(&x).map_err(at)?;
-                    let omega2 = k.frequency_bound(mat.modulus, rho);
-                    sys.lump(&nodes, &k.nodal_volume, rho);
-                    sys.tets.push(Tet { label, index, nodes, material: mi, kinematics: k, length: 2.0 * c / omega2.sqrt() });
+                    let mass = k.nodal_volume.map(|v| rho * v);
+                    let omega2 = frequency_bound(&k, k.gradient_norm(), mat.modulus, &mass);
+                    sys.lump(&nodes, &mass);
+                    sys.tets.push(Tet { label, index, nodes, material: mi, kinematics: k, mass, omega2, length: 2.0 * c / omega2.sqrt() });
                 }
                 other => return Err(Error::Internal(format!("{} element kind", other.name()))),
             }
@@ -270,10 +274,10 @@ impl System {
     }
 
     /// Add an element's row-sum lumped masses.
-    fn lump(&mut self, nodes: &[usize], nodal_volume: &[f64], density: f64) {
-        for (node, v) in nodes.iter().zip(nodal_volume) {
+    fn lump(&mut self, nodes: &[usize], mass: &[f64]) {
+        for (node, v) in nodes.iter().zip(mass) {
             if let Some(m) = self.mass.get_mut(*node) {
-                *m += density * v;
+                *m += v;
             }
         }
     }
diff --git a/crates/explicit/tests/large_deformation.rs b/crates/explicit/tests/large_deformation.rs
new file mode 100644
index 0000000..caf8afc
--- /dev/null
+++ b/crates/explicit/tests/large_deformation.rs
@@ -0,0 +1,482 @@
+//! Large deformation and plasticity in explicit dynamics against closed-form solutions
+//! (docs/solvers/literature.md, A4, A6, A10):
+//! - simple shear of a hypoelastic material (J. K. Dienes, *On the analysis of rotation and
+//!   stress rate in deforming bodies*, Acta Mech. 32:217–232, 1979): the Jaumann rate gives the
+//!   oscillating σ12 = G sin γ, σ11 = −σ22 = G(1 − cos γ), the Green–Naghdi rate the monotonic
+//!   solution in β = arctan(γ/2);
+//! - a stressed element turned rigidly by 90° in many increments: the stress turns with it and
+//!   keeps its magnitude (objectivity), where the small-strain update does not;
+//! - uniaxial stress under large stretch with J2 plasticity, loaded and unloaded: Cauchy stress
+//!   against the logarithmic strain follows the bilinear flow curve and unloads elastically;
+//! - the Taylor bar (G. I. Taylor, Proc. R. Soc. A 194:289–299, 1948): the yield stress Taylor's
+//!   one-dimensional estimate recovers from the computed final shape.
+//!
+//! Run with `--nocapture` to see the errors.
+
+use modelcraft_explicit::elements::{Motion, State, Workspace, internal_forces};
+use modelcraft_explicit::system::System;
+use modelcraft_explicit::{ExplicitResult, Failure, Limits, solve_explicit};
+use modelcraft_model::*;
+
+/// A block [0, size] of `n` cells: one hex8 per cell, or six tet4 (Kuhn's split).
+fn block(n: [usize; 3], size: [f64; 3], kind: ElementKind) -> Model {
+    let mut m = Model::default();
+    let label = |i: usize, j: usize, k: usize| (1 + i + (n[0] + 1) * (j + (n[1] + 1) * k)) as u32;
+    for k in 0..=n[2] {
+        for j in 0..=n[1] {
+            for i in 0..=n[0] {
+                let x = [i as f64 * size[0] / n[0] as f64, j as f64 * size[1] / n[1] as f64, k as f64 * size[2] / n[2] as f64];
+                m.set_node(label(i, j, k), x).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[2] {
+        for j in 0..n[1] {
+            for i in 0..n[0] {
+                let at = |c: [usize; 3]| label(i + c[0], j + c[1], k + c[2]);
+                if kind == ElementKind::Hex8 {
+                    m.set_element(next, kind, 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 p = v.map(|c| [0, 1, 2].map(|d| c[d] as f64 * size[d] / n[d] as f64));
+                    let e = |a: usize| [0, 1, 2].map(|d| p[a][d] - p[0][d]);
+                    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, kind, 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();
+    m
+}
+
+fn node_set(m: &mut Model, name: &str, f: impl Fn([f64; 3]) -> bool) {
+    let nodes: Vec<u32> = m.nodes().iter().filter(|(_, n)| f(n.xyz)).map(|(l, _)| *l).collect();
+    m.add_set(name, Set { nodes, ..Set::default() }).unwrap();
+}
+
+fn define(m: &mut Model, kind: ElementKind, material: Material, step: ExplicitDynamic) {
+    m.add_material("Mat", material).unwrap();
+    let integration = if kind == ElementKind::Hex8 { Integration::Reduced } else { Integration::Full };
+    let formulation = Formulation { integration, ..Formulation::default() };
+    m.add_section("solid", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
+    m.add_step(Step::new("Dynamic", Procedure::ExplicitDynamic(step)), None).unwrap();
+}
+
+fn bc(m: &mut Model, name: &str, bc: BoundaryCondition) {
+    m.add_boundary_condition(name, bc).unwrap();
+    m.set_boundary_condition_active("Dynamic", name, true).unwrap();
+}
+
+fn elastic(young: f64, poisson: f64, density: f64) -> Material {
+    Material { elastic: Some(Elastic::Isotropic { young, poisson }), density: Some(density), ..Material::default() }
+}
+
+/// Dienes' closed forms: (σ11, σ12) / G at shear γ.
+fn dienes(rate: StressRate, gamma: f64) -> (f64, f64) {
+    match rate {
+        StressRate::Jaumann => (1.0 - gamma.cos(), gamma.sin()),
+        StressRate::GreenNaghdi => {
+            let beta = (gamma / 2.0).atan();
+            let lc = beta.cos().ln();
+            let s11 = 4.0 * ((2.0 * beta).cos() * lc + beta * (2.0 * beta).sin() - beta.sin().powi(2));
+            let s12 = 2.0 * (2.0 * beta).cos() * (2.0 * beta - 2.0 * (2.0 * beta).tan() * lc - beta.tan());
+            (s11, s12)
+        }
+    }
+}
+
+/// A unit cube sheared to γ = 6 by prescribing every node: x = X + γ(t) Y e_x.
+#[test]
+fn dienes_simple_shear_with_both_rates() {
+    let (young, poisson) = (1.0, 0.3);
+    let g = young / (2.0 * (1.0 + poisson));
+    let gamma_end = 6.0;
+    for kind in [ElementKind::Hex8, ElementKind::Tet4] {
+        for rate in [StressRate::Jaumann, StressRate::GreenNaghdi] {
+            let mut m = block([1, 1, 1], [1.0; 3], kind);
+            let end = 4000.0;
+            let mut step = ExplicitDynamic::new(end);
+            step.stress_rate = rate;
+            step.output_interval = Some(end / 60.0);
+            step.max_increments = 100_000;
+            define(&mut m, kind, elastic(young, poisson, 1.0), step);
+            node_set(&mut m, "bottom", |x| x[1] == 0.0);
+            node_set(&mut m, "top", |x| x[1] == 1.0);
+            bc(&mut m, "hold", BoundaryCondition::Fixed { target: Target::Set("bottom".into()) });
+            bc(
+                &mut m,
+                "shear",
+                BoundaryCondition::Displacement {
+                    target: Target::Set("top".into()),
+                    values: [Some(gamma_end), Some(0.0), Some(0.0), None, None, None],
+                },
+            );
+            let r = solve_explicit(&m, None, &Limits::default()).unwrap();
+            assert!(r.completed);
+            let mut worst: f64 = 0.0;
+            for f in &r.frames {
+                let gamma = gamma_end * f.time / end;
+                let (s11, s12) = dienes(rate, gamma);
+                for s in &f.stress {
+                    let e = [(s[0] / g - s11).abs(), (s[1] / g + s11).abs(), (s[3] / g - s12).abs(), s[2].abs() / g, s[4].abs() / g, s[5].abs() / g];
+                    worst = worst.max(e.into_iter().fold(0.0, f64::max));
+                }
+            }
+            let last = r.frames.last().unwrap();
+            eprintln!(
+                "Dienes {} {rate:?}: σ/G at γ = 6: ({:.5}, {:.5}) vs ({:.5}, {:.5}); largest error over the run {worst:.2e} G; {} increments",
+                kind.name(),
+                last.stress[0][0] / g,
+                last.stress[0][3] / g,
+                dienes(rate, gamma_end).0,
+                dienes(rate, gamma_end).1,
+                r.increments
+            );
+            assert!(worst < 2e-3, "{} {rate:?}: {worst}", kind.name());
+        }
+    }
+}
+
+/// Drive one element through a sequence of nodal displacements with the solver's own element
+/// update, returning its stress after each.
+struct Driver {
+    sys: System,
+    state: State,
+    ws: Workspace,
+    u: Vec<f64>,
+    time: f64,
+}
+
+impl Driver {
+    fn new(model: &Model) -> Driver {
+        let resolved = model.resolve().unwrap();
+        let step = model.step("Dynamic").unwrap();
+        let Procedure::ExplicitDynamic(settings) = &step.procedure else { panic!() };
+        let sys = System::new(model, &resolved, step, settings, 1000).unwrap();
+        let (state, ws) = (State::new(&sys), Workspace::new(&sys));
+        let u = vec![0.0; sys.dofs()];
+        Driver { sys, state, ws, u, time: 0.0 }
+    }
+
+    /// Move the nodes to `x(X)` over an increment of length `dt`.
+    fn go(&mut self, to: impl Fn([f64; 3]) -> [f64; 3], dt: f64) -> [f64; 6] {
+        self.try_go(to, dt).unwrap()
+    }
+
+    fn try_go(&mut self, to: impl Fn([f64; 3]) -> [f64; 3], dt: f64) -> modelcraft_explicit::Result<[f64; 6]> {
+        let mut du = vec![0.0; self.u.len()];
+        for (i, x0) in self.sys.xyz.iter().enumerate() {
+            let x = to(*x0);
+            for k in 0..3 {
+                du[3 * i + k] = x[k] - x0[k] - self.u[3 * i + k];
+                self.u[3 * i + k] = x[k] - x0[k];
+            }
+        }
+        self.time += dt;
+        let v: Vec<f64> = du.iter().map(|d| d / dt).collect();
+        let bulk = BulkViscosity { linear: 0.0, quadratic: 0.0 };
+        let m = Motion { u: &self.u, du: &du, v: &v, dt, time: self.time, bulk: &bulk };
+        internal_forces(&self.sys, &m, &mut self.state, &mut self.ws)?;
+        Ok(self.state.elements[0].point.stress)
+    }
+}
+
+fn rotation_z(angle: f64, x: [f64; 3]) -> [f64; 3] {
+    let (s, c) = angle.sin_cos();
+    [c * x[0] - s * x[1], s * x[0] + c * x[1], x[2]]
+}
+
+/// Stretch an element 1% along x, then turn it rigidly by 90° about z in 50 increments: the
+/// stress turns with it (σ_yy takes σ_xx's value), exactly for both rates, and the internal
+/// forces stay those of the turned stress. The small-strain update keeps σ_xx and adds a
+/// spurious stress from the rotation's "strain".
+#[test]
+fn rigid_rotation_turns_the_stress_with_the_body() {
+    for (rate, nonlinear) in [(StressRate::Jaumann, true), (StressRate::GreenNaghdi, true), (StressRate::Jaumann, false)] {
+        let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
+        let mut step = ExplicitDynamic::new(1.0);
+        step.stress_rate = rate;
+        step.nonlinear_geometry = nonlinear;
+        define(&mut m, ElementKind::Hex8, elastic(1000.0, 0.3, 1.0), step);
+        let mut d = Driver::new(&m);
+        let stretch = |x: [f64; 3]| [1.01 * x[0], x[1], x[2]];
+        let before = d.go(stretch, 1.0);
+        let n = 50;
+        let mut after = before;
+        for k in 1..=n {
+            let angle = std::f64::consts::FRAC_PI_2 * k as f64 / n as f64;
+            after = d.go(|x| rotation_z(angle, stretch(x)), 1.0);
+        }
+        // R σ Rᵀ for 90° about z: xx ↔ yy, xy → −xy, xz → yz, yz → −xz.
+        let want = [before[1], before[0], before[2], -before[3], -before[5], before[4]];
+        let error = after.iter().zip(want).map(|(a, b)| (a - b).abs()).fold(0.0, f64::max) / before[0].abs();
+        eprintln!("rigid rotation {rate:?} nonlinear {nonlinear}: σ before {before:.4?}, after {after:.4?}; error {error:.2e} of σ_xx");
+        if nonlinear {
+            assert!(error < 1e-12, "{rate:?}: {error}");
+            // Forces of the turned stress balance on the turned element: ∑ f = 0, ∑ x × f = 0.
+            let f = &d.ws.stress_force;
+            let total: f64 = (0..3).map(|k| f.iter().skip(k).step_by(3).sum::<f64>().abs()).sum();
+            assert!(total < 1e-10 * before[0].abs());
+        } else {
+            assert!(error > 0.5, "the small-strain update is not objective: {error}");
+        }
+    }
+}
+
+/// Uniaxial stress to a logarithmic strain of 0.3 and back by 0.006 (inside the elastic range) in a single hexahedron with
+/// bilinear isotropic hardening: the lateral stretch is iterated each increment so the lateral
+/// stresses vanish. The Cauchy stress follows σ = σ_y + H ε_p with ln λ = σ/E + ε_p (for a
+/// hypoelastic–plastic law under coaxial stretching the integral of D is the logarithmic
+/// strain), then unloads with slope E.
+#[test]
+fn j2_uniaxial_stress_strain_curve_and_unloading_at_large_strain() {
+    let (young, poisson, y0, h) = (200_000.0, 0.3, 250.0, 2000.0);
+    for kind in [ElementKind::Hex8, ElementKind::Tet4] {
+        let mut m = block([1, 1, 1], [1.0; 3], kind);
+        let plastic = Plastic::Bilinear { yield_stress: y0, hardening_modulus: h, hardening: Hardening::Isotropic };
+        define(&mut m, kind, Material { plastic: Some(plastic), ..elastic(young, poisson, 7.8e-9) }, ExplicitDynamic::new(1.0));
+        let mut d = Driver::new(&m);
+        let mut path: Vec<f64> = (1..=600).map(|k| 0.3 * k as f64 / 600.0).collect();
+        path.extend((1..=200).map(|k| 0.3 - 0.006 * k as f64 / 200.0));
+        let (mut lateral, mut worst, mut peak) = (0.0f64, 0.0f64, (0.0, 0.0, 0.0));
+        let mut last = [0.0; 6];
+        for (step, axial) in path.iter().enumerate() {
+            // Secant iteration on the lateral log strain for σ_yy = 0 (σ_zz follows by symmetry).
+            let saved = (d.state.clone(), d.u.clone(), d.time);
+            let trial = |lat: f64, d: &mut Driver| {
+                (d.state, d.u, d.time) = saved.clone();
+                let (a, l) = (axial.exp(), lat.exp());
+                d.go(|x| [a * x[0], l * x[1], l * x[2]], 1e-6)
+            };
+            let (mut a, mut b) = (lateral, lateral - 1e-4);
+            let (mut fa, mut fb) = (trial(a, &mut d)[1], trial(b, &mut d)[1]);
+            for _ in 0..60 {
+                if fb.abs() < 1e-9 || fb == fa {
+                    break;
+                }
+                let c = b - fb * (b - a) / (fb - fa);
+                (a, fa) = (b, fb);
+                b = c;
+                fb = trial(b, &mut d)[1];
+            }
+            lateral = b;
+            let s = trial(b, &mut d);
+            let p = d.state.elements.iter().map(|e| e.point.plastic_strain).fold(0.0, f64::max);
+            if step < 600 {
+                // Loading: σ = y0 + H ε_p and ε = σ/E + ε_p, so σ = y0 + E_t (ε − y0/E).
+                let et = young * h / (young + h);
+                let exact = if *axial * young <= y0 { young * axial } else { y0 + et * (axial - y0 / young) };
+                worst = worst.max((s[0] - exact).abs() / exact);
+                peak = (*axial, s[0], p);
+            } else {
+                // Unloading: elastic from the peak.
+                let exact = peak.1 - young * (peak.0 - axial);
+                worst = worst.max((s[0] - exact).abs() / peak.1);
+                assert!((p - peak.2).abs() < 1e-12, "no plastic flow while unloading");
+            }
+            last = s;
+        }
+        eprintln!(
+            "J2 uniaxial {}: peak σ {:.4} at ln λ = 0.3 (exact {:.4}); after unloading to 0.294 σ = {:.4}; largest relative error {worst:.2e}",
+            kind.name(),
+            peak.1,
+            y0 + young * h / (young + h) * (0.3 - y0 / young),
+            last[0]
+        );
+        assert!(worst < 1e-6, "{}: {worst}", kind.name());
+    }
+}
+
+/// The Taylor bar: a copper-like bar of square section hits a rigid wall at 190 m/s (units mm,
+/// t, s, MPa) and mushrooms. Taylor's one-dimensional analysis of a rigid–perfectly plastic bar
+/// (Taylor 1948) relates the final length L₁ and the length X of the undeformed rear part to the
+/// yield stress: σ_y = ρ U² (L − X) / (2 (L − L₁) ln(L / X)). Applied to the computed shape it
+/// should give back the yield stress the material has, within the accuracy of the estimate
+/// itself (which neglects elastic waves and radial inertia): 25%.
+#[test]
+fn taylor_bar_recovers_its_yield_stress() {
+    let (rho, young, poisson, yield_stress) = (8.93e-9, 117_000.0, 0.35, 400.0);
+    let (length, width, speed) = (25.4, 7.6, 190_000.0);
+    let end = 80e-6;
+    let mut step = ExplicitDynamic::new(end);
+    step.output_interval = Some(end / 4.0);
+    step.history_interval = Some(end / 400.0);
+    step.initial_velocities = vec![InitialVelocity { set: "nodes".into(), velocity: [0.0, 0.0, -speed] }];
+    step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0] }];
+    step.hourglass = Hourglass { kind: HourglassKind::Stiffness, coefficient: 0.1 };
+    let mut m = block([4, 4, 30], [width, width, length], ElementKind::Hex8);
+    node_set(&mut m, "nodes", |_| true);
+    let plastic = Plastic::Bilinear { yield_stress, hardening_modulus: 0.0, hardening: Hardening::Isotropic };
+    define(&mut m, ElementKind::Hex8, Material { plastic: Some(plastic), ..elastic(young, poisson, rho) }, step);
+    let r = solve_explicit(&m, None, &Limits::default()).unwrap();
+    let f = r.frames.last().unwrap();
+    // Final length: the top face's mean height; undeformed part: the axis above the last
+    // element with plastic strain, along the centre column of elements.
+    let z = |label: u32| m.node(label).unwrap().xyz[2] + f.displacement[r.nodes.binary_search(&label).unwrap()][2];
+    let top: Vec<u32> = m.nodes().iter().filter(|(_, n)| n.xyz[2] == length).map(|(l, _)| *l).collect();
+    let l1 = top.iter().map(|l| z(*l)).sum::<f64>() / top.len() as f64;
+    let column: Vec<usize> = (0..30).map(|k| 4 * 4 * k + 4 + 1).collect();
+    // The rear part counts as undeformed from the first element (from the impact end) whose
+    // plastic strain is below 1%, small against the mushroom's.
+    let deformed = column.iter().position(|&i| f.plastic_strain[i] < 0.01).unwrap_or(30);
+    let x = length * (30 - deformed) as f64 / 30.0;
+    let estimate = rho * speed * speed * (length - x) / (2.0 * (length - l1) * (length / x).ln());
+    let kinetic = r.energy.kinetic / r.initial_kinetic;
+    eprintln!(
+        "Taylor bar: L₁/L = {:.4}, X/L = {:.4}; Taylor's estimate σ_y = {estimate:.1} vs {yield_stress} ({:+.1}%); final KE/KE₀ {kinetic:.4}; plastic work {:.3e} of {:.3e}; balance {:.2e}; {} increments",
+        l1 / length,
+        x / length,
+        100.0 * (estimate / yield_stress - 1.0),
+        r.energy.plastic,
+        r.initial_kinetic,
+        r.energy.relative_error(r.initial_kinetic),
+        r.increments
+    );
+    assert!(r.completed && l1 < 0.95 * length && l1 > 0.5 * length);
+    assert!((estimate / yield_stress - 1.0).abs() < 0.25, "{estimate}");
+    assert!(r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy);
+    assert!(r.energy.plastic > 0.5 * r.initial_kinetic, "most of the impact energy goes into plastic work: {:?}", r.energy);
+}
+
+fn taylor_with(material: Material) -> (Model, ExplicitResult) {
+    let end = 40e-6;
+    let mut step = ExplicitDynamic::new(end);
+    step.output_interval = Some(end / 2.0);
+    step.initial_velocities = vec![InitialVelocity { set: "nodes".into(), velocity: [0.0, 0.0, -300_000.0] }];
+    step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0] }];
+    step.hourglass = Hourglass { kind: HourglassKind::Stiffness, coefficient: 0.1 };
+    let mut m = block([2, 2, 12], [7.6, 7.6, 25.4], ElementKind::Hex8);
+    node_set(&mut m, "nodes", |_| true);
+    define(&mut m, ElementKind::Hex8, material, step);
+    let r = solve_explicit(&m, None, &Limits::default()).unwrap();
+    (m, r)
+}
+
+/// Johnson–Cook strength, adiabatic heating, damage and erosion in a fast impact: elements at
+/// the impact face reach the erosion strain (compression keeps their damage low, as the
+/// fracture strain grows with confinement), are removed and reported; the run stays balanced
+/// and repeatable bit for bit.
+#[test]
+fn johnson_cook_impact_erodes_elements() {
+    // OFHC copper constants as published by Johnson & Cook (1983, 1985), units mm, t, s, MPa, K.
+    let plastic = Plastic::JohnsonCook {
+        a: 90.0,
+        b: 292.0,
+        n: 0.31,
+        c: 0.025,
+        m: 1.09,
+        reference_strain_rate: 1.0,
+        room_temperature: 294.0,
+        melt_temperature: 1356.0,
+        specific_heat: Some(3.83e8),
+        taylor_quinney: 0.9,
+    };
+    let damage = Damage::JohnsonCook { d1: 0.54, d2: 4.89, d3: -3.03, d4: 0.014, d5: 1.12 };
+    let material = Material {
+        plastic: Some(plastic),
+        damage: Some(damage),
+        erosion: Some(Erosion { plastic_strain: 0.8 }),
+        ..elastic(124_000.0, 0.34, 8.96e-9)
+    };
+    let (_, r) = taylor_with(material.clone());
+    let f = r.frames.last().unwrap();
+    let hottest = f.temperature.iter().copied().fold(0.0, f64::max);
+    eprintln!(
+        "Johnson–Cook impact: {} elements eroded ({:?}…), eroded energy {:.3e}; hottest element {hottest:.1} K; largest damage {:.3}; balance {:.2e}",
+        r.eroded.len(),
+        r.eroded.first(),
+        r.energy.eroded,
+        f.damage.iter().copied().fold(0.0, f64::max),
+        r.energy.relative_error(r.initial_kinetic)
+    );
+    assert!(!r.eroded.is_empty() && r.eroded.iter().all(|e| matches!(e.failure, Failure::Damage | Failure::PlasticStrain)));
+    assert_eq!(f.failed.iter().filter(|x| **x).count(), r.eroded.len());
+    assert!(hottest > 294.0 + 10.0);
+    assert!(r.energy.eroded > 0.0 && r.warnings.iter().any(|w| w.contains("eroded")));
+    assert!(r.energy.relative_error(r.initial_kinetic) < 0.02, "{:?}", r.energy);
+    let (_, again) = taylor_with(material);
+    assert!(r == again, "bit-identical runs");
+}
+
+/// Damage with a constant fracture strain D1 (D2 = D4 = D5 = 0): an element stretched in
+/// tension fails when its plastic strain reaches D1, and from then on carries no stress.
+#[test]
+fn damage_fails_an_element_at_its_fracture_strain() {
+    let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
+    let plastic = Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Isotropic };
+    let damage = Damage::JohnsonCook { d1: 0.2, d2: 0.0, d3: 0.0, d4: 0.0, d5: 0.0 };
+    define(
+        &mut m,
+        ElementKind::Hex8,
+        Material { plastic: Some(plastic), damage: Some(damage), ..elastic(200_000.0, 0.3, 7.8e-9) },
+        ExplicitDynamic::new(1.0),
+    );
+    let mut d = Driver::new(&m);
+    let mut failed_at = None;
+    for k in 1..=1000 {
+        let stretch = 1.0 + 0.5 * k as f64 / 1000.0;
+        let before = d.state.elements[0].point.plastic_strain;
+        d.go(|x| [stretch * x[0], x[1] / stretch.sqrt(), x[2] / stretch.sqrt()], 1e-6);
+        let e = &d.state.elements[0];
+        if e.failed.is_some() {
+            failed_at = Some((before, e.point.plastic_strain, e.failed));
+            break;
+        }
+    }
+    let (before, at, why) = failed_at.expect("the element fails");
+    eprintln!("damage: failed between plastic strain {before:.5} and {at:.5} (D1 = 0.2)");
+    assert_eq!(why, Some(Failure::Damage));
+    assert!(before < 0.2 && at >= 0.2 - 1e-12, "{before} {at}");
+    // Failed: no forces from it.
+    d.go(|x| [1.6 * x[0], x[1], x[2]], 1e-6);
+    assert!(d.ws.stress_force.iter().all(|f| *f == 0.0));
+}
+
+/// Inelastic models that cannot run are refused by name, never panics.
+#[test]
+fn inelastic_hostile_models_are_errors() {
+    // Plasticity with orthotropic elasticity is refused by the model itself.
+    let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
+    let ortho = Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [0.5; 3] };
+    let plastic = Plastic::Bilinear { yield_stress: 1.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic };
+    assert!(m.add_material("Bad", Material { elastic: Some(ortho), density: Some(1.0), plastic: Some(plastic), ..Material::default() }).is_err());
+    // An element turned inside out stops the run with an error naming it.
+    define(&mut m, ElementKind::Hex8, elastic(1.0, 0.3, 1.0), ExplicitDynamic::new(1.0));
+    let mut d = Driver::new(&m);
+    d.go(|x| [x[0], x[1], 0.5 * x[2]], 1.0);
+    let e = d.try_go(|x| [x[0], x[1], 0.5 - x[2]], 1.0).unwrap_err().to_string();
+    assert!(e.contains("element 1") && e.contains("inside out"), "{e}");
+    // A crushed element shrinks the increment until the cap stops the run, with a warning.
+    let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
+    let mut step = ExplicitDynamic::new(1.0);
+    step.max_increments = 20_000;
+    define(&mut m, ElementKind::Hex8, elastic(1.0, 0.3, 1.0), step);
+    node_set(&mut m, "top", |x| x[2] == 1.0);
+    node_set(&mut m, "bottom", |x| x[2] == 0.0);
+    bc(&mut m, "hold", BoundaryCondition::Fixed { target: Target::Set("bottom".into()) });
+    bc(
+        &mut m,
+        "crush",
+        BoundaryCondition::Displacement { target: Target::Set("top".into()), values: [Some(0.0), Some(0.0), Some(-2.0), None, None, None] },
+    );
+    match solve_explicit(&m, None, &Limits::default()) {
+        Ok(r) => assert!(!r.completed && r.warnings.iter().any(|w| w.contains("cap"))),
+        Err(e) => assert!(e.to_string().contains("element 1"), "{e}"),
+    }
+}
diff --git a/crates/explicit/tests/verification.rs b/crates/explicit/tests/verification.rs
index 99a2e3b..a90c394 100644
--- a/crates/explicit/tests/verification.rs
+++ b/crates/explicit/tests/verification.rs
@@ -77,8 +77,11 @@ impl Bar {
             }
         }
         m.add_set("end_face", Set { faces, ..Set::default() }).unwrap();
-        m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: Some(DENSITY), fluid: None })
-            .unwrap();
+        m.add_material(
+            "Mat",
+            Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: Some(DENSITY), ..Material::default() },
+        )
+        .unwrap();
         let integration = if kind == ElementKind::Hex8 { Integration::Reduced } else { Integration::Full };
         let formulation = Formulation { integration, ..Formulation::default() };
         m.add_section("bar", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
@@ -334,12 +337,12 @@ fn caps_and_unsupported_models_are_errors() {
     assert!(solve_explicit(&full.model, None, &Limits::default()).unwrap_err().to_string().contains("not supported yet"));
     let mut ortho = Bar::new(ElementKind::Tet4, 4, ExplicitDynamic::new(end));
     let constants = Elastic::EngineeringConstants { young: [YOUNG; 3], poisson: [0.0; 3], shear: [YOUNG / 2.0; 3] };
-    ortho.model.replace_material("Mat", Material { elastic: Some(constants), density: Some(DENSITY), fluid: None }).unwrap();
+    ortho.model.replace_material("Mat", Material { elastic: Some(constants), density: Some(DENSITY), ..Material::default() }).unwrap();
     assert!(solve_explicit(&ortho.model, None, &Limits::default()).unwrap_err().to_string().contains("orthotropic"));
     let mut light = Bar::new(ElementKind::Tet4, 4, ExplicitDynamic::new(end));
     light
         .model
-        .replace_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: None, fluid: None })
+        .replace_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: None, ..Material::default() })
         .unwrap();
     assert!(solve_explicit(&light.model, None, &Limits::default()).is_err());
     // A node behind its wall at the start.
diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs
index b22806d..f9decff 100644
--- a/crates/fem/src/lib.rs
+++ b/crates/fem/src/lib.rs
@@ -20,6 +20,10 @@
 //! - [`face`]: consistent nodal forces of pressure on element faces;
 //! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass
 //!   shape vectors and element frequency bounds for explicit dynamics;
+//! - [`rotation`]: incrementally objective stress updates (Hughes–Winget Jaumann, Flanagan–Taylor
+//!   Green–Naghdi);
+//! - [`plasticity`]: J2 radial return with isotropic, kinematic and Johnson–Cook hardening,
+//!   Johnson–Cook damage;
 //! - [`dense`]: the small dense linear algebra inside an element.
 //!
 //! Conventions: element matrices are ordered node by node, each node's degrees of freedom in
@@ -38,7 +42,9 @@ pub mod face;
 pub mod jet;
 pub mod material;
 pub mod plane;
+pub mod plasticity;
 pub mod quadrature;
+pub mod rotation;
 pub mod shape;
 pub mod uniform;
 
diff --git a/crates/fem/src/plasticity.rs b/crates/fem/src/plasticity.rs
new file mode 100644
index 0000000..d667c99
--- /dev/null
+++ b/crates/fem/src/plasticity.rs
@@ -0,0 +1,484 @@
+//! Hypoelastic–plastic stress update of an isotropic material: von Mises (J2) plasticity with
+//! isotropic, linear kinematic or Johnson–Cook hardening, integrated by radial return, and
+//! Johnson–Cook damage.
+//!
+//! **Radial return** (J. C. Simo, T. J. R. Hughes, *Computational Inelasticity*, Springer 1998,
+//! §3.3, Box 3.1 and 3.2; R. D. Krieg, D. B. Krieg, J. Press. Vessel Technol. 99(4):510–515,
+//! 1977, for its accuracy): with shear modulus G and bulk modulus K, the pressure takes the
+//! volumetric increment elastically, p ← p + K Δθ, and the deviator the elastic trial
+//! s* = s + 2G Δe. With the back stress α (kinematic hardening), ξ* = s* − α and the trial
+//! von Mises stress q* = √(3/2)‖ξ*‖. If q* exceeds the current yield stress, the equivalent
+//! plastic strain increment Δε_p solves the consistency condition
+//! q* − 3G Δε_p − H_k Δε_p = σ_y(ε_p + Δε_p, Δε_p / Δt, T), and with the flow direction
+//! n = ξ*/‖ξ*‖: s = s* − 2G √(3/2) Δε_p n, α ← α + √(2/3) H_k Δε_p n (Prager's linear kinematic
+//! rule; H_k the kinematic modulus, so the uniaxial stress rises by H_k Δε_p either way).
+//! - linear hardening: Δε_p = (q* − σ_y)/(3G + H) in closed form;
+//! - a hardening table: piecewise linear, solved exactly segment by segment;
+//! - Johnson–Cook (rate and temperature dependent): a bracketed scalar Newton iteration with
+//!   bisection safeguard on [0, (q* − σ_y)/(3G)] (the root is in it because σ_y only grows with
+//!   Δε_p), capped at [`MAX_ITERATIONS`].
+//!
+//! **Johnson–Cook** flow stress σ_y = (A + B ε_pⁿ)(1 + C ln ε̇*)(1 − T*ᵐ) (G. R. Johnson,
+//! W. H. Cook, Proc. 7th Int. Symp. Ballistics, 1983) with ε̇* = ε̇_p / ε̇₀ taken as at least 1
+//! (the logarithm stays ≥ 0, which also removes its singularity at rest) and
+//! T* = (T − T_room)/(T_melt − T_room) clamped to [0, 1]. The temperature is that at the start
+//! of the increment; plastic work heats the material adiabatically, ΔT = χ q Δε_p / (ρ c_p).
+//!
+//! **Johnson–Cook damage** (G. R. Johnson, W. H. Cook, Eng. Fract. Mech. 21(1):31–48, 1985):
+//! D ← D + Δε_p / ε_f with ε_f = (D1 + D2 e^{D3 σ*})(1 + D4 ln ε̇*)(1 + D5 T*), σ* the mean
+//! stress over the von Mises stress after the increment. A fracture strain below
+//! [`MIN_FRACTURE_STRAIN`] counts as that (some published sets give negative values at high
+//! triaxiality, where the point fails at once). The point has failed at D ≥ 1.
+//!
+//! The strain increment the update takes is the **mechanical** one: a caller with initial or
+//! thermal strains subtracts their increment first.
+//!
+//! Stress 6-vectors hold tensor components (xx, yy, zz, xy, xz, yz), strain increments
+//! engineering shears.
+
+use modelcraft_model::{Damage, Elastic, Hardening, Plastic};
+
+use crate::material::von_mises;
+use crate::{Error, Result};
+
+/// Most iterations of the rate-dependent return mapping.
+pub const MAX_ITERATIONS: u32 = 100;
+
+/// The smallest Johnson–Cook fracture strain.
+pub const MIN_FRACTURE_STRAIN: f64 = 1e-6;
+
+/// The flow stress of a plastic material.
+#[derive(Debug, Clone, PartialEq)]
+pub enum Flow {
+    /// σ_y = y0 + h ε_p (isotropic) or a surface of radius y0 moving with modulus h (kinematic).
+    Linear {
+        y0: f64,
+        h: f64,
+        kinematic: bool,
+    },
+    /// Isotropic hardening by (plastic strain, yield stress) points, constant after the last.
+    Table(Vec<(f64, f64)>),
+    JohnsonCook {
+        a: f64,
+        b: f64,
+        n: f64,
+        c: f64,
+        m: f64,
+        rate0: f64,
+        room: f64,
+        melt: f64,
+    },
+}
+
+/// A hypoelastic, possibly plastic, isotropic material.
+#[derive(Debug, Clone, PartialEq)]
+pub struct J2 {
+    pub shear: f64,
+    pub bulk: f64,
+    pub flow: Option<Flow>,
+    /// Johnson–Cook damage constants D1–D5.
+    pub damage: Option<[f64; 5]>,
+    /// χ / c_p: the temperature rise per plastic work per mass (0 without heating).
+    pub heating: f64,
+}
+
+/// The state of a material point.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct PointState {
+    pub stress: [f64; 6],
+    pub back: [f64; 6],
+    pub plastic_strain: f64,
+    pub temperature: f64,
+    pub damage: f64,
+}
+
+impl PointState {
+    pub fn new(temperature: f64) -> PointState {
+        PointState { stress: [0.0; 6], back: [0.0; 6], plastic_strain: 0.0, temperature, damage: 0.0 }
+    }
+}
+
+/// What one update did.
+#[derive(Debug, Clone, Copy, Default, PartialEq)]
+pub struct Increment {
+    /// The equivalent plastic strain increment.
+    pub plastic_strain: f64,
+    /// Plastic work per volume, q Δε_p.
+    pub plastic_work: f64,
+}
+
+/// ‖s‖² of a symmetric tensor's 6-vector (tensor components).
+fn norm2(s: &[f64; 6]) -> f64 {
+    s[0] * s[0] + s[1] * s[1] + s[2] * s[2] + 2.0 * (s[3] * s[3] + s[4] * s[4] + s[5] * s[5])
+}
+
+impl Flow {
+    /// The flow stress at plastic strain `e`, plastic strain rate `rate` and temperature `t`, and
+    /// its derivative with respect to the plastic strain increment over an increment `dt` (for
+    /// the Newton iteration).
+    fn stress(&self, e: f64, rate: f64, t: f64, dt: f64) -> (f64, f64) {
+        match *self {
+            Flow::Linear { y0, h, kinematic } => {
+                if kinematic {
+                    (y0, 0.0)
+                } else {
+                    (y0 + h * e, h)
+                }
+            }
+            Flow::Table(ref points) => table(points, e),
+            Flow::JohnsonCook { a, b, n, c, m, rate0, room, melt } => {
+                let hardening = a + b * e.max(0.0).powf(n);
+                let d_hardening = if e > 0.0 {
+                    b * n * e.powf(n - 1.0)
+                } else if n < 1.0 && b > 0.0 {
+                    f64::INFINITY
+                } else {
+                    b * n
+                };
+                let ratio = rate / rate0;
+                let (rate_term, d_rate) = if ratio > 1.0 && dt > 0.0 { (1.0 + c * ratio.ln(), c / (rate * dt)) } else { (1.0, 0.0) };
+                let homologous = ((t - room) / (melt - room)).clamp(0.0, 1.0);
+                let thermal = if homologous > 0.0 { 1.0 - homologous.powf(m) } else { 1.0 };
+                (hardening * rate_term * thermal, (d_hardening * rate_term + hardening * d_rate) * thermal)
+            }
+        }
+    }
+
+    /// The kinematic hardening modulus.
+    fn kinematic(&self) -> f64 {
+        match *self {
+            Flow::Linear { h, kinematic: true, .. } => h,
+            _ => 0.0,
+        }
+    }
+
+    /// The rate sensitivity's reference rate and the homologous temperature (for damage).
+    fn rate_and_temperature(&self, t: f64) -> (f64, f64) {
+        match *self {
+            Flow::JohnsonCook { rate0, room, melt, .. } => (rate0, ((t - room) / (melt - room)).clamp(0.0, 1.0)),
+            _ => (1.0, 0.0),
+        }
+    }
+}
+
+/// The yield stress of a hardening table at plastic strain `e`, and its slope.
+fn table(points: &[(f64, f64)], e: f64) -> (f64, f64) {
+    for w in points.windows(2) {
+        if let [(e0, s0), (e1, s1)] = *w
+            && e < e1
+        {
+            let slope = (s1 - s0) / (e1 - e0);
+            return (s0 + slope * (e - e0).max(0.0), slope);
+        }
+    }
+    (points.last().map_or(0.0, |p| p.1), 0.0)
+}
+
+impl J2 {
+    /// The material of a model's isotropic elastic constants with its plasticity and damage.
+    pub fn from_model(elastic: &Elastic, plastic: Option<&Plastic>, damage: Option<&Damage>) -> Result<J2> {
+        let Elastic::Isotropic { young, poisson } = *elastic else {
+            return Err(Error::NotSupported("plasticity of an orthotropic material".into()));
+        };
+        let shear = young / (2.0 * (1.0 + poisson));
+        let bulk = young / (3.0 * (1.0 - 2.0 * poisson));
+        let mut heating = 0.0;
+        let flow = plastic.map(|p| match p {
+            Plastic::Bilinear { yield_stress, hardening_modulus, hardening } => {
+                Flow::Linear { y0: *yield_stress, h: *hardening_modulus, kinematic: *hardening == Hardening::Kinematic }
+            }
+            Plastic::Table { points } => Flow::Table(points.iter().map(|p| (p.strain, p.stress)).collect()),
+            Plastic::JohnsonCook { a, b, n, c, m, reference_strain_rate, room_temperature, melt_temperature, specific_heat, taylor_quinney } => {
+                heating = specific_heat.map_or(0.0, |cp| taylor_quinney / cp);
+                Flow::JohnsonCook {
+                    a: *a,
+                    b: *b,
+                    n: *n,
+                    c: *c,
+                    m: *m,
+                    rate0: *reference_strain_rate,
+                    room: *room_temperature,
+                    melt: *melt_temperature,
+                }
+            }
+        });
+        let damage = damage.map(|Damage::JohnsonCook { d1, d2, d3, d4, d5 }| [*d1, *d2, *d3, *d4, *d5]);
+        if !(shear > 0.0 && bulk > 0.0 && shear.is_finite() && bulk.is_finite()) {
+            return Err(Error::Invalid(format!("E = {young}, ν = {poisson} are not a stable material")));
+        }
+        Ok(J2 { shear, bulk, flow, damage, heating })
+    }
+
+    /// The temperature a point starts at.
+    pub fn initial_temperature(&self) -> f64 {
+        match self.flow {
+            Some(Flow::JohnsonCook { room, .. }) => room,
+            _ => 0.0,
+        }
+    }
+
+    /// Update `state` over an increment of length `dt` with the deviatoric strain increment
+    /// `dev` (engineering shears; its trace is ignored) and the volumetric strain increment
+    /// `vol`, at the current density `density` (for heating).
+    pub fn update(&self, state: &mut PointState, dev: &[f64; 6], vol: f64, dt: f64, density: f64) -> Increment {
+        let s = &mut state.stress;
+        let mean = (s[0] + s[1] + s[2]) / 3.0;
+        let tr = (dev[0] + dev[1] + dev[2]) / 3.0;
+        let g2 = 2.0 * self.shear;
+        // Trial deviator.
+        let mut trial = [s[0] - mean + g2 * (dev[0] - tr), s[1] - mean + g2 * (dev[1] - tr), s[2] - mean + g2 * (dev[2] - tr), 0.0, 0.0, 0.0];
+        for k in 3..6 {
+            trial[k] = s[k] + self.shear * dev[k];
+        }
+        let pressure = mean + self.bulk * vol;
+        let mut out = Increment::default();
+        if let Some(flow) = &self.flow {
+            let xi: [f64; 6] = std::array::from_fn(|k| trial[k] - state.back[k]);
+            let norm = norm2(&xi).sqrt();
+            let q = (1.5f64).sqrt() * norm;
+            let (yield_now, _) = flow.stress(state.plastic_strain, 0.0, state.temperature, dt);
+            if q > yield_now && norm > 0.0 {
+                let hk = flow.kinematic();
+                let de = self.plastic_increment(flow, state, q, yield_now, hk, dt);
+                let n: [f64; 6] = xi.map(|x| x / norm);
+                let radial = g2 * (1.5f64).sqrt() * de;
+                let shift = (2.0f64 / 3.0).sqrt() * hk * de;
+                for k in 0..6 {
+                    trial[k] -= radial * n[k];
+                    state.back[k] += shift * n[k];
+                }
+                state.plastic_strain += de;
+                let q_new = q - (3.0 * self.shear + hk) * de;
+                out = Increment { plastic_strain: de, plastic_work: q_new * de };
+                if self.heating > 0.0 && density > 0.0 {
+                    state.temperature += self.heating * out.plastic_work / density;
+                }
+            }
+        }
+        for k in 0..3 {
+            trial[k] += pressure;
+        }
+        state.stress = trial;
+        if let (Some(d), Some(flow)) = (self.damage, &self.flow)
+            && out.plastic_strain > 0.0
+        {
+            let q = von_mises(&state.stress);
+            let triaxiality = if q > 0.0 { pressure / q } else { 0.0 };
+            let (rate0, homologous) = flow.rate_and_temperature(state.temperature);
+            let rate = if dt > 0.0 { out.plastic_strain / dt / rate0 } else { 0.0 };
+            let log_rate = if rate > 1.0 { rate.ln() } else { 0.0 };
+            let fracture = (d[0] + d[1] * (d[2] * triaxiality).exp()) * (1.0 + d[3] * log_rate) * (1.0 + d[4] * homologous);
+            state.damage += out.plastic_strain / fracture.max(MIN_FRACTURE_STRAIN);
+        }
+        out
+    }
+
+    /// The plastic strain increment of a trial state that yields.
+    fn plastic_increment(&self, flow: &Flow, state: &PointState, q: f64, yield_now: f64, hk: f64, dt: f64) -> f64 {
+        let g3 = 3.0 * self.shear;
+        let (e0, t) = (state.plastic_strain, state.temperature);
+        match flow {
+            Flow::Linear { h, kinematic, .. } => (q - yield_now) / (g3 + if *kinematic { hk } else { *h }),
+            Flow::Table(points) => {
+                // Walk the segments from the current one: on a segment of slope h starting at
+                // (e_a, s_a), q − 3G Δ = s_a + h (e0 + Δ − e_a).
+                for w in points.windows(2) {
+                    let [(ea, sa), (eb, sb)] = *w else { continue };
+                    if eb <= e0 {
+                        continue;
+                    }
+                    let h = (sb - sa) / (eb - ea);
+                    let at = sa + h * (e0 - ea);
+                    let de = (q - at) / (g3 + h);
+                    if e0 + de <= eb {
+                        return de.max(0.0);
+                    }
+                }
+                // Beyond the table: perfectly plastic at the last stress.
+                let last = points.last().map_or(yield_now, |p| p.1);
+                ((q - last) / g3).max(0.0)
+            }
+            Flow::JohnsonCook { .. } => {
+                let residual = |de: f64| {
+                    let rate = if dt > 0.0 { de / dt } else { 0.0 };
+                    let (y, dy) = flow.stress(e0 + de, rate, t, dt);
+                    (q - g3 * de - y, -g3 - dy)
+                };
+                let (mut lo, mut hi) = (0.0, ((q - yield_now) / g3).max(0.0));
+                let mut de = 0.5 * hi;
+                for _ in 0..MAX_ITERATIONS {
+                    let (r, dr) = residual(de);
+                    if r > 0.0 {
+                        lo = de;
+                    } else {
+                        hi = de;
+                    }
+                    if r.abs() <= 1e-13 * q || hi - lo <= 1e-15 * hi {
+                        break;
+                    }
+                    let newton = de - r / dr;
+                    de = if newton.is_finite() && newton > lo && newton < hi { newton } else { 0.5 * (lo + hi) };
+                }
+                de
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use modelcraft_model::HardeningPoint;
+
+    const E: f64 = 200_000.0;
+    const NU: f64 = 0.3;
+
+    fn material(plastic: Plastic) -> J2 {
+        J2::from_model(&Elastic::Isotropic { young: E, poisson: NU }, Some(&plastic), None).unwrap()
+    }
+
+    /// Uniaxial stress at a constant strain rate: the axial strain increments by `de` while the
+    /// lateral strain increment is found by a secant iteration on the point update itself so
+    /// the lateral stresses vanish. Returns (axial strain, axial stress, plastic strain) per step.
+    fn uniaxial(m: &J2, steps: &[f64], dt: f64) -> Vec<(f64, f64, f64)> {
+        let mut s = PointState::new(m.initial_temperature());
+        let mut out = Vec::new();
+        let mut strain = 0.0;
+        for de in steps {
+            let run = |lateral: f64| {
+                let mut t = s;
+                let d = [*de, lateral * de, lateral * de, 0.0, 0.0, 0.0];
+                m.update(&mut t, &d, d[0] + d[1] + d[2], dt, 7.8e-9);
+                t
+            };
+            let (mut a, mut b) = (-NU, -0.5);
+            let (mut fa, mut fb) = (run(a).stress[1], run(b).stress[1]);
+            for _ in 0..100 {
+                if fb == fa || fb.abs() < 1e-12 {
+                    break;
+                }
+                let c = b - fb * (b - a) / (fb - fa);
+                (a, fa) = (b, fb);
+                b = c;
+                fb = run(b).stress[1];
+            }
+            s = run(b);
+            strain += de;
+            out.push((strain, s.stress[0], s.plastic_strain));
+        }
+        out
+    }
+
+    #[test]
+    fn bilinear_uniaxial_curve_and_unloading() {
+        let (y0, h) = (250.0, 2000.0);
+        for kinematic in [false, true] {
+            let hardening = if kinematic { Hardening::Kinematic } else { Hardening::Isotropic };
+            let m = material(Plastic::Bilinear { yield_stress: y0, hardening_modulus: h, hardening });
+            // Load to 1% strain, unload to −1%.
+            let mut steps = vec![1e-4; 100];
+            steps.extend(vec![-1e-4; 200]);
+            let curve = uniaxial(&m, &steps, 1e-3);
+            // Loading: σ = E ε below yield, then σ = y0 + H ε_p with ε = σ/E + ε_p, that is
+            // σ = y0 + E_t (ε − y0/E), E_t = E H/(E + H).
+            let et = E * h / (E + h);
+            let (e1, s1, p1) = curve[99];
+            let exact = y0 + et * (e1 - y0 / E);
+            assert!((s1 - exact).abs() < 1e-9 * exact, "{kinematic}: {s1} vs {exact}");
+            assert!((p1 - (s1 - y0) / h).abs() < 1e-12, "{p1}");
+            // Unloading is elastic until the reversed yield: isotropic at −σ1, kinematic at
+            // σ1 − 2 y0 (the surface keeps its size and moves).
+            let reverse = if kinematic { s1 - 2.0 * y0 } else { -s1 };
+            for (e, s, _) in &curve[100..] {
+                let elastic = s1 - E * (e1 - e);
+                if elastic > reverse + 1.0 {
+                    assert!((s - elastic).abs() < 1e-8 * s1, "{kinematic}: {s} vs elastic {elastic}");
+                }
+            }
+            let (_, s_end, _) = curve[299];
+            assert!(s_end < reverse, "{kinematic}: {s_end} vs {reverse}");
+            eprintln!("bilinear kinematic = {kinematic}: σ(1%) = {s1:.6} (exact {exact:.6}); reversed yield at {reverse:.3}; σ(−1%) = {s_end:.4}");
+        }
+    }
+
+    #[test]
+    fn table_and_johnson_cook_follow_their_flow_curves() {
+        let points = vec![
+            HardeningPoint { strain: 0.0, stress: 200.0 },
+            HardeningPoint { strain: 0.01, stress: 300.0 },
+            HardeningPoint { strain: 0.05, stress: 340.0 },
+        ];
+        let m = material(Plastic::Table { points });
+        let curve = uniaxial(&m, &vec![2e-4; 400], 1e-3);
+        for (_, s, p) in &curve {
+            if *p > 0.0 {
+                let want = table(&[(0.0, 200.0), (0.01, 300.0), (0.05, 340.0)], *p).0;
+                assert!((s - want).abs() < 1e-8 * want, "{s} vs {want} at {p}");
+            }
+        }
+        assert!(curve.last().unwrap().2 > 0.05, "past the table's end");
+        // Johnson–Cook at a constant plastic rate: the stress is the flow stress at that rate.
+        let (a, b, n, c) = (350.0, 275.0, 0.36, 0.022);
+        let jc = Plastic::JohnsonCook {
+            a,
+            b,
+            n,
+            c,
+            m: 1.0,
+            reference_strain_rate: 1.0,
+            room_temperature: 293.0,
+            melt_temperature: 1356.0,
+            specific_heat: None,
+            taylor_quinney: 0.9,
+        };
+        let m = material(jc);
+        let dt = 1e-6;
+        let curve = uniaxial(&m, &vec![1e-3; 200], dt);
+        let mut worst: f64 = 0.0;
+        for w in curve.windows(2).skip(20) {
+            let (_, s, p) = w[1];
+            let rate = (p - w[0].2) / dt;
+            let want = (a + b * p.powf(n)) * (1.0 + c * rate.ln());
+            worst = worst.max((s - want).abs() / want);
+        }
+        eprintln!("Johnson–Cook uniaxial: largest relative error of the flow stress {worst:.2e}");
+        assert!(worst < 1e-9, "{worst}");
+    }
+
+    #[test]
+    fn heating_softens_and_damage_accumulates() {
+        let jc = |cp: Option<f64>| Plastic::JohnsonCook {
+            a: 350.0,
+            b: 275.0,
+            n: 0.36,
+            c: 0.0,
+            m: 1.0,
+            reference_strain_rate: 1.0,
+            room_temperature: 293.0,
+            melt_temperature: 1356.0,
+            specific_heat: cp,
+            taylor_quinney: 0.9,
+        };
+        let el = Elastic::Isotropic { young: E, poisson: NU };
+        let damage = Damage::JohnsonCook { d1: 0.5, d2: 0.0, d3: 0.0, d4: 0.0, d5: 0.0 };
+        let cold = J2::from_model(&el, Some(&jc(None)), Some(&damage)).unwrap();
+        let hot = J2::from_model(&el, Some(&jc(Some(3.83e8))), Some(&damage)).unwrap();
+        let a = uniaxial(&cold, &vec![1e-3; 300], 1e-3);
+        let b = uniaxial(&hot, &vec![1e-3; 300], 1e-3);
+        assert!(b.last().unwrap().1 < a.last().unwrap().1);
+        // Constant fracture strain 0.5: D = ε_p / 0.5.
+        let mut s = PointState::new(293.0);
+        for _ in 0..300 {
+            cold.update(&mut s, &[1e-3, -5e-4, -5e-4, 0.0, 0.0, 0.0], 0.0, 1e-3, 7.8e-9);
+        }
+        assert!((s.damage - s.plastic_strain / 0.5).abs() < 1e-12, "{} vs {}", s.damage, s.plastic_strain);
+    }
+
+    #[test]
+    fn hostile_increments_stay_finite() {
+        let m = material(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic });
+        let mut s = PointState::new(0.0);
+        m.update(&mut s, &[1e3, -1e3, 0.0, 5e2, 0.0, 0.0], 0.0, 0.0, 0.0);
+        assert!(s.stress.iter().all(|x| x.is_finite()) && (von_mises(&s.stress) - 250.0).abs() < 1e-9 * 250.0);
+    }
+}
diff --git a/crates/fem/src/rotation.rs b/crates/fem/src/rotation.rs
new file mode 100644
index 0000000..414dcfe
--- /dev/null
+++ b/crates/fem/src/rotation.rs
@@ -0,0 +1,263 @@
+//! Incrementally objective stress updates for large deformation: the rotation that carries a
+//! stress (and every other tensor of the material state) from one configuration to the next.
+//!
+//! Increments are given as the gradient of the displacement increment with respect to the
+//! **midpoint configuration**, ΔG = ∂Δu/∂x_{n+½} (= Δt L at the midpoint). Its symmetric part is
+//! the strain increment ΔD and its skew part the spin increment ΔW.
+//!
+//! **Jaumann rate, Hughes–Winget integration** (T. J. R. Hughes, J. Winget, *Finite rotation
+//! effects in numerical integration of rate constitutive equations arising in large-deformation
+//! analysis*, IJNME 15(12):1862–1867, 1980): the stress is rotated with
+//! Q = (I − ½ΔW)⁻¹(I + ½ΔW) and the constitutive increment added: σ_{n+1} = Q σ_n Qᵀ + Δσ(ΔD).
+//! Because ΔG is taken on the midpoint configuration, a rigid rotation x_{n+1} = R x_n gives
+//! ΔG = 2(R − I)(R + I)⁻¹, which is skew (ΔD = 0), and its Cayley transform Q is R itself: a
+//! rigid rotation of any size rotates the stress exactly and strains nothing.
+//!
+//! **Green–Naghdi rate, Flanagan–Taylor integration** (D. P. Flanagan, L. M. Taylor, *An
+//! accurate numerical algorithm for stress integration with finite rotations*, CMAME
+//! 62(3):305–320, 1987): the rotation R of the polar decomposition F = V R is integrated
+//! incrementally instead of decomposing F. From V̇ = L V − V Ω with Ω = Ṙ Rᵀ and V symmetric,
+//! Ω = W + X where the skew X solves V X + X V = D V − V D, that is, with axial vectors,
+//! [(tr V) I − V] x = axial(D V − V D). R is advanced with the Cayley transform of Ω Δt, V with
+//! V̇, and the material is updated in the unrotated frame: σ̂ = Rᵀ σ R, d̂ = Rᵀ D R. A material
+//! whose axes are those of the undeformed body (orthotropic elasticity) is updated correctly in
+//! that frame.
+//!
+//! Tensors are symmetric 6-vectors in the order xx, yy, zz, xy, xz, yz; strain-like vectors carry
+//! engineering shears (γ = 2ε), stress-like ones tensor components.
+
+use crate::dense::inverse3;
+
+/// A 3 × 3 matrix, row-major.
+pub type Mat3 = [[f64; 3]; 3];
+
+pub const IDENTITY: Mat3 = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
+
+pub fn mul(a: &Mat3, b: &Mat3) -> Mat3 {
+    let mut c = [[0.0; 3]; 3];
+    for i in 0..3 {
+        for j in 0..3 {
+            c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
+        }
+    }
+    c
+}
+
+pub fn transpose(a: &Mat3) -> Mat3 {
+    [[a[0][0], a[1][0], a[2][0]], [a[0][1], a[1][1], a[2][1]], [a[0][2], a[1][2], a[2][2]]]
+}
+
+/// A symmetric tensor from its 6-vector of tensor components.
+pub fn tensor(s: &[f64; 6]) -> Mat3 {
+    [[s[0], s[3], s[4]], [s[3], s[1], s[5]], [s[4], s[5], s[2]]]
+}
+
+/// The 6-vector of tensor components of (the symmetric part of) a matrix.
+pub fn voigt(a: &Mat3) -> [f64; 6] {
+    [a[0][0], a[1][1], a[2][2], 0.5 * (a[0][1] + a[1][0]), 0.5 * (a[0][2] + a[2][0]), 0.5 * (a[1][2] + a[2][1])]
+}
+
+/// Q s Qᵀ for a stress-like 6-vector.
+pub fn rotate(q: &Mat3, s: &[f64; 6]) -> [f64; 6] {
+    voigt(&mul(&mul(q, &tensor(s)), &transpose(q)))
+}
+
+/// Q e Qᵀ for a strain-like 6-vector (engineering shears).
+pub fn rotate_strain(q: &Mat3, e: &[f64; 6]) -> [f64; 6] {
+    let half = [e[0], e[1], e[2], 0.5 * e[3], 0.5 * e[4], 0.5 * e[5]];
+    let r = rotate(q, &half);
+    [r[0], r[1], r[2], 2.0 * r[3], 2.0 * r[4], 2.0 * r[5]]
+}
+
+/// Q v for a vector.
+pub fn apply(q: &Mat3, v: &[f64; 3]) -> [f64; 3] {
+    [0, 1, 2].map(|i| q[i][0] * v[0] + q[i][1] * v[1] + q[i][2] * v[2])
+}
+
+/// The strain increment ΔD (engineering shears) and the spin increment ΔW of a displacement
+/// increment gradient.
+pub fn split(g: &Mat3) -> ([f64; 6], Mat3) {
+    let d = [g[0][0], g[1][1], g[2][2], g[0][1] + g[1][0], g[0][2] + g[2][0], g[1][2] + g[2][1]];
+    let mut w = [[0.0; 3]; 3];
+    for i in 0..3 {
+        for j in 0..3 {
+            w[i][j] = 0.5 * (g[i][j] - g[j][i]);
+        }
+    }
+    (d, w)
+}
+
+/// The Cayley transform (I − ½A)⁻¹(I + ½A) of a skew matrix: an exact rotation. `None` only
+/// for non-finite input (I − ½A is invertible for every real skew A).
+pub fn cayley(a: &Mat3) -> Option<Mat3> {
+    let mut minus = IDENTITY;
+    let mut plus = IDENTITY;
+    for i in 0..3 {
+        for j in 0..3 {
+            minus[i][j] -= 0.5 * a[i][j];
+            plus[i][j] += 0.5 * a[i][j];
+        }
+    }
+    let (inv, _) = inverse3(&minus)?;
+    let q = mul(&inv, &plus);
+    q.iter().flatten().all(|v| v.is_finite()).then_some(q)
+}
+
+/// The axial vector a of a skew matrix A (A v = a × v).
+pub fn axial(a: &Mat3) -> [f64; 3] {
+    [0.5 * (a[2][1] - a[1][2]), 0.5 * (a[0][2] - a[2][0]), 0.5 * (a[1][0] - a[0][1])]
+}
+
+/// The skew matrix of an axial vector.
+pub fn skew(a: &[f64; 3]) -> Mat3 {
+    [[0.0, -a[2], a[1]], [a[2], 0.0, -a[0]], [-a[1], a[0], 0.0]]
+}
+
+/// The polar rotation R and left stretch V (a 6-vector of tensor components) of the
+/// Flanagan–Taylor integration.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Polar {
+    pub rotation: Mat3,
+    pub stretch: [f64; 6],
+}
+
+impl Default for Polar {
+    fn default() -> Polar {
+        Polar { rotation: IDENTITY, stretch: [1.0, 1.0, 1.0, 0.0, 0.0, 0.0] }
+    }
+}
+
+impl Polar {
+    /// Advance R and V over an increment with displacement increment gradient `g` (on the
+    /// midpoint configuration). V̇ = L V − V Ω is integrated with the midpoint rule (a predictor
+    /// for V at the midpoint, then Ω and the increment from it), so V and R are second-order
+    /// accurate in the increment. Returns the increment's rotation (R_{n+1} = q R_n); `None`
+    /// when the increment is not finite or V degenerates.
+    pub fn advance(&mut self, g: &Mat3) -> Option<Mat3> {
+        let v0 = tensor(&self.stretch);
+        let (dv, omega) = Polar::rate(g, &v0)?;
+        let mut mid = v0;
+        for i in 0..3 {
+            for j in 0..3 {
+                mid[i][j] += 0.5 * dv[i][j];
+            }
+        }
+        let (dv, omega_mid) = Polar::rate(g, &mid).unwrap_or((dv, omega));
+        let mut vn = v0;
+        for i in 0..3 {
+            for j in 0..3 {
+                vn[i][j] += dv[i][j];
+            }
+        }
+        let q = cayley(&omega_mid)?;
+        self.rotation = mul(&q, &self.rotation);
+        self.stretch = voigt(&vn);
+        self.stretch.iter().all(|s| s.is_finite()).then_some(q)
+    }
+
+    /// The increments of V (L V − V Ω, times Δt) and the polar spin Ω Δt at stretch `v`.
+    fn rate(g: &Mat3, v: &Mat3) -> Option<(Mat3, Mat3)> {
+        let (d6, w) = split(g);
+        let d = tensor(&[d6[0], d6[1], d6[2], 0.5 * d6[3], 0.5 * d6[4], 0.5 * d6[5]]);
+        let (dv, vd) = (mul(&d, v), mul(v, &d));
+        let mut rhs = [[0.0; 3]; 3];
+        for i in 0..3 {
+            for j in 0..3 {
+                rhs[i][j] = dv[i][j] - vd[i][j];
+            }
+        }
+        let trace = v[0][0] + v[1][1] + v[2][2];
+        let mut m = v.map(|r| r.map(|x| -x));
+        for (i, row) in m.iter_mut().enumerate() {
+            row[i] += trace;
+        }
+        let (inv, _) = inverse3(&m)?;
+        let xs = skew(&apply(&inv, &axial(&rhs)));
+        let mut omega = w;
+        for i in 0..3 {
+            for j in 0..3 {
+                omega[i][j] += xs[i][j];
+            }
+        }
+        let (lv, vo) = (mul(g, v), mul(v, &omega));
+        let mut rate = [[0.0; 3]; 3];
+        for i in 0..3 {
+            for j in 0..3 {
+                rate[i][j] = lv[i][j] - vo[i][j];
+            }
+        }
+        Some((rate, omega))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn rotation_about(axis: [f64; 3], angle: f64) -> Mat3 {
+        let n = (axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]).sqrt();
+        let k = skew(&axis.map(|c| c / n));
+        let k2 = mul(&k, &k);
+        let mut r = IDENTITY;
+        for i in 0..3 {
+            for j in 0..3 {
+                r[i][j] += angle.sin() * k[i][j] + (1.0 - angle.cos()) * k2[i][j];
+            }
+        }
+        r
+    }
+
+    /// The midpoint gradient of the increment x → R x: 2 (R − I)(R + I)⁻¹.
+    fn rigid_increment(r: &Mat3) -> Mat3 {
+        let mut a = *r;
+        let mut b = *r;
+        for i in 0..3 {
+            a[i][i] -= 1.0;
+            b[i][i] += 1.0;
+        }
+        let (inv, _) = inverse3(&b).unwrap();
+        mul(&a, &inv).map(|row| row.map(|v| 2.0 * v))
+    }
+
+    #[test]
+    fn a_rigid_rotation_is_exact_with_both_rates() {
+        let r = rotation_about([0.3, -0.5, 0.8], 1.1);
+        let g = rigid_increment(&r);
+        let (d, w) = split(&g);
+        assert!(d.iter().all(|x| x.abs() < 1e-15), "{d:?}");
+        let q = cayley(&w).unwrap();
+        assert!(q.iter().flatten().zip(r.iter().flatten()).all(|(a, b)| (a - b).abs() < 1e-15), "{q:?} vs {r:?}");
+        let mut p = Polar::default();
+        p.advance(&g).unwrap();
+        assert!(p.rotation.iter().flatten().zip(r.iter().flatten()).all(|(a, b)| (a - b).abs() < 1e-14));
+        assert!(p.stretch.iter().zip([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]).all(|(a, b)| (a - b).abs() < 1e-14), "{:?}", p.stretch);
+        let s = [1.0, 2.0, 3.0, 0.4, 0.5, 0.6];
+        let rs = rotate(&q, &s);
+        // Invariants: trace and the norm.
+        let norm = |s: &[f64; 6]| s[..3].iter().map(|x| x * x).sum::<f64>() + 2.0 * s[3..].iter().map(|x| x * x).sum::<f64>();
+        assert!((rs[0] + rs[1] + rs[2] - 6.0).abs() < 1e-14 && (norm(&rs) - norm(&s)).abs() < 1e-13);
+    }
+
+    #[test]
+    fn polar_rotation_of_a_stretch_then_rotation() {
+        // Stretch along x by 1.5 in many small steps, then rotate 0.7 rad about z in many small
+        // steps: V = R diag(1.5, 1, 1) Rᵀ, R = rotation.
+        let mut p = Polar::default();
+        let n = 2000;
+        let step = 1.5f64.powf(1.0 / n as f64);
+        for _ in 0..n {
+            // x_{n+1} = diag(step) x_n; midpoint gradient 2 (s − 1)/(s + 1).
+            let mut g = [[0.0; 3]; 3];
+            g[0][0] = 2.0 * (step - 1.0) / (step + 1.0);
+            p.advance(&g).unwrap();
+        }
+        let r = rotation_about([0.0, 0.0, 1.0], 0.7 / n as f64);
+        for _ in 0..n {
+            p.advance(&rigid_increment(&r)).unwrap();
+        }
+        let want = rotation_about([0.0, 0.0, 1.0], 0.7);
+        let v = rotate(&want, &[1.5, 1.0, 1.0, 0.0, 0.0, 0.0]);
+        assert!(p.rotation.iter().flatten().zip(want.iter().flatten()).all(|(a, b)| (a - b).abs() < 1e-9), "{:?}", p.rotation);
+        assert!(p.stretch.iter().zip(v).all(|(a, b)| (a - b).abs() < 1e-6), "{:?} vs {v:?}", p.stretch);
+    }
+}
diff --git a/crates/fem/src/uniform.rs b/crates/fem/src/uniform.rs
index a577fb1..8cbbfa4 100644
--- a/crates/fem/src/uniform.rs
+++ b/crates/fem/src/uniform.rs
@@ -109,6 +109,65 @@ impl UniformStrain<8> {
     }
 }
 
+/// The hexahedron's volume as a form in its corner positions (in the model's hex8 node order):
+/// with C_IJK = ∫ ε_abc ∂_aN_I ∂_bN_J ∂_cN_K dξ (the determinant of the Jacobian expanded over
+/// the corners), which is antisymmetric in I, J, K, V = Σ_{I<J<K} C_IJK x_I · (x_J × x_K). The 24
+/// nonzero coefficients with I < J < K are ±1/12; this table holds them times 12 (the tests
+/// integrate them again).
+const VOLUME_FORM: [(usize, usize, usize, f64); 24] = [
+    (0, 1, 2, -1.0),
+    (0, 1, 3, -1.0),
+    (0, 1, 4, 1.0),
+    (0, 1, 5, 1.0),
+    (0, 2, 3, -1.0),
+    (0, 3, 4, -1.0),
+    (0, 3, 7, -1.0),
+    (0, 4, 5, -1.0),
+    (0, 4, 7, 1.0),
+    (1, 2, 3, -1.0),
+    (1, 2, 5, 1.0),
+    (1, 2, 6, 1.0),
+    (1, 4, 5, -1.0),
+    (1, 5, 6, -1.0),
+    (2, 3, 6, 1.0),
+    (2, 3, 7, 1.0),
+    (2, 5, 6, -1.0),
+    (2, 6, 7, -1.0),
+    (3, 4, 7, 1.0),
+    (3, 6, 7, -1.0),
+    (4, 5, 6, 1.0),
+    (4, 5, 7, 1.0),
+    (4, 6, 7, 1.0),
+    (5, 6, 7, 1.0),
+];
+
+impl UniformStrain<8> {
+    /// The integrated gradients and volume of a hexahedron in closed form: b_I = ∂V/∂x_I of the
+    /// volume form ([`VOLUME_FORM`]; D. P. Flanagan, T. Belytschko 1981 give it written out),
+    /// so b_I gathers C_IJK x_J × x_K, and V = ⅓ Σ_I b_I · x_I. A few hundred operations instead of a Gauss loop, for the
+    /// current configuration of large deformation. The element is refused when its volume is not
+    /// positive. `nodal_volume` is left at the mean V/8: an explicit solver keeps the lumped
+    /// masses of the undeformed element.
+    pub fn hex8_deformed(xyz: &[[f64; 3]; 8]) -> Result<Hex8> {
+        let mut b = [[0.0f64; 3]; 8];
+        let cross = |a: [f64; 3], b: [f64; 3]| [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
+        for (i, j, k, c) in VOLUME_FORM {
+            let (xi, xj, xk) = (xyz[i], xyz[j], xyz[k]);
+            for (target, v) in [(i, cross(xj, xk)), (j, cross(xk, xi)), (k, cross(xi, xj))] {
+                for d in 0..3 {
+                    b[target][d] += c * v[d];
+                }
+            }
+        }
+        let b = b.map(|v| v.map(|c| c / 12.0));
+        let volume: f64 = b.iter().zip(xyz).map(|(bi, x)| bi[0] * x[0] + bi[1] * x[1] + bi[2] * x[2]).sum::<f64>() / 3.0;
+        if !(volume > 0.0 && volume.is_finite() && b.iter().flatten().all(|v| v.is_finite())) {
+            return Err(degenerate("hex8"));
+        }
+        Ok(UniformStrain { b, volume, nodal_volume: [volume / 8.0; 8] })
+    }
+}
+
 impl UniformStrain<4> {
     /// The linear tetrahedron of corners `xyz` (the model's tet4 node order).
     pub fn tet4(xyz: &[[f64; 3]; 4]) -> Result<Tet4> {
@@ -281,6 +340,52 @@ mod tests {
         assert!(Hex8::hex8(&[[f64::NAN; 3]; 8]).is_err());
     }
 
+    #[test]
+    fn the_volume_form_is_the_integrated_one() {
+        let mut c = [[[0.0f64; 8]; 8]; 8];
+        for p in Rule::Hex(2).points() {
+            let s = shape::hex8(p.xi);
+            let d = s.derivatives();
+            for i in 0..8 {
+                for j in 0..8 {
+                    for k in 0..8 {
+                        let m = [[d[i][0], d[j][0], d[k][0]], [d[i][1], d[j][1], d[k][1]], [d[i][2], d[j][2], d[k][2]]];
+                        c[i][j][k] += p.weight * crate::dense::det3(&m);
+                    }
+                }
+            }
+        }
+        let mut table = Vec::new();
+        for i in 0..8 {
+            for j in i + 1..8 {
+                for k in j + 1..8 {
+                    if (12.0 * c[i][j][k]).abs() > 1e-9 {
+                        assert!((12.0 * c[i][j][k] - (12.0 * c[i][j][k]).round()).abs() < 1e-12);
+                        table.push((i, j, k, (12.0 * c[i][j][k]).round()));
+                    }
+                }
+            }
+        }
+        assert_eq!(table, VOLUME_FORM.to_vec());
+    }
+
+    #[test]
+    fn closed_form_gradients_equal_the_integrated_ones() {
+        let hex = distorted();
+        let (a, b) = (Hex8::hex8(&hex).unwrap(), Hex8::hex8_deformed(&hex).unwrap());
+        assert!((a.volume - b.volume).abs() < 1e-14 * a.volume, "{} vs {}", a.volume, b.volume);
+        for (x, y) in a.b.iter().flatten().zip(b.b.iter().flatten()) {
+            assert!((x - y).abs() < 1e-14, "{x} vs {y}");
+        }
+        let mut inverted = hex;
+        inverted.swap(0, 1);
+        inverted.swap(2, 3);
+        inverted.swap(4, 5);
+        inverted.swap(6, 7);
+        assert!(Hex8::hex8_deformed(&inverted).is_err());
+        assert!(Hex8::hex8_deformed(&[[f64::NAN; 3]; 8]).is_err());
+    }
+
     #[test]
     fn hourglass_vectors_are_orthogonal_to_linear_fields() {
         let hex = distorted();
diff --git a/crates/model/src/fluids.rs b/crates/model/src/fluids.rs
index e1a05d7..7932241 100644
--- a/crates/model/src/fluids.rs
+++ b/crates/model/src/fluids.rs
@@ -510,7 +510,7 @@ mod tests {
         }
         let mut solid = water();
         solid.fluid = None;
-        m.add_material("Steel", Material { elastic: None, density: Some(7.85e-9), fluid: None }).unwrap();
+        m.add_material("Steel", Material { elastic: None, density: Some(7.85e-9), ..Material::default() }).unwrap();
         assert!(m.add_fluid_domain("Dry", FluidDomain { material: "Steel".into(), ..channel() }).is_err());
         assert!(m.replace_material("Water", solid).is_err(), "a domain needs the fluid behaviour");
         assert!(matches!(m.remove_material("Water"), Err(Error::InUse(_))));
diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs
index 9d19fc0..fa64148 100644
--- a/crates/model/src/lib.rs
+++ b/crates/model/src/lib.rs
@@ -54,11 +54,11 @@ pub use constraints::{Coupling, CouplingKind, ReferencePoint};
 pub use fluids::{BoxFace, BoxFaces, FaceCondition, FluidBoundary, FluidDomain, Profile, SurfaceBoundary, SurfaceCondition};
 pub use geometry::{Classification, Entity, GeometryRef};
 pub use loads::{Distribution, Load, centrifugal_acceleration};
-pub use materials::{Elastic, Fluid, Material};
+pub use materials::{Damage, Elastic, Erosion, Fluid, Hardening, HardeningPoint, Material, Plastic};
 pub use mesh::{Element, ElementEdge, ElementFace, ElementKind, Node, check_coordinates};
 pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, SectionKind};
 pub use sets::Set;
-pub use steps::{BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, OutputVariable, Procedure, RigidWall, Step};
+pub use steps::{BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, OutputVariable, Procedure, RigidWall, Step, StressRate};
 pub use units::UnitSystem;
 
 /// The largest node or element label. Nastran's GRID and element ids stop at 99 999 999 (eight
diff --git a/crates/model/src/materials.rs b/crates/model/src/materials.rs
index d72713b..c3d7ef6 100644
--- a/crates/model/src/materials.rs
+++ b/crates/model/src/materials.rs
@@ -1,8 +1,21 @@
-//! Materials: linear elasticity, mass density and fluid behaviour.
+//! Materials: linear elasticity, mass density, fluid behaviour, and the inelastic behaviours of
+//! explicit dynamics: von Mises (J2) plasticity with a bilinear, tabulated or Johnson–Cook flow stress,
+//! Johnson–Cook damage, and element erosion criteria.
+//!
+//! The records are solver-neutral and named for what they describe, never after a product's
+//! keyword or material number; the deck writers map them to each format.
+//!
+//! - Johnson–Cook flow stress: G. R. Johnson, W. H. Cook, *A constitutive model and data for
+//!   metals subjected to large strains, high strain rates and high temperatures*, Proc. 7th Int.
+//!   Symp. on Ballistics, The Hague, 541–547, 1983:
+//!   σ_y = (A + B ε_pⁿ)(1 + C ln ε̇*)(1 − T*ᵐ), ε̇* = ε̇_p / ε̇₀, T* = (T − T_room)/(T_melt − T_room).
+//! - Johnson–Cook fracture: G. R. Johnson, W. H. Cook, *Fracture characteristics of three metals
+//!   subjected to various strains, strain rates, temperatures and pressures*, Eng. Fract. Mech.
+//!   21(1):31–48, 1985: ε_f = (D1 + D2 e^{D3 σ*})(1 + D4 ln ε̇*)(1 + D5 T*), D = Σ Δε_p / ε_f.
 
 use serde::{Deserialize, Serialize};
 
-use crate::{Error, Model, Result, context, insert_new, not_in_use, positive, remove_existing, replace_existing};
+use crate::{Error, Model, Result, context, finite, insert_new, not_in_use, positive, remove_existing, replace_existing};
 
 /// Names of the engineering constants' shear moduli, in order.
 pub(crate) const SHEAR_MODULI: [&str; 3] = ["G12", "G13", "G23"];
@@ -78,6 +91,179 @@ impl Fluid {
     }
 }
 
+/// Most points of a plastic hardening table.
+pub const MAX_TABLE_POINTS: usize = 1000;
+
+/// How the yield surface moves with plastic flow.
+#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)]
+#[serde(rename_all = "snake_case")]
+pub enum Hardening {
+    /// The yield surface grows (the yield stress rises with the equivalent plastic strain).
+    #[default]
+    Isotropic,
+    /// The yield surface translates (linear kinematic hardening of the back stress, Prager), so
+    /// reversed loading yields early (Bauschinger effect).
+    Kinematic,
+}
+
+/// One point of a hardening table: the yield stress at an equivalent plastic strain.
+#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
+#[serde(deny_unknown_fields)]
+pub struct HardeningPoint {
+    /// Equivalent plastic strain.
+    pub strain: f64,
+    /// Yield stress (von Mises equivalent).
+    pub stress: f64,
+}
+
+/// Rate-independent or Johnson–Cook von Mises plasticity (needs isotropic elasticity).
+#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
+#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)]
+pub enum Plastic {
+    /// Linear hardening: the yield stress `yield_stress + hardening_modulus · ε_p` (isotropic), or
+    /// a yield surface of radius `yield_stress` whose centre moves with modulus
+    /// `hardening_modulus` (kinematic). The hardening modulus is the slope of the uniaxial stress
+    /// against the plastic strain (0: perfectly plastic).
+    Bilinear {
+        yield_stress: f64,
+        #[serde(default)]
+        hardening_modulus: f64,
+        #[serde(default)]
+        hardening: Hardening,
+    },
+    /// Isotropic hardening by a table of yield stress against equivalent plastic strain, linear
+    /// between the points and constant after the last. The first point is at plastic strain 0;
+    /// the strains increase.
+    Table { points: Vec<HardeningPoint> },
+    /// Johnson–Cook flow stress (isotropic hardening): `a`, `b`, `n` the strain hardening, `c`
+    /// the strain-rate sensitivity relative to `reference_strain_rate` (rates below it count as
+    /// it), `m` the thermal softening between `room_temperature` and `melt_temperature`. With a
+    /// `specific_heat`, plastic work heats the material adiabatically (the fraction
+    /// `taylor_quinney` of it); without one the temperature stays at room temperature.
+    JohnsonCook {
+        a: f64,
+        b: f64,
+        n: f64,
+        #[serde(default)]
+        c: f64,
+        #[serde(default = "one")]
+        m: f64,
+        #[serde(default = "one")]
+        reference_strain_rate: f64,
+        room_temperature: f64,
+        melt_temperature: f64,
+        #[serde(default, skip_serializing_if = "Option::is_none")]
+        specific_heat: Option<f64>,
+        #[serde(default = "default_taylor_quinney")]
+        taylor_quinney: f64,
+    },
+}
+
+fn one() -> f64 {
+    1.0
+}
+
+fn default_taylor_quinney() -> f64 {
+    0.9
+}
+
+impl Plastic {
+    /// The name of the kind, for messages.
+    pub fn name(&self) -> &'static str {
+        match self {
+            Plastic::Bilinear { .. } => "bilinear plasticity",
+            Plastic::Table { .. } => "tabulated plasticity",
+            Plastic::JohnsonCook { .. } => "Johnson–Cook plasticity",
+        }
+    }
+
+    /// The initial yield stress (at zero plastic strain, the reference rate and room
+    /// temperature).
+    pub fn initial_yield(&self) -> f64 {
+        match self {
+            Plastic::Bilinear { yield_stress, .. } => *yield_stress,
+            Plastic::Table { points } => points.first().map_or(0.0, |p| p.stress),
+            Plastic::JohnsonCook { a, .. } => *a,
+        }
+    }
+
+    pub fn validate(&self) -> Result<()> {
+        let nonnegative = |what: &str, v: f64| {
+            if v >= 0.0 && v.is_finite() { Ok(v) } else { Err(Error::Invalid(format!("{what} must be a finite number ≥ 0, got {v}"))) }
+        };
+        match self {
+            Plastic::Bilinear { yield_stress, hardening_modulus, .. } => {
+                positive("yield stress", *yield_stress)?;
+                nonnegative("hardening modulus", *hardening_modulus)?;
+            }
+            Plastic::Table { points } => {
+                if points.is_empty() || points.len() > MAX_TABLE_POINTS {
+                    return Err(Error::Invalid(format!("a hardening table has 1 to {MAX_TABLE_POINTS} points, got {}", points.len())));
+                }
+                if points.first().is_some_and(|p| p.strain != 0.0) {
+                    return Err(Error::Invalid("a hardening table starts at plastic strain 0".into()));
+                }
+                for p in points {
+                    nonnegative("plastic strain", p.strain)?;
+                    positive("yield stress", p.stress)?;
+                }
+                if points.windows(2).any(|w| w.get(1).map(|b| b.strain) <= w.first().map(|a| a.strain)) {
+                    return Err(Error::Invalid("the plastic strains of a hardening table must increase".into()));
+                }
+            }
+            Plastic::JohnsonCook { a, b, n, c, m, reference_strain_rate, room_temperature, melt_temperature, specific_heat, taylor_quinney } => {
+                positive("A", *a)?;
+                nonnegative("B", *b)?;
+                nonnegative("n", *n)?;
+                nonnegative("C", *c)?;
+                positive("m", *m)?;
+                positive("reference strain rate", *reference_strain_rate)?;
+                finite("room temperature", *room_temperature)?;
+                finite("melt temperature", *melt_temperature)?;
+                if !(melt_temperature > room_temperature) {
+                    return Err(Error::Invalid(format!(
+                        "the melt temperature {melt_temperature} must be above the room temperature {room_temperature}"
+                    )));
+                }
+                if let Some(cp) = specific_heat {
+                    positive("specific heat", *cp)?;
+                }
+                if !(0.0..=1.0).contains(taylor_quinney) {
+                    return Err(Error::Invalid(format!("the Taylor–Quinney fraction is in [0, 1], got {taylor_quinney}")));
+                }
+            }
+        }
+        Ok(())
+    }
+}
+
+/// Damage that accumulates with plastic strain; a material point fails when it reaches 1.
+#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
+#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)]
+pub enum Damage {
+    /// Johnson–Cook fracture strain ε_f = (d1 + d2 e^{d3 σ*})(1 + d4 ln ε̇*)(1 + d5 T*), with σ*
+    /// the mean stress over the von Mises stress; ε̇* and T* as in Johnson–Cook plasticity (`d4`
+    /// and `d5` need it).
+    JohnsonCook {
+        d1: f64,
+        d2: f64,
+        d3: f64,
+        #[serde(default)]
+        d4: f64,
+        #[serde(default)]
+        d5: f64,
+    },
+}
+
+/// When elements of the material are removed from the analysis (eroded), besides failing by
+/// [`Damage`].
+#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
+#[serde(deny_unknown_fields)]
+pub struct Erosion {
+    /// The equivalent plastic strain at which an element erodes.
+    pub plastic_strain: f64,
+}
+
 /// A material.
 #[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
 #[serde(deny_unknown_fields)]
@@ -90,6 +276,16 @@ pub struct Material {
     /// How it flows, for fluid domains (needs a density).
     #[serde(default, skip_serializing_if = "Option::is_none")]
     pub fluid: Option<Fluid>,
+    /// Von Mises plasticity (explicit dynamics; the linear implicit solver uses the elastic
+    /// constants alone).
+    #[serde(default, skip_serializing_if = "Option::is_none")]
+    pub plastic: Option<Plastic>,
+    /// Damage to failure (explicit dynamics: a failed element erodes).
+    #[serde(default, skip_serializing_if = "Option::is_none")]
+    pub damage: Option<Damage>,
+    /// Erosion criteria (explicit dynamics).
+    #[serde(default, skip_serializing_if = "Option::is_none")]
+    pub erosion: Option<Erosion>,
 }
 
 impl Material {
@@ -109,6 +305,33 @@ impl Material {
         if let Some(d) = self.density {
             positive("density", d)?;
         }
+        if let Some(p) = &self.plastic {
+            if !matches!(self.elastic, Some(Elastic::Isotropic { .. })) {
+                return Err(Error::Invalid(format!("{} needs isotropic elastic constants", p.name())));
+            }
+            p.validate()?;
+        }
+        if let Some(Damage::JohnsonCook { d1, d2, d3, d4, d5 }) = self.damage {
+            for (what, v) in [("D1", d1), ("D2", d2), ("D3", d3), ("D4", d4), ("D5", d5)] {
+                finite(what, v)?;
+            }
+            match &self.plastic {
+                None => return Err(Error::Invalid("damage needs plasticity (it accumulates with the plastic strain)".into())),
+                Some(Plastic::JohnsonCook { .. }) => {}
+                Some(_) if d4 != 0.0 || d5 != 0.0 => {
+                    return Err(Error::Invalid(
+                        "the damage terms D4 and D5 need Johnson–Cook plasticity (its reference strain rate and temperatures)".into(),
+                    ));
+                }
+                Some(_) => {}
+            }
+        }
+        if let Some(e) = self.erosion {
+            positive("erosion plastic strain", e.plastic_strain)?;
+            if self.plastic.is_none() {
+                return Err(Error::Invalid("a plastic-strain erosion criterion needs plasticity".into()));
+            }
+        }
         Ok(())
     }
 }
@@ -188,4 +411,69 @@ mod tests {
         let zero_shear = Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0, 0.0, 1.0] };
         assert!(zero_shear.validate().unwrap_err().to_string().contains("G13"));
     }
+
+    fn steel() -> Material {
+        Material { density: Some(7.85e-9), ..iso(2.1e5, 0.3) }
+    }
+
+    fn jc() -> Plastic {
+        Plastic::JohnsonCook {
+            a: 792.0,
+            b: 510.0,
+            n: 0.26,
+            c: 0.014,
+            m: 1.03,
+            reference_strain_rate: 1.0,
+            room_temperature: 294.0,
+            melt_temperature: 1793.0,
+            specific_heat: Some(4.77e8),
+            taylor_quinney: 0.9,
+        }
+    }
+
+    #[test]
+    fn plasticity_damage_and_erosion_validate() {
+        let ok = |p: Plastic| Material { plastic: Some(p), ..steel() }.validate();
+        ok(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Kinematic }).unwrap();
+        ok(Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 250.0 }, HardeningPoint { strain: 0.1, stress: 300.0 }] }).unwrap();
+        ok(jc()).unwrap();
+        for bad in [
+            Plastic::Bilinear { yield_stress: 0.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic },
+            Plastic::Bilinear { yield_stress: 1.0, hardening_modulus: -1.0, hardening: Hardening::Isotropic },
+            Plastic::Table { points: vec![] },
+            Plastic::Table { points: vec![HardeningPoint { strain: 0.1, stress: 250.0 }] },
+            Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 250.0 }, HardeningPoint { strain: 0.0, stress: 300.0 }] },
+            Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: f64::NAN }] },
+            Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 1.0 }; MAX_TABLE_POINTS + 1] },
+        ] {
+            assert!(ok(bad.clone()).is_err(), "{bad:?}");
+        }
+        let mut cold = jc();
+        if let Plastic::JohnsonCook { melt_temperature, room_temperature, .. } = &mut cold {
+            *melt_temperature = *room_temperature;
+        }
+        assert!(ok(cold).is_err(), "melt above room");
+        // Plasticity needs isotropic elasticity; damage and erosion need plasticity.
+        let ortho = Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0; 3] };
+        assert!(Material { elastic: Some(ortho), plastic: Some(jc()), ..steel() }.validate().is_err());
+        let damage = Damage::JohnsonCook { d1: 0.05, d2: 3.44, d3: -2.12, d4: 0.002, d5: 0.61 };
+        Material { plastic: Some(jc()), damage: Some(damage), erosion: Some(Erosion { plastic_strain: 2.0 }), ..steel() }.validate().unwrap();
+        assert!(Material { damage: Some(damage), ..steel() }.validate().is_err());
+        let bilinear = Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic };
+        assert!(Material { plastic: Some(bilinear.clone()), damage: Some(damage), ..steel() }.validate().is_err(), "D4, D5 need Johnson–Cook");
+        let plain = Damage::JohnsonCook { d1: 0.5, d2: 0.0, d3: 0.0, d4: 0.0, d5: 0.0 };
+        Material { plastic: Some(bilinear), damage: Some(plain), ..steel() }.validate().unwrap();
+        assert!(Material { erosion: Some(Erosion { plastic_strain: 1.0 }), ..steel() }.validate().is_err());
+        // JSON: defaults fill in, unknown fields are refused.
+        let m: Material = serde_json::from_str(
+            r#"{"elastic":{"type":"isotropic","young":1,"poisson":0.3},"plastic":{"type":"johnson_cook","a":1,"b":2,"n":0.5,"room_temperature":0,"melt_temperature":100}}"#,
+        )
+        .unwrap();
+        assert!(
+            matches!(m.plastic, Some(Plastic::JohnsonCook { m, reference_strain_rate, taylor_quinney, .. }) if m == 1.0 && reference_strain_rate == 1.0 && taylor_quinney == 0.9)
+        );
+        assert!(serde_json::from_str::<Material>(r#"{"density":1,"plastic":{"type":"bilinear","yield_stress":1,"oops":2}}"#).is_err());
+        let text = serde_json::to_string(&Material { plastic: Some(jc()), damage: Some(damage), ..steel() }).unwrap();
+        assert_eq!(serde_json::from_str::<Material>(&text).unwrap(), Material { plastic: Some(jc()), damage: Some(damage), ..steel() });
+    }
 }
diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs
index 63ac4e3..72a260a 100644
--- a/crates/model/src/steps.rs
+++ b/crates/model/src/steps.rs
@@ -118,6 +118,30 @@ pub struct ExplicitDynamic {
     /// Artificial bulk viscosity.
     #[serde(default)]
     pub bulk_viscosity: BulkViscosity,
+    /// Large deformation (updated Lagrangian: gradients of the current configuration, stress
+    /// updated in rate form with an objective stress rate); `false` keeps the gradients of the
+    /// undeformed mesh (small displacements and rotations).
+    #[serde(default = "yes")]
+    pub nonlinear_geometry: bool,
+    /// The objective stress rate of large deformation.
+    #[serde(default)]
+    pub stress_rate: StressRate,
+}
+
+fn yes() -> bool {
+    true
+}
+
+/// The objective stress rate a large-deformation stress update uses.
+#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)]
+#[serde(rename_all = "snake_case")]
+pub enum StressRate {
+    /// Jaumann (co-rotational with the spin), integrated with the Hughes–Winget rotation.
+    #[default]
+    Jaumann,
+    /// Green–Naghdi (co-rotational with the rotation of the polar decomposition), integrated as
+    /// Flanagan and Taylor do.
+    GreenNaghdi,
 }
 
 fn default_dt_scale() -> f64 {
@@ -142,6 +166,8 @@ impl ExplicitDynamic {
             rigid_walls: Vec::new(),
             hourglass: Hourglass::default(),
             bulk_viscosity: BulkViscosity::default(),
+            nonlinear_geometry: true,
+            stress_rate: StressRate::default(),
         }
     }
 
diff --git a/crates/model/tests/definition.rs b/crates/model/tests/definition.rs
index 596815d..0178af5 100644
--- a/crates/model/tests/definition.rs
+++ b/crates/model/tests/definition.rs
@@ -56,8 +56,11 @@ fn defined() -> Model {
     m.add_set("all", Set::of_geometry(vec![GeometryRef::body("Bar").unwrap()])).unwrap();
     m.add_set("root", Set::of_geometry(vec![GeometryRef::face("Bar", "left").unwrap()])).unwrap();
     m.add_set("end", Set::of_geometry(vec![GeometryRef::face("Bar", "right").unwrap()])).unwrap();
-    m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9), fluid: None })
-        .unwrap();
+    m.add_material(
+        "Steel",
+        Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() },
+    )
+    .unwrap();
     m.add_section("bar", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
         .unwrap();
     m.add_reference_point("RP", ReferencePoint { xyz: [3.0, 0.5, 0.5], label: 0 }).unwrap();
@@ -121,7 +124,8 @@ fn the_check_names_what_is_missing() {
     m.add_set("side", Set { faces: vec![ElementFace { element: 1, face: 3 }], ..Set::default() }).unwrap();
     m.add_load("ghost-load", Load::Pressure { set: "ghost".into(), magnitude: 1.0 }).unwrap();
     m.add_load("twist", Load::Moment { target: Target::Set("end".into()), vector: [1.0, 0.0, 0.0], distribution: Distribution::Total }).unwrap();
-    m.add_material("Soft", Material { elastic: Some(Elastic::Isotropic { young: 2.1e11, poisson: 0.3 }), density: None, fluid: None }).unwrap();
+    m.add_material("Soft", Material { elastic: Some(Elastic::Isotropic { young: 2.1e11, poisson: 0.3 }), density: None, ..Material::default() })
+        .unwrap();
     m.add_section("again", Section { set: "side".into(), material: "Soft".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
         .unwrap();
     m.add_section(
@@ -327,8 +331,11 @@ fn a_total_force_on_a_face_is_the_same_on_every_mesh() {
     m.set_units(UnitSystem::MmTonneS);
     m.add_set("all", Set::of_geometry(vec![GeometryRef::body("Block").unwrap()])).unwrap();
     m.add_set("end", Set::of_geometry(vec![GeometryRef::face("Block", "end").unwrap()])).unwrap();
-    m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9), fluid: None })
-        .unwrap();
+    m.add_material(
+        "Steel",
+        Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() },
+    )
+    .unwrap();
     m.add_section("block", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
         .unwrap();
     let force = |distribution| Load::Force { target: Target::Set("end".into()), vector: [100.0, 0.0, -50.0], distribution };
@@ -417,8 +424,11 @@ fn plate(kind: SectionKind, x0: f64) -> Model {
     m.add_set("all", Set::of_geometry(vec![GeometryRef::face("Plate", "area").unwrap()])).unwrap();
     m.add_set("left", Set::of_geometry(vec![GeometryRef::new("Plate", Entity::Edge("left".into())).unwrap()])).unwrap();
     m.add_set("right", Set::of_geometry(vec![GeometryRef::new("Plate", Entity::Edge("right".into())).unwrap()])).unwrap();
-    m.add_material("Steel", Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9), fluid: None })
-        .unwrap();
+    m.add_material(
+        "Steel",
+        Material { elastic: Some(Elastic::Isotropic { young: 2.1e5, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() },
+    )
+    .unwrap();
     m.add_section("plate", Section { set: "all".into(), material: "Steel".into(), kind, formulation: Formulation::default() }).unwrap();
     m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("left".into()) }).unwrap();
     m.add_load("pull", Load::Traction { set: "right".into(), vector: [10.0, 0.0, 0.0] }).unwrap();
diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs
index 700c0f4..5f20529 100644
--- a/crates/solver/src/analysis.rs
+++ b/crates/solver/src/analysis.rs
@@ -12,7 +12,7 @@ use crate::constraints::{Budget, Reduction, coupling_rows};
 use crate::eigen::{Operator, Wanted, lanczos};
 use crate::ordering::nested_dissection;
 use crate::sparse::{Csr, NodeGraph, add_block, add_diagonal, add_scalar};
-use crate::system::{Mesh, loads, prescribed};
+use crate::system::{Mesh, inelastic_warnings, loads, prescribed};
 use crate::{Error, Limits, Result};
 
 /// Size and effort of a solution.
@@ -209,6 +209,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul
     let step = model.select_step(step, |p| matches!(p, Procedure::Static {}), "static")?;
     let mesh = Mesh::new(model, &resolved, limits.max_dofs)?;
     let mut warnings: Vec<String> = resolved.warnings.iter().map(ToString::to_string).collect();
+    warnings.extend(inelastic_warnings(model));
     let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice()));
     let order = nested_dissection(&mesh.xyz, &graph.adjacency);
     let setup = Setup { model, resolved: &resolved, mesh: &mesh, step };
@@ -336,6 +337,7 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt
     }
     let mesh = Mesh::new(model, &resolved, limits.max_dofs)?;
     let mut warnings: Vec<String> = resolved.warnings.iter().map(ToString::to_string).collect();
+    warnings.extend(inelastic_warnings(model));
     let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice()));
     let order = nested_dissection(&mesh.xyz, &graph.adjacency);
     let setup = Setup { model, resolved: &resolved, mesh: &mesh, step };
diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs
index 8dbac28..1179394 100644
--- a/crates/solver/src/system.rs
+++ b/crates/solver/src/system.rs
@@ -468,3 +468,17 @@ fn body(f: &mut [f64], mesh: &Mesh, e: &Elem, b: [f64; 3]) -> Result<()> {
     }
     Ok(())
 }
+
+/// The linear solver uses the elastic constants alone: a warning for each material in use whose
+/// plasticity, damage or erosion it leaves out.
+pub fn inelastic_warnings(model: &Model) -> Vec<String> {
+    let used: BTreeSet<&str> = model.sections().values().map(|s| s.material.as_str()).collect();
+    used.into_iter()
+        .filter(|name| model.material(name).is_some_and(|m| m.plastic.is_some() || m.damage.is_some() || m.erosion.is_some()))
+        .map(|name| {
+            format!(
+                "material `{name}`: the linear solver uses its elastic constants only (its plasticity, damage and erosion act in explicit dynamics)"
+            )
+        })
+        .collect()
+}
diff --git a/crates/solver/tests/common/mod.rs b/crates/solver/tests/common/mod.rs
index df5467a..5db0eb9 100644
--- a/crates/solver/tests/common/mod.rs
+++ b/crates/solver/tests/common/mod.rs
@@ -123,7 +123,7 @@ impl Grid {
     pub fn define(&mut self, elastic: Elastic, density: Option<f64>, formulation: Formulation, procedure: Procedure) {
         let all: Vec<u32> = self.model.elements().keys().copied().collect();
         self.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
-        self.model.add_material("Mat", Material { elastic: Some(elastic), density, fluid: None }).unwrap();
+        self.model.add_material("Mat", Material { elastic: Some(elastic), density, ..Material::default() }).unwrap();
         self.model.add_section("solid", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
         self.model.add_step(Step::new("Step", procedure), None).unwrap();
     }
diff --git a/crates/solver/tests/plane.rs b/crates/solver/tests/plane.rs
index c6b0121..d50f2fb 100644
--- a/crates/solver/tests/plane.rs
+++ b/crates/solver/tests/plane.rs
@@ -81,7 +81,9 @@ impl Grid {
     fn define(&mut self, kind: SectionKind, young: f64, nu: f64, density: Option<f64>, formulation: Formulation, procedure: Procedure) {
         let all: Vec<u32> = self.model.elements().keys().copied().collect();
         self.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
-        self.model.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young, poisson: nu }), density, fluid: None }).unwrap();
+        self.model
+            .add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young, poisson: nu }), density, ..Material::default() })
+            .unwrap();
         self.model.add_section("plane", Section { set: "all".into(), material: "Mat".into(), kind, formulation }).unwrap();
         self.model.add_step(Step::new("Step", procedure), None).unwrap();
     }
diff --git a/crates/solver/tests/verification.rs b/crates/solver/tests/verification.rs
index 32715ac..bb2ca1e 100644
--- a/crates/solver/tests/verification.rs
+++ b/crates/solver/tests/verification.rs
@@ -366,3 +366,26 @@ fn errors_name_the_problem() {
     assert!(solve_frequency(&g.model, Some("Step"), &FrequencyOptions::default(), &Limits::default()).is_err(), "a static step");
     assert!(solve_static(&g.model, Some("Step"), &Limits::default()).is_ok());
 }
+
+/// Plasticity is explicit dynamics' business: the linear static solver gives the elastic answer
+/// for a plastic material, bit for bit, and says it left the plasticity out.
+#[test]
+fn plastic_materials_solve_elastically_with_a_warning() {
+    let solve = |plastic: bool| {
+        let mut g = grid([4, 1, 1], ElementKind::Hex20, |q| [4.0 * q[0], q[1], q[2]]);
+        g.define(iso(1000.0, 0.3), None, full(), Procedure::Static {});
+        if plastic {
+            let mut m = g.model.material("Mat").unwrap().clone();
+            m.plastic = Some(modelcraft_model::Plastic::Bilinear { yield_stress: 1.0, hardening_modulus: 0.0, hardening: Default::default() });
+            g.model.replace_material("Mat", m).unwrap();
+        }
+        g.fixed("root", |q| near(q[0], 0.0));
+        g.face_set("tip", |q| near(q[0], 1.0));
+        g.load("p", Load::Pressure { set: "tip".into(), magnitude: -10.0 });
+        static_solve(&g)
+    };
+    let (a, b) = (solve(false), solve(true));
+    assert_eq!(a.displacement, b.displacement);
+    assert!(b.warnings.iter().any(|w| w.contains("material `Mat`") && w.contains("elastic constants only")), "{:?}", b.warnings);
+    assert!(!a.warnings.iter().any(|w| w.contains("elastic constants only")));
+}
diff --git a/docs/architecture.md b/docs/architecture.md
index 923f025..80c0278 100644
--- a/docs/architecture.md
+++ b/docs/architecture.md
@@ -21,11 +21,11 @@ crates in strictly lower layers.
 | L0 | `modelcraft-model` | Solver-neutral FE model: nodes and elements by label, their geometry tags, sets, materials, sections, couplings, BCs, loads and steps, the model check and the resolved view the writers read. Knows no file format and no geometry type. | started |
 | L1 | `modelcraft-geometry` | The only door to SolidCraft: opens `.solidcraft` designs (rebuilt from their timeline, bodies with persistent face names), B-rep bodies, measures, tessellation into a face-tagged `SurfaceMesh`; face-tagged 3MF files read into the same type without the kernel (docs/formats/tagged-mesh.md). No SolidCraft type in its API. | started |
 | L1 | `modelcraft-mesh` | Surface preparation (weld, watertightness, feature edges), robust predicates, Delaunay tetrahedralization and refinement, Voronoi-based optimization, hex meshing, quality metrics. | started |
-| L1 | `modelcraft-fem` | Element kernels shared by every structural solver: the element table (tet4, tet10, wedge6, hex8 full and incompatible modes, hex20 full and reduced), shape functions, quadrature, element stiffness, consistent and lumped mass, body and pressure loads, stresses and their extrapolation, isotropic and orthotropic elasticity, and the uniform-strain (one-point) hex8 and tet4 kinematics, hourglass shape vectors and element frequency bounds of the explicit solver. Pure functions, no solver loop. | started |
+| L1 | `modelcraft-fem` | Element kernels shared by every structural solver: the element table (tet4, tet10, wedge6, hex8 full and incompatible modes, hex20 full and reduced), shape functions, quadrature, element stiffness, consistent and lumped mass, body and pressure loads, stresses and their extrapolation, isotropic and orthotropic elasticity, the uniform-strain (one-point) hex8 and tet4 kinematics, hourglass shape vectors and element frequency bounds of the explicit solver, objective stress rates and J2/Johnson–Cook return mapping. Pure functions, no solver loop. | started |
 | L1 | `modelcraft-decks` | Abaqus input, Nastran bulk data and LS-DYNA keyword readers and writers, to and from `model`. Text out, no file system; `deck.write` saves it. | started |
 | L1 | `modelcraft-render` | Display data without a GPU dependency: the visible faces, element edges and nodes of a model, the orbit camera, a CPU rasterizer for headless snapshots, geometry surfaces with their face boundaries, deformed shapes and banded contour plots of nodal results, picking. | started |
 | L2 | `modelcraft-solver` | Implicit solver: assembly from the resolved model, boundary conditions, sparse direct solution, linear static and modal analysis, stresses and reactions. | started |
-| L2 | `modelcraft-explicit` | Explicit dynamics: central-difference time integration with lumped mass, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, rigid walls, energy account, histories and field frames. | started |
+| L2 | `modelcraft-explicit` | Explicit dynamics: central-difference time integration with lumped mass, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, rigid walls, large deformation with objective rates, plasticity, damage and erosion, energy account, histories and field frames. | started |
 | L2 | `modelcraft-lbm` | Lattice Boltzmann flow solver on the CPU: D3Q19 two-relaxation-time scheme, bounce-back and slip walls, Zou–He inlets and outlets, obstacles voxelised from the face-tagged surface with exact predicates, momentum-exchange forces per face. Uses `mesh` for the predicates and the boundary surface. | started |
 | L3 | `modelcraft-fsi` | Fluid–structure coupling between the L2 solvers (later). | planned |
 | L4 | `modelcraft-doc` | The document: geometry, mesh, model, results, history. | planned |
@@ -187,7 +187,8 @@ a kinematically coupled point (twist and shear stress).
 
 `modelcraft-explicit` integrates a model's `explicit_dynamic` steps in time (`solve.explicit`;
 `solve.history` and `solve.probe` read the stored histories and frames). Methods and papers are
-cited in the module docs and in docs/solvers/literature.md (Part H).
+cited in the module docs and in docs/solvers/literature.md (Parts H and J). Both structural solvers
+number their degrees of freedom with `modelcraft_fem::DofLayout` (docs/fem.md).
 
 1. **System.** `Model::resolve()`; hex8 sections with reduced integration become one-point
    uniform-strain hexahedra with Flanagan–Belytschko hourglass control, tet4 sections Bonet–Burton
@@ -210,11 +211,24 @@ cited in the module docs and in docs/solvers/literature.md (Part H).
    interval and field frames (displacement, velocity, element stress and strain) at the output
    interval, capped in size; kept in the session like the implicit results.
 
-The first version is geometrically linear with linear elastic materials. Verification
-(`crates/explicit/tests/verification.rs`, analytic only): stress waves in hex8 and tet4 bars
-(front at c·t, σ = ρ c v, work input), a bar striking a rigid wall (contact time 2L/c, rebound
-speed, plateau force ρ c v₀ A), the energy balance, bit-for-bit repeatability; in the engine
-(`solve.rs`), a cantilever's free-vibration period against `solve.frequency`.
+6. **Large deformation and materials (X2).** With `nonlinear_geometry` (the default) the
+   elements are updated Lagrangian: strain and spin increments on the midpoint configuration,
+   forces on the current one (hex8 gradients in closed form, `UniformStrain::hex8_deformed`). The
+   material state (stress, back stress, strain, stiffness hourglass forces) is rotated with the
+   Jaumann rate (Hughes–Winget) or the Green–Naghdi rate (Flanagan–Taylor), `fem::rotation`;
+   materials are hypoelastic or J2 plastic (bilinear isotropic/kinematic, hardening table,
+   Johnson–Cook with adiabatic heating) by radial return with Johnson–Cook damage,
+   `fem::plasticity`. Elements whose damage reaches 1 or whose plastic strain reaches the
+   material's erosion strain are eroded (reported with their energy; their nodes keep mass and
+   momentum). The energy account adds the plastic work and the eroded energy.
+
+Verification (`crates/explicit/tests/verification.rs` and `large_deformation.rs`, analytic
+only): stress waves in hex8 and tet4 bars (front at c·t, σ = ρ c v, work input), a bar striking a
+rigid wall (contact time 2L/c, rebound speed, plateau force ρ c v₀ A), the energy balance,
+bit-for-bit repeatability; Dienes' simple shear for both rates, a stressed element turned rigidly
+by 90°, the J2 uniaxial curve and unloading at large strain, damage failure at the fracture
+strain, and the Taylor bar against Taylor's 1948 estimate; in the engine (`solve.rs`), a
+cantilever's free-vibration period against `solve.frequency`.
 
 ## Lattice Boltzmann solver
 
diff --git a/docs/features.md b/docs/features.md
index eb2e996..85421da 100644
--- a/docs/features.md
+++ b/docs/features.md
@@ -5,7 +5,7 @@ Generated from the capability catalogue (`crates/engine/src/catalogue.rs`) by
 Every entry is a toolbar button. **Ready** entries run their command; **coming soon** entries are
 shown greyed with the roadmap milestone that brings them (ROADMAP.md).
 
-46 ready, 55 coming soon.
+48 ready, 53 coming soon.
 
 ## Geometry
 
@@ -56,9 +56,9 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md).
 | Panel | Capability | Status | Command |
 |---|---|---|---|
 | Materials | Material — Elastic (isotropic, orthotropic) material with density | ready | `model.material.create`, `model.material.edit`, `model.material.delete` |
-| Materials | Plasticity — Plastic hardening data | coming soon (M8) |  |
+| Materials | Plasticity — Von Mises plasticity: bilinear, hardening table or Johnson–Cook | ready | `model.material.plasticity` |
 | Materials | Hyperelastic — Hyperelastic material models | coming soon (M8) |  |
-| Materials | Damage — Damage initiation and evolution | coming soon (M8) |  |
+| Materials | Damage — Johnson–Cook damage and element erosion | ready | `model.material.damage` |
 | Materials | Thermal — Conductivity, specific heat and expansion | coming soon (M8) |  |
 | Materials | Equation of State — Polynomial, ideal gas, Mie–Grüneisen, JWL and tabulated equations of state | coming soon (M11) |  |
 | Sections | Section — Solid, shell, beam, rod or 2D (plane stress, plane strain, axisymmetric) section assigned to a set, with its element formulation | ready | `model.section.create`, `model.section.edit`, `model.section.delete` |
diff --git a/docs/parity.md b/docs/parity.md
index ae545bc..4352e52 100644
--- a/docs/parity.md
+++ b/docs/parity.md
@@ -44,7 +44,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
 | Mesh | Mesh editing (move, merge, split, collapse, remesh region), imported orphan meshes | M6 | yes | started (node/element add/remove) |
 | Model | Sets and surfaces bound to geometry names | M3 | yes | done |
 | Model | Materials: elastic (isotropic, orthotropic), density, unit-system check | M3 | | done |
-| Model | Materials: plasticity, hyperelasticity, damage, thermal | M8 | | planned |
+| Model | Materials: plasticity, hyperelasticity, damage, thermal | M8 | | started (von Mises plasticity: bilinear, table, Johnson–Cook; Johnson–Cook damage, erosion; written to Abaqus `*PLASTIC`, Nastran `MATS1`, LS-DYNA plastic materials where expressible) |
 | Model | Sections: solid, shell, beam profiles; assignment by set | M3 | | done |
 | Model | 2D sections: plane stress, plane strain, axisymmetric (XY plane, Y the axis); element edges in sets; centrifugal load | M4 | | done |
 | Model | Composite layups, material orientations, continuum shells | M8 | yes | planned |
@@ -62,7 +62,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
 | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) |
 | Decks | Readers keeping unknown keywords | M9 | | planned |
 | Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); reference points, kinematic and distributing couplings, moments; centrifugal loads) |
-| Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core: small-strain elastic, rigid walls; large deformation, plasticity and contact planned) |
+| Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; contact planned) |
 | Results | Deformed shape, contours, mode shapes, probes, reaction sums | M4 | | started (results in the session, `solve.probe`, reaction sums; display planned) |
 | Results | XY and path plots, cuts, display groups, free-body forces, animation | M9 | | planned |
 | Results | Readers for public result formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`) | M9 | | planned |
diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md
index b0c305f..68a3537 100644
--- a/docs/solvers/literature.md
+++ b/docs/solvers/literature.md
@@ -1061,3 +1061,33 @@ lattice Boltzmann package's code was opened.
   5(2):317–328, 1959, with the higher terms of A. S. Sangani, A. Acrivos, Int. J. Multiphase Flow
   8(4):343–360, 1982 (DOIs not checked). Benchmark (published numerical reference values, tagged,
   not part of the analytic gate): Schäfer & Turek 1996, 2D-1 (C8).
+
+## Part J: Large deformation and inelastic materials (X2), as implemented
+
+The papers and books `modelcraft_fem::{rotation, plasticity}` and `modelcraft-explicit`'s X2
+update are written from (cited again in the module docs):
+
+- Updated Lagrangian one-point elements, strain increments on the midpoint configuration:
+  T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, *Nonlinear Finite Elements for Continua and
+  Structures*, 2nd ed., Wiley 2014 (ch. 4–6); the closed-form integrated gradients b = ∂V/∂x of
+  the hexahedron, D. P. Flanagan, T. Belytschko, IJNME 17(5):679–706, 1981 (we tabulate the 24
+  coefficients of the volume form and check them against their exact integration).
+- Jaumann rate, incrementally objective: T. J. R. Hughes, J. Winget, IJNME 15(12):1862–1867,
+  1980 (Cayley rotation of the midpoint spin increment; a rigid rotation of any size is exact).
+- Green–Naghdi rate: D. P. Flanagan, L. M. Taylor, CMAME 62(3):305–320, 1987 (incremental polar
+  rotation from V̇ = L V − V Ω; we solve [(tr V) I − V] x = axial(D V − V D) and integrate V with
+  the midpoint rule).
+- Radial return for J2 with isotropic and linear kinematic (Prager) hardening: J. C. Simo, T. J.
+  R. Hughes, *Computational Inelasticity*, Springer 1998, §3.3; R. D. Krieg, D. B. Krieg, J.
+  Press. Vessel Technol. 99(4):510–515, 1977.
+- Johnson–Cook flow stress (rate clamped at ε̇* ≥ 1, T* in [0, 1], adiabatic heating with a
+  Taylor–Quinney fraction) and fracture strain D1–D5 with D = Σ Δε_p/ε_f: G. R. Johnson, W. H.
+  Cook 1983 (7th Int. Symp. Ballistics) and Eng. Fract. Mech. 21(1):31–48, 1985. The fracture
+  strain is floored at 10⁻⁶ (our choice).
+- Erosion: G. R. Johnson, R. A. Stryk, Int. J. Impact Eng. 5:411–421, 1987; T. Belytschko, J. I.
+  Lin, ibid. 5:111–127, 1987 (delete, keep nodal mass and momentum, report the energy).
+- Verification: J. K. Dienes, Acta Mech. 32:217–232, 1979 (simple shear, closed forms for both
+  rates); G. I. Taylor, Proc. R. Soc. A 194:289–299, 1948 (the yield stress from the final length
+  and undeformed length of a flat-ended projectile). Published experimental Taylor-bar lengths
+  (Wilkins & Guinan 1973) are not used: that needs the owner's sign-off for a tagged benchmark
+  suite.