diff --git a/crates/fem/src/constitutive.rs b/crates/fem/src/constitutive.rs
index d801387..529bd2e 100644
--- a/crates/fem/src/constitutive.rs
+++ b/crates/fem/src/constitutive.rs
@@ -36,11 +36,19 @@
 //!   §6.6, eqs. 6.165–6.169 and 6.193, in Kirchhoff form). The split is what the mixed (hybrid)
 //!   elements need: they take the isochoric part from [`MaterialModel::isochoric`] and the
 //!   volumetric one from [`MaterialModel::volumetric`] with the element's own dilatation.
+//! - **J2 plasticity** ([`crate::elastoplastic`]): Simo's multiplicative plasticity with the
+//!   exponential return map in the principal logarithmic strains (finite strain), or the
+//!   additive radial return (small strain, [`MaterialModel::small_strain`]); decoupled too, the
+//!   return deviatoric and the Hencky volume K/2 (ln J)².
+
+use std::sync::Arc;
 
 use fieldcraft_model::{Elastic, Hyperelastic};
 
 use crate::dense::det3;
+use crate::elastoplastic;
 use crate::material::{Elasticity, elasticity, stress};
+use crate::plasticity::J2;
 use crate::rotation::{IDENTITY, Mat3, mul, transpose, voigt};
 use crate::{Error, Result};
 
@@ -61,11 +69,16 @@ pub struct Response {
 #[derive(Debug, Clone, PartialEq)]
 pub enum MaterialModel {
     StVenantKirchhoff(Box<Elasticity>),
-    NeoHookean { shear: f64, bulk: f64 },
+    NeoHookean {
+        shear: f64,
+        bulk: f64,
+    },
+    /// Von Mises plasticity (shared: one per material, whatever the element count).
+    Plastic(Arc<J2>),
 }
 
 /// The 6 × 6 matrix of the symmetric fourth-order identity 𝕀 (engineering-shear columns).
-fn identity4() -> Tangent {
+pub(crate) fn identity4() -> Tangent {
     let mut i4 = [[0.0; 6]; 6];
     for k in 0..6 {
         i4[k][k] = if k < 3 { 1.0 } else { 0.5 };
@@ -138,23 +151,42 @@ impl MaterialModel {
         match self {
             MaterialModel::StVenantKirchhoff(_) => "St Venant–Kirchhoff",
             MaterialModel::NeoHookean { .. } => "neo-Hookean",
+            MaterialModel::Plastic(_) => "J2 plastic",
         }
     }
 
     /// History variables per integration point.
     pub fn history_len(&self) -> usize {
-        0
+        match self {
+            MaterialModel::Plastic(_) => elastoplastic::HISTORY,
+            _ => 0,
+        }
+    }
+
+    /// Whether the response depends on the path (the material keeps history), so the solver
+    /// treats it with the inelastic kinematics and predictors.
+    pub fn is_path_dependent(&self) -> bool {
+        self.history_len() > 0
+    }
+
+    /// The equivalent plastic strain in one point's history; `None` for an elastic material.
+    pub fn plastic_strain(&self, history: &[f64]) -> Option<f64> {
+        match self {
+            MaterialModel::Plastic(_) => Some(elastoplastic::value(history, elastoplastic::PLASTIC_STRAIN)),
+            _ => None,
+        }
     }
 
     /// Whether the strain energy splits into isochoric and volumetric parts (what mixed
     /// elements need).
     pub fn is_decoupled(&self) -> bool {
-        matches!(self, MaterialModel::NeoHookean { .. })
+        matches!(self, MaterialModel::NeoHookean { .. } | MaterialModel::Plastic(_))
     }
 
-    /// The Kirchhoff stress and spatial tangent at deformation gradient `f`, from the history
+    /// The Kirchhoff stress and spatial tangent at deformation gradient `f` and (for
+    /// temperature-dependent plasticity) temperature `temperature`, from the history
     /// `committed` (written, updated, to `updated`).
-    pub fn response(&self, f: &Mat3, committed: &[f64], updated: &mut [f64]) -> Result<Response> {
+    pub fn response(&self, f: &Mat3, temperature: Option<f64>, committed: &[f64], updated: &mut [f64]) -> Result<Response> {
         updated.iter_mut().zip(committed).for_each(|(u, c)| *u = *c);
         let j = jacobian(f)?;
         match self {
@@ -165,18 +197,24 @@ impl MaterialModel {
                 let t = push_forward(f);
                 Ok(Response { stress: std::array::from_fn(|a| (0..6).map(|b| t[a][b] * s[b]).sum()), tangent: sandwich(&t, d) })
             }
-            MaterialModel::NeoHookean { .. } => {
-                let iso = self.isochoric(f)?;
-                let (p, dp) = self.volumetric(j).ok_or_else(|| Error::Invalid("volumetric response".into()))?;
+            MaterialModel::NeoHookean { .. } | MaterialModel::Plastic(_) => {
+                let iso = self.isochoric(f, temperature, committed, updated)?;
+                let (p, dp) = self.volumetric(j, false).ok_or_else(|| Error::Invalid("volumetric response".into()))?;
                 Ok(add_pressure(&iso, j, p, j * dp))
             }
         }
     }
 
-    /// The isochoric part of a decoupled material (τ_iso = dev τ̄ and its tangent).
-    pub fn isochoric(&self, f: &Mat3) -> Result<Response> {
-        let &MaterialModel::NeoHookean { shear: mu, .. } = self else {
-            return Err(Error::NotSupported(format!("the isochoric–volumetric split of a {} material", self.name())));
+    /// The isochoric part of a decoupled material (τ_iso = dev τ̄ and its tangent), from the
+    /// history `committed`; the history that changes is written to `updated`, which the caller
+    /// starts as a copy of `committed`.
+    pub fn isochoric(&self, f: &Mat3, temperature: Option<f64>, committed: &[f64], updated: &mut [f64]) -> Result<Response> {
+        let mu = match self {
+            MaterialModel::NeoHookean { shear, .. } => *shear,
+            MaterialModel::Plastic(j2) => return elastoplastic::finite_isochoric(j2, f, temperature, committed, updated),
+            MaterialModel::StVenantKirchhoff(_) => {
+                return Err(Error::NotSupported(format!("the isochoric–volumetric split of a {} material", self.name())));
+            }
         };
         let j = jacobian(f)?;
         let scale = j.powf(-2.0 / 3.0);
@@ -199,13 +237,29 @@ impl MaterialModel {
         Ok(Response { stress: tau, tangent: c })
     }
 
-    /// U′(J) and U″(J) of a decoupled material's volumetric energy.
-    pub fn volumetric(&self, j: f64) -> Option<(f64, f64)> {
-        match *self {
-            MaterialModel::NeoHookean { bulk, .. } => Some((bulk * (j - 1.0), bulk)),
+    /// U′(J) and U″(J) of a decoupled material's volumetric energy at dilatation J (with
+    /// `small`, J = 1 + tr ε of small strain, where every material's is K/2 (J − 1)²).
+    pub fn volumetric(&self, j: f64, small: bool) -> Option<(f64, f64)> {
+        match self {
+            MaterialModel::NeoHookean { bulk, .. } => Some((bulk * (j - 1.0), *bulk)),
+            MaterialModel::Plastic(j2) => elastoplastic::volumetric(j2, j, small),
             MaterialModel::StVenantKirchhoff(_) => None,
         }
     }
+
+    /// Small strain: the stress and its tangent at the (mechanical) strain `strain`
+    /// (engineering shears), from the history `committed` (written, updated, to `updated`);
+    /// `bulk` false keeps the deviatoric part only (what a mixed element takes). St
+    /// Venant–Kirchhoff is linear elasticity here; hyperelastic materials have no small-strain
+    /// form.
+    pub fn small_strain(&self, strain: &[f64; 6], temperature: Option<f64>, committed: &[f64], updated: &mut [f64], bulk: bool) -> Result<Response> {
+        updated.iter_mut().zip(committed).for_each(|(u, c)| *u = *c);
+        match self {
+            MaterialModel::Plastic(j2) => elastoplastic::small(j2, strain, temperature, committed, updated, bulk),
+            MaterialModel::StVenantKirchhoff(d) if bulk => Ok(Response { stress: stress(d, strain), tangent: **d }),
+            _ => Err(Error::NotSupported(format!("the {} material at small strain{}", self.name(), if bulk { "" } else { " in mixed elements" }))),
+        }
+    }
 }
 
 /// An isochoric response with a pressure p acting on dilatation J: τ + J p I, and the tangent
@@ -269,7 +323,7 @@ mod tests {
     /// The tangent's consistency: dS = C : dE with C the pull-back of c, checked by central
     /// differences of S along several directions dF (dE = sym(Fᵀ dF), engineering shears).
     fn check_tangent(m: &MaterialModel, f: &Mat3) {
-        let r = m.response(f, &[], &mut []).unwrap();
+        let r = m.response(f, None, &[], &mut []).unwrap();
         let (inv, _) = inverse3(f).unwrap();
         let back = sandwich(&push_forward(&inv), &r.tangent);
         let scale = back.iter().flatten().fold(0.0f64, |a, v| a.max(v.abs()));
@@ -280,7 +334,7 @@ mod tests {
             df[(k + 1) % 3][(k + 2) % 3] += 0.3;
             let at = |s: f64| {
                 let fp: Mat3 = std::array::from_fn(|i| std::array::from_fn(|j| f[i][j] + s * df[i][j]));
-                pull_back_stress(&fp, &m.response(&fp, &[], &mut []).unwrap().stress)
+                pull_back_stress(&fp, &m.response(&fp, None, &[], &mut []).unwrap().stress)
             };
             let (sp, sm) = (at(h), at(-h));
             let ftdf = mul(&transpose(f), &df);
@@ -300,7 +354,7 @@ mod tests {
             check_tangent(&m, &deformation());
             check_tangent(&m, &IDENTITY);
             // Symmetric (hyperelastic): major symmetry of c.
-            let c = m.response(&deformation(), &[], &mut []).unwrap().tangent;
+            let c = m.response(&deformation(), None, &[], &mut []).unwrap().tangent;
             for a in 0..6 {
                 for b in 0..6 {
                     assert!((c[a][b] - c[b][a]).abs() < 1e-9 * 5000.0, "{a} {b}");
@@ -317,7 +371,7 @@ mod tests {
         let d = isotropic(e, nu).unwrap();
         let neo = MaterialModel::from_model(&Elastic::Isotropic { young: e, poisson: nu }, Some(&Hyperelastic::NeoHookean {})).unwrap();
         for m in [MaterialModel::StVenantKirchhoff(Box::new(d)), neo] {
-            let c = m.response(&IDENTITY, &[], &mut []).unwrap().tangent;
+            let c = m.response(&IDENTITY, None, &[], &mut []).unwrap().tangent;
             for a in 0..6 {
                 for b in 0..6 {
                     assert!((c[a][b] - d[a][b]).abs() < 1e-9 * e, "{}: {a} {b}", m.name());
@@ -329,7 +383,7 @@ mod tests {
             let k = crate::rotation::skew(&axis);
             let k2 = mul(&k, &k);
             let rot: Mat3 = std::array::from_fn(|i| std::array::from_fn(|j| IDENTITY[i][j] + s * k[i][j] + (1.0 - cth) * k2[i][j]));
-            let tau = m.response(&rot, &[], &mut []).unwrap().stress;
+            let tau = m.response(&rot, None, &[], &mut []).unwrap().stress;
             assert!(tau.iter().all(|v| v.abs() < 1e-12 * e), "{}: {tau:?}", m.name());
         }
     }
@@ -343,14 +397,14 @@ mod tests {
         let (lam, l) = (1.5f64, 0.85f64);
         let f = [[lam, 0.0, 0.0], [0.0, l, 0.0], [0.0, 0.0, l]];
         let j = lam * l * l;
-        let tau = m.response(&f, &[], &mut []).unwrap().stress;
+        let tau = m.response(&f, None, &[], &mut []).unwrap().stress;
         let i1 = lam * lam + 2.0 * l * l;
         let sxx = mu / j.powf(5.0 / 3.0) * (lam * lam - i1 / 3.0) + kappa * (j - 1.0);
         let syy = mu / j.powf(5.0 / 3.0) * (l * l - i1 / 3.0) + kappa * (j - 1.0);
         assert!((tau[0] / j - sxx).abs() < 1e-10 * mu && (tau[1] / j - syy).abs() < 1e-10 * mu, "{tau:?}");
         let e = logarithmic_strain(&f).unwrap();
         assert!((e[0] - lam.ln()).abs() < 1e-14 && (e[1] - l.ln()).abs() < 1e-14 && e[3].abs() < 1e-14, "{e:?}");
-        assert!(m.response(&[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0]], &[], &mut []).is_err(), "inverted");
-        assert!(MaterialModel::StVenantKirchhoff(Box::new(isotropic(1.0, 0.3).unwrap())).isochoric(&IDENTITY).is_err());
+        assert!(m.response(&[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0]], None, &[], &mut []).is_err(), "inverted");
+        assert!(MaterialModel::StVenantKirchhoff(Box::new(isotropic(1.0, 0.3).unwrap())).isochoric(&IDENTITY, None, &[], &mut []).is_err());
     }
 }
diff --git a/crates/fem/src/elastoplastic.rs b/crates/fem/src/elastoplastic.rs
new file mode 100644
index 0000000..e9af4f0
--- /dev/null
+++ b/crates/fem/src/elastoplastic.rs
@@ -0,0 +1,407 @@
+//! Von Mises (J2) elastoplasticity at a material point of implicit (Newton) analysis: the
+//! stress after an increment and its consistent tangent, at finite strain and at small strain.
+//! Both kinematics call the radial return of [`crate::plasticity`], the update explicit
+//! dynamics uses, so the yield surface, the hardening (bilinear isotropic or linear kinematic,
+//! tabulated, Johnson–Cook) and the plastic strain increment of a trial state are one
+//! implementation for both solvers.
+//!
+//! **Finite strain** (J. C. Simo, *Algorithms for static and dynamic multiplicative plasticity
+//! that preserve the classical return mapping schemes of the infinitesimal theory*, CMAME
+//! 99:61–112, 1992; Simo & Hughes 1998, ch. 9): the multiplicative split F = Fₑ Fₚ with the
+//! elastic left Cauchy–Green tensor bₑ = F Cₚ⁻¹ Fᵀ and a Hencky (quadratic in the logarithmic
+//! strain) elastic energy, so the Kirchhoff stress is τ = K ln J I + 2G dev εₑ with
+//! εₑ = ½ ln bₑ. The history is Cₚ⁻¹ (stored as Cₚ⁻¹ − I, so the virgin state is zeros) and the
+//! equivalent plastic strain. An increment freezes Cₚ⁻¹ at its committed value: the trial
+//! bₑᵗʳ = F Cₚ⁻¹ Fᵀ, its principal stretches and directions, the trial logarithmic strains
+//! ε_A = ½ ln λ_A², and then the **return map of the infinitesimal theory on the principal
+//! values** (the exponential map integrates the flow rule exactly in the principal frame, which
+//! the isotropy of the model keeps fixed during the return). The update writes back
+//! bₑ = Σ exp(2 εₑ,A) n_A ⊗ n_A and Cₚ⁻¹ = F⁻¹ bₑ F⁻ᵀ, which keeps det Cₚ = 1 (plastic flow
+//! is isochoric).
+//!
+//! **Its tangent.** With Cₚ⁻¹ frozen, the Lie derivative of bₑᵗʳ vanishes, so for a rate of
+//! deformation d the trial strain changes by dεᵗʳ = ½ (∂ ln b/∂b) : (d b + b d), and
+//! L_v τ = 𝔻ᵉᵖ : dεᵗʳ − (d τ + τ d), with 𝔻ᵉᵖ the algorithmic tangent of the return map
+//! ([`crate::plasticity::J2::tangent`]). In the principal frame of b the derivative of the
+//! logarithm is the Daleckii–Krein divided difference, so dεᵗʳ_AB = θ_AB d_AB with θ_AA = 1 and
+//! θ_AB = ½ (λ²_A + λ²_B) ln(λ²_A/λ²_B)/(λ²_A − λ²_B) = artanh(r)/r, r = (λ²_A − λ²_B)/(λ²_A + λ²_B),
+//! smooth through equal stretches (its series 1 + r²/3 + r⁴/5 + … near r = 0), so uniaxial and
+//! undeformed states need no special case (Simo 1992, §3.3, gives the same tangent in spectral
+//! form; R. Bhatia, *Matrix Analysis*, Springer 1997, §V.3, for the divided differences). The
+//! tangent is symmetric for isotropic hardening: 𝔻ᵉᵖ and the θ scaling share the principal
+//! frame.
+//!
+//! **Kinematic hardening at finite strain** is refused: a back stress is not coaxial with bₑ, so
+//! the principal-frame return no longer holds, and it needs a formulation of its own (for
+//! example the Lagrangian logarithmic strain space of C. Miehe, N. Apel, M. Lambrecht, CMAME
+//! 191:5383–5425, 2002). It works with small strain.
+//!
+//! **Small strain** (Simo & Hughes 1998, §3.3): ε = εₑ + εₚ additively, the history the plastic
+//! strain εₚ, the back stress α and the equivalent plastic strain; the trial εₑ = ε − εₚ,ₙ goes
+//! through the same return, and εₚ advances by √(3/2) Δε̄ₚ n.
+//!
+//! **Mixed (hybrid) elements** take the isochoric part (the deviatoric return, without K) and the
+//! volumetric energy U(J) = K/2 (ln J)², U′ = K ln J / J (finite strain), or U = K/2 (J − 1)² with
+//! J the 1 + tr ε of small strain: J2 plasticity is pressure insensitive and isochoric, so the
+//! split is exact.
+//!
+//! **History layout** per integration point ([`HISTORY`] values): 0–5 Cₚ⁻¹ − I (finite strain,
+//! tensor components) or εₚ (small strain, engineering shears), 6–11 the back stress α, 12 the
+//! equivalent plastic strain.
+
+use crate::constitutive::{Response, Tangent};
+use crate::dense::{inverse3, symmetric_eigen};
+use crate::plasticity::{J2, PointState, Return};
+use crate::rotation::{IDENTITY, Mat3, mul, tensor, transpose, voigt};
+use crate::{Error, Result};
+
+/// History values per integration point.
+pub const HISTORY: usize = 13;
+
+/// Where the equivalent plastic strain is in a point's history.
+pub const PLASTIC_STRAIN: usize = 12;
+
+/// History value `i` (0 when the slice is short).
+pub(crate) fn value(h: &[f64], i: usize) -> f64 {
+    h.get(i).copied().unwrap_or(0.0)
+}
+
+/// Cₚ⁻¹ of a finite-strain history (stored as Cₚ⁻¹ − I).
+fn plastic_stretch(h: &[f64]) -> Mat3 {
+    let stored = tensor(&std::array::from_fn(|k| value(h, k)));
+    std::array::from_fn(|i| std::array::from_fn(|j| IDENTITY[i][j] + stored[i][j]))
+}
+
+/// Store Cₚ⁻¹ as Cₚ⁻¹ − I in the first six history values.
+fn store_plastic_stretch(cp: &Mat3, h: &mut [f64]) {
+    let v = voigt(cp);
+    for (k, slot) in h.iter_mut().take(6).enumerate() {
+        *slot = v[k] - if k < 3 { 1.0 } else { 0.0 };
+    }
+}
+
+/// artanh(r)/r, smooth through r = 0.
+fn log_divided_difference(r: f64) -> f64 {
+    if r.abs() < 1e-4 {
+        let r2 = r * r;
+        1.0 + r2 * (1.0 / 3.0 + r2 * (1.0 / 5.0 + r2 / 7.0))
+    } else {
+        r.atanh() / r
+    }
+}
+
+/// The state a return starts from: the committed back stress and plastic strain, at the point's
+/// temperature (Johnson–Cook softening; room temperature without one).
+fn start(j2: &J2, committed: &[f64], temperature: Option<f64>) -> PointState {
+    let mut s = PointState::new(temperature.unwrap_or_else(|| j2.initial_temperature()));
+    s.back = std::array::from_fn(|k| value(committed, 6 + k));
+    s.plastic_strain = value(committed, PLASTIC_STRAIN);
+    s
+}
+
+fn singular() -> Error {
+    Error::Degenerate("the plastic material point has no tangent: its hardening falls faster than −3G (a softening material localises)".into())
+}
+
+/// The deviatoric trial stress 2G dev ε of a strain 6-vector (engineering shears), as tensor
+/// components.
+fn trial_deviator(j2: &J2, e: &[f64; 6]) -> [f64; 6] {
+    let mean = (e[0] + e[1] + e[2]) / 3.0;
+    let g2 = 2.0 * j2.shear;
+    [g2 * (e[0] - mean), g2 * (e[1] - mean), g2 * (e[2] - mean), j2.shear * e[3], j2.shear * e[4], j2.shear * e[5]]
+}
+
+/// Finite strain: the isochoric (deviatoric) Kirchhoff stress and its spatial tangent at
+/// deformation gradient `f`, from the history `committed` (written, updated, to `updated`). Add
+/// the volumetric part with [`volumetric`] (or the mixed element's pressure) for the whole
+/// response.
+pub fn finite_isochoric(j2: &J2, f: &Mat3, temperature: Option<f64>, committed: &[f64], updated: &mut [f64]) -> Result<Response> {
+    if let Some(crate::plasticity::Flow::Linear { kinematic: true, .. }) = j2.flow {
+        return Err(Error::NotSupported(
+            "kinematic hardening at finite strain (its back stress is not coaxial with the elastic stretch; use nonlinear_geometry false or isotropic hardening)"
+                .into(),
+        ));
+    }
+    let b = mul(&mul(f, &plastic_stretch(committed)), &transpose(f));
+    let (l, q) = symmetric_eigen(&b)
+        .filter(|(l, _)| l.iter().all(|x| *x > 0.0))
+        .ok_or_else(|| Error::Degenerate("the elastic stretch of a plastic material point is not positive (the material is inverted)".into()))?;
+    let eps = l.map(|x| 0.5 * x.ln());
+    let mut state = start(j2, committed, temperature);
+    let mut trial = trial_deviator(j2, &[eps[0], eps[1], eps[2], 0.0, 0.0, 0.0]);
+    let ret: Option<Return> = j2.radial_return(&mut state, &mut trial, 0.0);
+    // Back to the global frame: τ = Σ τ_A q_A ⊗ q_A, n likewise.
+    let spectral = |v: [f64; 3]| -> Mat3 { std::array::from_fn(|i| std::array::from_fn(|j| (0..3).map(|a| v[a] * q[i][a] * q[j][a]).sum())) };
+    let tau = voigt(&spectral([trial[0], trial[1], trial[2]]));
+    let global = ret.map(|r| Return { normal: voigt(&spectral([r.normal[0], r.normal[1], r.normal[2]])), ..r });
+    let d = j2.tangent(global.as_ref(), false).ok_or_else(singular)?;
+    // θ_AB of the logarithm's derivative.
+    let theta: Mat3 =
+        std::array::from_fn(|a| std::array::from_fn(|bb| if a == bb { 1.0 } else { log_divided_difference((l[a] - l[bb]) / (l[a] + l[bb])) }));
+    let tau_m = tensor(&tau);
+    let qt = transpose(&q);
+    let mut c: Tangent = [[0.0; 6]; 6];
+    for (col, &(i, j)) in crate::constitutive::PAIRS.iter().enumerate() {
+        // The unit rate of deformation of this column (engineering shear 1: d_ij = d_ji = ½).
+        let mut dm = [[0.0; 3]; 3];
+        let half = if i == j { 1.0 } else { 0.5 };
+        dm[i][j] = half;
+        dm[j][i] = half;
+        let mut dp = mul(&mul(&qt, &dm), &q);
+        for a in 0..3 {
+            for bb in 0..3 {
+                dp[a][bb] *= theta[a][bb];
+            }
+        }
+        let de = voigt(&mul(&mul(&q, &dp), &qt));
+        let de = [de[0], de[1], de[2], 2.0 * de[3], 2.0 * de[4], 2.0 * de[5]];
+        let spin = voigt(&{
+            let (x, y) = (mul(&dm, &tau_m), mul(&tau_m, &dm));
+            std::array::from_fn(|r| std::array::from_fn(|s| x[r][s] + y[r][s]))
+        });
+        for row in 0..6 {
+            c[row][col] = (0..6).map(|k| d[row][k] * de[k]).sum::<f64>() - spin[row];
+        }
+    }
+    if let Some(r) = ret {
+        // bₑ = Σ exp(2 εₑ,A) q_A ⊗ q_A with εₑ = ε − √(3/2) Δε̄ₚ n, and Cₚ⁻¹ = F⁻¹ bₑ F⁻ᵀ.
+        let dg = (1.5f64).sqrt() * r.plastic_strain;
+        let be = spectral(std::array::from_fn(|a| (2.0 * (eps[a] - dg * r.normal[a])).exp()));
+        let (inv, _) = inverse3(f).ok_or_else(|| Error::Degenerate("the deformation gradient is singular".into()))?;
+        store_plastic_stretch(&mul(&mul(&inv, &be), &transpose(&inv)), updated);
+        if let Some(slot) = updated.get_mut(PLASTIC_STRAIN) {
+            *slot = state.plastic_strain;
+        }
+    }
+    Ok(Response { stress: tau, tangent: c })
+}
+
+/// U′(J) and U″(J) of the Hencky volumetric energy K/2 (ln J)² (finite strain), or of
+/// K/2 (J − 1)² for the dilatation J = 1 + tr ε of small strain.
+pub fn volumetric(j2: &J2, j: f64, small: bool) -> Option<(f64, f64)> {
+    if !(j > 0.0 && j.is_finite()) {
+        return None;
+    }
+    let k = j2.bulk;
+    Some(if small { (k * (j - 1.0), k) } else { (k * j.ln() / j, k * (1.0 - j.ln()) / (j * j)) })
+}
+
+/// Small strain: the stress and the consistent tangent at total (mechanical) strain `strain`
+/// (engineering shears), from the history `committed` (written, updated, to `updated`); with
+/// `bulk` false only the deviatoric part (what a mixed element takes).
+pub fn small(j2: &J2, strain: &[f64; 6], temperature: Option<f64>, committed: &[f64], updated: &mut [f64], bulk: bool) -> Result<Response> {
+    let elastic: [f64; 6] = std::array::from_fn(|k| strain[k] - value(committed, k));
+    let mut state = start(j2, committed, temperature);
+    let mut stress = trial_deviator(j2, &elastic);
+    let ret = j2.radial_return(&mut state, &mut stress, 0.0);
+    if bulk {
+        let p = j2.bulk * (elastic[0] + elastic[1] + elastic[2]);
+        for s in stress.iter_mut().take(3) {
+            *s += p;
+        }
+    }
+    let tangent = j2.tangent(ret.as_ref(), bulk).ok_or_else(singular)?;
+    if let Some(r) = ret {
+        let dg = (1.5f64).sqrt() * r.plastic_strain;
+        for (k, slot) in updated.iter_mut().take(6).enumerate() {
+            *slot = value(committed, k) + dg * r.normal[k] * if k < 3 { 1.0 } else { 2.0 };
+        }
+        for (k, slot) in updated.iter_mut().skip(6).take(6).enumerate() {
+            *slot = state.back[k];
+        }
+        if let Some(slot) = updated.get_mut(PLASTIC_STRAIN) {
+            *slot = state.plastic_strain;
+        }
+    }
+    Ok(Response { stress, tangent })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::constitutive::{add_pressure, push_forward};
+    use fieldcraft_model::{Elastic, Hardening, HardeningPoint, Plastic};
+
+    const E: f64 = 200_000.0;
+    const NU: f64 = 0.3;
+
+    fn j2(plastic: Plastic) -> J2 {
+        J2::from_model(&Elastic::Isotropic { young: E, poisson: NU }, Some(&plastic), None).unwrap()
+    }
+
+    fn bilinear(h: f64, hardening: Hardening) -> J2 {
+        j2(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: h, hardening })
+    }
+
+    /// The whole finite-strain response, isochoric plus the Hencky volume.
+    fn finite(m: &J2, f: &Mat3, committed: &[f64], updated: &mut [f64]) -> Response {
+        let iso = finite_isochoric(m, f, None, committed, updated).unwrap();
+        let j = crate::dense::det3(f);
+        let (p, dp) = volumetric(m, j, false).unwrap();
+        add_pressure(&iso, j, p, j * dp)
+    }
+
+    /// The pulled-back tangent C (dS = C : dE) against central differences of S along several
+    /// directions of F, as `constitutive::tests::check_tangent` does, from fixed history.
+    fn check_finite_tangent(m: &J2, f: &Mat3, committed: &[f64]) -> f64 {
+        let mut scratch = vec![0.0; HISTORY];
+        let r = finite(m, f, committed, &mut scratch);
+        let (inv, _) = inverse3(f).unwrap();
+        let t = push_forward(&inv);
+        let back: Tangent =
+            std::array::from_fn(|a| std::array::from_fn(|b| (0..6).map(|k| (0..6).map(|l| t[a][k] * r.tangent[k][l] * t[b][l]).sum::<f64>()).sum()));
+        let scale = back.iter().flatten().fold(0.0f64, |a, v| a.max(v.abs()));
+        let h = 1e-7;
+        let mut worst: f64 = 0.0;
+        for k in 0..9 {
+            let mut df = [[0.0; 3]; 3];
+            df[k / 3][k % 3] = 1.0;
+            df[(k + 1) % 3][(k + 2) % 3] += 0.3;
+            let at = |s: f64| {
+                let fp: Mat3 = std::array::from_fn(|i| std::array::from_fn(|j| f[i][j] + s * df[i][j]));
+                let mut scratch = vec![0.0; HISTORY];
+                let tau = finite(m, &fp, committed, &mut scratch).stress;
+                crate::rotation::rotate(&inverse3(&fp).unwrap().0, &tau)
+            };
+            let (sp, sm) = (at(h), at(-h));
+            let ftdf = mul(&transpose(f), &df);
+            let de = [ftdf[0][0], ftdf[1][1], ftdf[2][2], ftdf[0][1] + ftdf[1][0], ftdf[0][2] + ftdf[2][0], ftdf[1][2] + ftdf[2][1]];
+            for a in 0..6 {
+                let numeric = (sp[a] - sm[a]) / (2.0 * h);
+                let analytic: f64 = (0..6).map(|b| back[a][b] * de[b]).sum();
+                worst = worst.max((numeric - analytic).abs() / scale);
+            }
+        }
+        worst
+    }
+
+    #[test]
+    fn finite_strain_tangent_is_consistent_elastic_and_plastic() {
+        let deformations: [Mat3; 3] = [
+            [[1.02, 0.015, -0.005], [-0.01, 0.99, 0.02], [0.007, -0.012, 1.005]],
+            // Large: well into plasticity.
+            [[1.3, 0.15, -0.05], [-0.1, 0.85, 0.2], [0.07, -0.12, 1.05]],
+            // Uniaxial (two equal principal stretches).
+            [[1.2, 0.0, 0.0], [0.0, 0.92, 0.0], [0.0, 0.0, 0.92]],
+        ];
+        let tables = Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: 250.0 }, HardeningPoint { strain: 0.5, stress: 600.0 }] };
+        let jc = 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: None,
+            taylor_quinney: 0.9,
+        };
+        // A plastic history to start from: a prior stretch committed.
+        let mut prior = vec![0.0; HISTORY];
+        finite(&bilinear(2000.0, Hardening::Isotropic), &[[1.1, 0.02, 0.0], [0.0, 0.95, 0.03], [0.0, 0.0, 0.96]], &[0.0; HISTORY], &mut prior);
+        assert!(prior[PLASTIC_STRAIN] > 0.0);
+        for m in [bilinear(2000.0, Hardening::Isotropic), bilinear(0.0, Hardening::Isotropic), j2(tables), j2(jc)] {
+            for f in &deformations {
+                for committed in [&[0.0; HISTORY][..], &prior[..]] {
+                    let worst = check_finite_tangent(&m, f, committed);
+                    assert!(worst < 1e-7, "{:?} at {f:?}: {worst}", m.flow);
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn finite_strain_tangent_is_symmetric_and_elastic_at_rest() {
+        let m = bilinear(2000.0, Hardening::Isotropic);
+        let mut scratch = vec![0.0; HISTORY];
+        let r = finite(&m, &[[1.3, 0.15, -0.05], [-0.1, 0.85, 0.2], [0.07, -0.12, 1.05]], &[0.0; HISTORY], &mut scratch);
+        for a in 0..6 {
+            for b in 0..6 {
+                assert!((r.tangent[a][b] - r.tangent[b][a]).abs() < 1e-9 * E, "{a} {b}");
+            }
+        }
+        let rest = finite(&m, &IDENTITY, &[0.0; HISTORY], &mut scratch);
+        let d = crate::material::isotropic(E, NU).unwrap();
+        for a in 0..6 {
+            for b in 0..6 {
+                assert!((rest.tangent[a][b] - d[a][b]).abs() < 1e-9 * E);
+            }
+        }
+        assert!(rest.stress.iter().all(|v| *v == 0.0));
+    }
+
+    /// Uniaxial stress along x at finite strain, from a virgin point to a logarithmic strain ε
+    /// in one increment (the lateral stretch found so σ_yy = 0): the Kirchhoff stress is the
+    /// small-strain curve in the logarithmic strain, τ = σ_y + E_t (ε − σ_y/E), because the
+    /// exponential map integrates a proportional path exactly.
+    #[test]
+    fn finite_uniaxial_matches_the_logarithmic_closed_form() {
+        let (y0, h) = (250.0, 2000.0);
+        let m = bilinear(h, Hardening::Isotropic);
+        let strain: f64 = 0.2;
+        let lam = strain.exp();
+        let sigma = |l: f64| {
+            let f = [[lam, 0.0, 0.0], [0.0, l, 0.0], [0.0, 0.0, l]];
+            let mut u = vec![0.0; HISTORY];
+            let r = finite(&m, &f, &[0.0; HISTORY], &mut u);
+            let j = lam * l * l;
+            (r.stress[0] / j, r.stress[1] / j, u)
+        };
+        // Secant on the lateral stretch.
+        let (mut a, mut b) = (1.0, 0.9);
+        for _ in 0..60 {
+            let (fa, fb) = (sigma(a).1, sigma(b).1);
+            if fb.abs() < 1e-11 || fb == fa {
+                break;
+            }
+            (a, b) = (b, b - fb * (b - a) / (fb - fa));
+        }
+        let (sxx, _, u) = sigma(b);
+        let j = lam * b * b;
+        // In Kirchhoff stress τ = J σ: τ = y0 + E_t (ε − y0/E) along the log strain path.
+        let et = E * h / (E + h);
+        let exact = y0 + et * (strain - y0 / E);
+        assert!((sxx * j - exact).abs() < 1e-8 * exact, "{} vs {exact}", sxx * j);
+        assert!((u[PLASTIC_STRAIN] - (exact - y0) / h).abs() < 1e-10);
+        // The plastic deformation is isochoric: det Cₚ⁻¹ = 1.
+        let cp = plastic_stretch(&u);
+        assert!((crate::dense::det3(&cp) - 1.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn small_strain_tangent_is_consistent_and_cyclic_kinematic_shows_bauschinger() {
+        let m = bilinear(10_000.0, Hardening::Kinematic);
+        let strain = [3e-3, -1e-3, 5e-4, 2e-3, -1e-3, 7e-4];
+        let mut prior = vec![0.0; HISTORY];
+        small(&m, &[2e-3, -1e-3, -1e-3, 0.0, 0.0, 0.0], None, &[0.0; HISTORY], &mut prior, true).unwrap();
+        for committed in [&[0.0; HISTORY][..], &prior[..]] {
+            let mut scratch = vec![0.0; HISTORY];
+            let r = small(&m, &strain, None, committed, &mut scratch, true).unwrap();
+            let h = 1e-9;
+            for b in 0..6 {
+                let at = |s: f64| {
+                    let mut e = strain;
+                    e[b] += s;
+                    small(&m, &e, None, committed, &mut [0.0; HISTORY], true).unwrap().stress
+                };
+                let (p, q) = (at(h), at(-h));
+                for a in 0..6 {
+                    let numeric = (p[a] - q[a]) / (2.0 * h);
+                    assert!((numeric - r.tangent[a][b]).abs() < 1e-5 * E, "{a} {b}: {numeric} vs {}", r.tangent[a][b]);
+                }
+            }
+        }
+        // Kinematic hardening is refused at finite strain, named.
+        assert!(matches!(finite_isochoric(&m, &IDENTITY, None, &[0.0; HISTORY], &mut [0.0; HISTORY]), Err(Error::NotSupported(_))));
+    }
+
+    #[test]
+    fn softening_steeper_than_3g_has_no_tangent() {
+        let g = E / (2.0 * (1.0 + NU));
+        let points = vec![HardeningPoint { strain: 0.0, stress: 250.0 }, HardeningPoint { strain: 1e-4, stress: 250.0 - 4.0 * g * 1e-4 }];
+        let m = j2(Plastic::Table { points });
+        let r = small(&m, &[5e-3, -2.5e-3, -2.5e-3, 0.0, 0.0, 0.0], None, &[0.0; HISTORY], &mut [0.0; HISTORY], true);
+        assert!(matches!(r, Err(Error::Degenerate(_))), "{r:?}");
+    }
+}
diff --git a/crates/fem/src/finite.rs b/crates/fem/src/finite.rs
index fad60c1..2ba5d08 100644
--- a/crates/fem/src/finite.rs
+++ b/crates/fem/src/finite.rs
@@ -1,6 +1,6 @@
-//! Finite-strain (total Lagrangian) solid and 2D elements: internal forces and the consistent
-//! tangent stiffness for nodal displacements u, and the Cauchy stress and logarithmic strain at
-//! the integration points.
+//! Solid and 2D elements of nonlinear steps: total Lagrangian (finite strain) or with small
+//! strain and an inelastic material; internal forces and the consistent tangent stiffness for
+//! nodal displacements u, and the Cauchy stress and the strain at the integration points.
 //!
 //! **Formulation.** Total Lagrangian: the integrals run over the reference (undeformed) element
 //! with its own quadrature, F = I + Σ u_a ⊗ ∇₀N_a at each point, and the material returns the
@@ -26,7 +26,15 @@
 //! symmetric; at u = 0 it is the linear hybrid stiffness K_dev + κ G M⁻¹ Gᵀ. The constant
 //! pressure case is the mean dilatation method of J. C. Nagtegaal, D. M. Parks, J. R. Rice
 //! (CMAME 4:153–177, 1974) at finite strain. St Venant–Kirchhoff has no such split, so hybrid
-//! elements of an elastic material are refused in nonlinear steps.
+//! elements of an elastic material are refused in nonlinear steps; J2 plasticity has it exactly
+//! (the return is deviatoric, the volume Hencky's K/2 (ln J)², [`crate::elastoplastic`]).
+//!
+//! **Small strain** ([`Kinematics::Small`], a step without nonlinear geometry with an inelastic
+//! material): the same loops on the linear kinematics, ε = sym ∇u with the reference gradients,
+//! f = ∫ Bᵀ σ dV₀ and K = ∫ Bᵀ 𝔻 B dV₀ with the material's consistent tangent 𝔻, and no
+//! geometric stiffness. Mixed elements project the dilatation 1 + tr ε and the pressure K tr ε
+//! the same way, which is the linear hybrid element of [`crate::hybrid`] while the material is
+//! elastic.
 //!
 //! **2D elements**: plane strain (F₃₃ = 1) and axisymmetric (F₃₃ = r/R, the hoop stretch, with
 //! the hoop column N/r of the current radius in Bₓ and a geometric term τ_θθ N_a N_b / r² on the
@@ -37,7 +45,8 @@
 //! (enhanced assumed strain, J. C. Simo, F. Armero, IJNME 33:1413–1449, 1992) needs the internal
 //! modes as element unknowns carried between iterations, and it develops spurious hourglass
 //! modes in large compression (P. Wriggers, S. Reese, CMAME 135:201–209, 1996), which want a
-//! stabilised variant of their own.
+//! stabilised variant of their own. With small strain they are refused for inelastic materials
+//! too: their internal modes would carry the plastic history between iterations.
 
 use crate::constitutive::{MaterialModel, Response, add_pressure, logarithmic_strain};
 use crate::dense::{Cholesky, Matrix, inverse3};
@@ -54,12 +63,25 @@ pub struct Finite {
     pub tangent: Option<Matrix>,
 }
 
+/// The kinematics of an evaluation.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Kinematics {
+    /// Total Lagrangian, large displacements, rotations and strains.
+    Finite,
+    /// Small displacements and strains (ε = sym ∇u, the reference gradients, no geometric
+    /// stiffness): only the material is nonlinear.
+    Small,
+}
+
 /// What one integration point gives the element.
 struct PointKinematics {
     weight: f64,
     f: Mat3,
+    /// det F, or the dilatation 1 + tr ε of small strain (what a mixed element projects).
     j: f64,
-    /// ∇ₓN of every node.
+    /// The small strain ε = sym ∇u (engineering shears); zero at finite strain.
+    strain: [f64; 6],
+    /// ∇ₓN of every node (∇₀N at small strain).
     gradients: Vec<[f64; 3]>,
     /// Reference position.
     x0: [f64; 3],
@@ -70,8 +92,22 @@ fn point_history(history: &[f64], len: usize, point: usize) -> &[f64] {
     history.get(point * len..(point + 1) * len).unwrap_or_default()
 }
 
-/// Add f_a += w Bₓ,aᵀ τ and, with `k`, the material and geometric tangent.
-fn add_point(force: &mut [f64], k: Option<&mut Matrix>, gradients: &[[f64; 3]], r: &Response, w: f64) {
+/// The small strain of a displacement gradient F − I (engineering shears).
+fn small_strain(f: &Mat3) -> [f64; 6] {
+    [f[0][0] - 1.0, f[1][1] - 1.0, f[2][2] - 1.0, f[0][1] + f[1][0], f[0][2] + f[2][0], f[1][2] + f[2][1]]
+}
+
+/// The response of one point's material with its history slices.
+fn point_response(m: &MaterialModel, kin: Kinematics, p: &PointKinematics, committed: &[f64], updated: &mut [f64]) -> Result<Response> {
+    match kin {
+        Kinematics::Finite => m.response(&p.f, None, committed, updated),
+        Kinematics::Small => m.small_strain(&p.strain, None, committed, updated, true),
+    }
+}
+
+/// Add f_a += w Bₓ,aᵀ τ and, with `k`, the material tangent and (with `geometric`) the
+/// geometric one.
+fn add_point(force: &mut [f64], k: Option<&mut Matrix>, gradients: &[[f64; 3]], r: &Response, w: f64, geometric: bool) {
     for (a, g) in gradients.iter().enumerate() {
         let b = b_block(g);
         for c in 0..3 {
@@ -82,6 +118,9 @@ fn add_point(force: &mut [f64], k: Option<&mut Matrix>, gradients: &[[f64; 3]],
     }
     let Some(k) = k else { return };
     accumulate_symmetric(&mut k.data, gradients, &r.tangent, w);
+    if !geometric {
+        return;
+    }
     let tau = crate::rotation::tensor(&r.stress);
     let n = k.n;
     for (a, ga) in gradients.iter().enumerate() {
@@ -122,14 +161,17 @@ fn check_history(m: &MaterialModel, points: usize, committed: &[f64], updated: &
 }
 
 impl Solid<'_> {
-    /// Refuse what has no finite-strain form here.
-    fn check_finite(&self, m: &MaterialModel, u: &[[f64; 3]]) -> Result<()> {
+    /// Refuse what has no form here.
+    fn check_nonlinear(&self, m: &MaterialModel, u: &[[f64; 3]], kin: Kinematics) -> Result<()> {
         check_displacements(u, self.ty.nodes)?;
         match self.formulation {
-            Formulation::IncompatibleModes => Err(Error::NotSupported(format!(
-                "{} incompatible modes in nonlinear steps (their finite-strain form, enhanced assumed strain, is not implemented; use full integration, hybrid or hex20)",
-                self.ty.kind.name()
-            ))),
+            Formulation::IncompatibleModes => {
+                let why = match kin {
+                    Kinematics::Finite => "in nonlinear steps (their finite-strain form, enhanced assumed strain, is not implemented",
+                    Kinematics::Small => "with an inelastic material (their internal modes would carry the plastic history",
+                };
+                Err(Error::NotSupported(format!("{} incompatible modes {why}; use full integration, hybrid or hex20)", self.ty.kind.name())))
+            }
             Formulation::Hybrid { .. } if !m.is_decoupled() => Err(Error::NotSupported(format!(
                 "hybrid elements of a {} material in nonlinear steps (the mixed formulation needs a volumetric–isochoric split of the strain energy; give the material the neo-Hookean option)",
                 m.name()
@@ -138,59 +180,78 @@ impl Solid<'_> {
         }
     }
 
-    /// F, J and the current gradients at every stiffness point.
-    fn kinematics(&self, u: &[[f64; 3]]) -> Result<Vec<PointKinematics>> {
+    /// F, J and the gradients at every stiffness point: the current ones at finite strain, the
+    /// reference ones at small strain.
+    fn kinematics(&self, u: &[[f64; 3]], kin: Kinematics) -> Result<Vec<PointKinematics>> {
         let n = self.ty.nodes;
         let mut out = Vec::new();
         for p in self.rule().points() {
             let g = self.gradients(p.xi)?;
             let reference = g.dx.get(..n).unwrap_or_default();
             let f = crate::constitutive::deformation_gradient(u, reference);
+            let (weight, x0) = (g.det * p.weight, self.position_of(g.shape.values()));
+            if kin == Kinematics::Small {
+                let strain = small_strain(&f);
+                let j = 1.0 + strain[0] + strain[1] + strain[2];
+                out.push(PointKinematics { weight, f, j, strain, gradients: reference.to_vec(), x0 });
+                continue;
+            }
             let (inv, j) = inverse3(&f).ok_or_else(|| inverted(self.ty.kind.name(), 0.0))?;
             if !(j > 0.0 && j.is_finite()) {
                 return Err(inverted(self.ty.kind.name(), j));
             }
             let gradients = reference.iter().map(|d| std::array::from_fn(|i| (0..3).map(|jj| d[jj] * inv[jj][i]).sum())).collect();
-            out.push(PointKinematics { weight: g.det * p.weight, f, j, gradients, x0: self.position_of(g.shape.values()) });
+            out.push(PointKinematics { weight, f, j, strain: [0.0; 6], gradients, x0 });
         }
         Ok(out)
     }
 
-    /// The effective response (Kirchhoff stress and tangent) at every point, the forces, and the
-    /// tangent when `tangent`.
-    fn finite_points(
+    /// The effective response (Kirchhoff stress and tangent, or the stress of small strain) at
+    /// every point, the forces, and the tangent when `tangent`.
+    fn nonlinear_points(
         &self,
         m: &MaterialModel,
         u: &[[f64; 3]],
+        kin: Kinematics,
         committed: &[f64],
         updated: &mut [f64],
         tangent: bool,
     ) -> Result<(Finite, Vec<(PointKinematics, Response)>)> {
-        self.check_finite(m, u)?;
-        let points = self.kinematics(u)?;
+        self.check_nonlinear(m, u, kin)?;
+        let points = self.kinematics(u, kin)?;
         let len = check_history(m, points.len(), committed, updated)?;
         let ndof = 3 * self.ty.nodes;
         let mut force = vec![0.0; ndof];
         let mut k = tangent.then(|| Matrix::zeros(ndof));
         let responses: Vec<Response> = if let Formulation::Hybrid { .. } = self.formulation {
-            self.mixed_responses(m, &points, k.as_mut())?
+            self.mixed_responses(m, &points, kin, committed, updated, k.as_mut())?
         } else {
             let mut out = Vec::with_capacity(points.len());
             for (i, p) in points.iter().enumerate() {
                 let updated_point = updated.get_mut(i * len..(i + 1) * len).unwrap_or_default();
-                out.push(m.response(&p.f, point_history(committed, len, i), updated_point)?);
+                out.push(point_response(m, kin, p, point_history(committed, len, i), updated_point)?);
             }
             out
         };
         for (p, r) in points.iter().zip(&responses) {
-            add_point(&mut force, k.as_mut(), &p.gradients, r, p.weight);
+            add_point(&mut force, k.as_mut(), &p.gradients, r, p.weight, kin == Kinematics::Finite);
         }
         Ok((Finite { force, tangent: k }, points.into_iter().zip(responses).collect()))
     }
 
-    /// The three-field mixed responses τ_iso + J p I with their tangents, and the dilatation
-    /// stiffness Gᵀ M⁻¹ H M⁻¹ G added to `k`.
-    fn mixed_responses(&self, m: &MaterialModel, points: &[PointKinematics], k: Option<&mut Matrix>) -> Result<Vec<Response>> {
+    /// The three-field mixed responses τ_iso + J p I with their tangents (small strain:
+    /// s + p I), and the dilatation stiffness Gᵀ M⁻¹ H M⁻¹ G added to `k`.
+    fn mixed_responses(
+        &self,
+        m: &MaterialModel,
+        points: &[PointKinematics],
+        kin: Kinematics,
+        committed: &[f64],
+        updated: &mut [f64],
+        k: Option<&mut Matrix>,
+    ) -> Result<Vec<Response>> {
+        let small = kin == Kinematics::Small;
+        let len = m.history_len();
         let (np, frame) = (self.np(), self.frame());
         let basis: Vec<[f64; 4]> = points.iter().map(|p| pressure_basis(np, p.x0, &frame)).collect();
         let mut mm = Matrix::zeros(np);
@@ -212,7 +273,7 @@ impl Solid<'_> {
             if !(th > 0.0 && th.is_finite()) {
                 return Err(inverted(self.ty.kind.name(), th));
             }
-            let (dp, ddp) = m.volumetric(th).ok_or_else(|| Error::Invalid("volumetric response".into()))?;
+            let (dp, ddp) = m.volumetric(th, small).ok_or_else(|| Error::Invalid("volumetric response".into()))?;
             for a in 0..np {
                 pressure[a] += p.weight * nb[a] * dp;
                 for b in 0..np {
@@ -222,20 +283,28 @@ impl Solid<'_> {
         }
         chol.solve(&mut pressure);
         let mut out = Vec::with_capacity(points.len());
-        for (p, nb) in points.iter().zip(&basis) {
+        for (i, (p, nb)) in points.iter().zip(&basis).enumerate() {
             let pr: f64 = (0..np).map(|a| nb[a] * pressure[a]).sum();
-            out.push(add_pressure(&m.isochoric(&p.f)?, p.j, pr, 0.0));
+            let (c, u) = (point_history(committed, len, i), updated.get_mut(i * len..(i + 1) * len).unwrap_or_default());
+            out.push(if small {
+                let mut r = m.small_strain(&p.strain, None, c, u, false)?;
+                r.stress.iter_mut().take(3).for_each(|s| *s += pr);
+                r
+            } else {
+                add_pressure(&m.isochoric(&p.f, None, c, u)?, p.j, pr, 0.0)
+            });
         }
         if let Some(k) = k {
-            // A = M⁻¹ G, np × ndof; K += Aᵀ H A.
+            // A = M⁻¹ G, np × ndof; K += Aᵀ H A. At small strain G = ∫ N_p ∇₀N dV₀.
             let ndof = k.n;
             let mut a_cols = vec![vec![0.0; np]; ndof];
             for (p, nb) in points.iter().zip(&basis) {
+                let jw = p.weight * if small { 1.0 } else { p.j };
                 for (node, g) in p.gradients.iter().enumerate() {
                     for c in 0..3 {
                         if let Some(col) = a_cols.get_mut(3 * node + c) {
                             for (x, n) in col.iter_mut().zip(nb) {
-                                *x += p.weight * n * p.j * g[c];
+                                *x += jw * n * g[c];
                             }
                         }
                     }
@@ -262,17 +331,35 @@ impl Solid<'_> {
     /// with the material history `committed` at the start of the increment (points × the
     /// material's history length) written, updated, to `updated`.
     pub fn finite_strain(&self, m: &MaterialModel, u: &[[f64; 3]], committed: &[f64], updated: &mut [f64], tangent: bool) -> Result<Finite> {
-        Ok(self.finite_points(m, u, committed, updated, tangent)?.0)
+        Ok(self.nonlinear_points(m, u, Kinematics::Finite, committed, updated, tangent)?.0)
     }
 
-    /// The Cauchy stress σ = τ/J and the logarithmic strain at every integration point (at its
-    /// reference position).
-    pub fn finite_strain_stresses(&self, m: &MaterialModel, u: &[[f64; 3]], committed: &[f64]) -> Result<Vec<PointStress>> {
+    /// As [`Self::finite_strain`] with either kinematics; with small strain an inelastic
+    /// material on the linear kinematics, f = ∫ Bᵀ σ dV and K = ∫ Bᵀ 𝔻 B dV with the
+    /// material's consistent tangent 𝔻.
+    pub fn nonlinear(
+        &self,
+        m: &MaterialModel,
+        u: &[[f64; 3]],
+        kin: Kinematics,
+        committed: &[f64],
+        updated: &mut [f64],
+        tangent: bool,
+    ) -> Result<Finite> {
+        Ok(self.nonlinear_points(m, u, kin, committed, updated, tangent)?.0)
+    }
+
+    /// The Cauchy stress and the strain (logarithmic at finite strain, ε at small strain) at
+    /// every integration point (at its reference position).
+    pub fn nonlinear_stresses(&self, m: &MaterialModel, u: &[[f64; 3]], kin: Kinematics, committed: &[f64]) -> Result<Vec<PointStress>> {
         let mut scratch = committed.to_vec();
-        let (_, points) = self.finite_points(m, u, committed, &mut scratch, false)?;
+        let (_, points) = self.nonlinear_points(m, u, kin, committed, &mut scratch, false)?;
         points
             .iter()
             .map(|(p, r)| {
+                if kin == Kinematics::Small {
+                    return Ok(PointStress { xyz: p.x0, strain: p.strain, stress: r.stress });
+                }
                 let strain = logarithmic_strain(&p.f).ok_or_else(|| inverted(self.ty.kind.name(), p.j))?;
                 Ok(PointStress { xyz: p.x0, strain, stress: r.stress.map(|v| v / p.j) })
             })
@@ -280,14 +367,26 @@ impl Solid<'_> {
     }
 }
 
+/// One integration point of a 2D element in a nonlinear step.
+struct PlanePoint {
+    f: Mat3,
+    j: f64,
+    x: [f64; 3],
+    strain: [f64; 6],
+    response: Response,
+}
+
 impl Plane<'_> {
-    fn check_finite(&self, u: &[[f64; 2]]) -> Result<()> {
+    fn check_nonlinear(&self, u: &[[f64; 2]], kin: Kinematics) -> Result<()> {
         check_displacements(u, self.ty.nodes)?;
         if let Theory::PlaneStress { .. } = self.theory {
-            return Err(Error::NotSupported(
-                "plane stress elements in nonlinear steps (the thickness change needs a local iteration; use plane strain, axisymmetric or solid elements)"
-                    .into(),
-            ));
+            let why = match kin {
+                Kinematics::Finite => "the thickness change needs a local iteration",
+                Kinematics::Small => "σ_zz = 0 needs a local iteration of the inelastic material",
+            };
+            return Err(Error::NotSupported(format!(
+                "plane stress elements in nonlinear steps ({why}; use plane strain, axisymmetric or solid elements)"
+            )));
         }
         if self.formulation == Formulation::IncompatibleModes {
             return Err(Error::NotSupported(format!(
@@ -298,17 +397,18 @@ impl Plane<'_> {
         Ok(())
     }
 
-    /// Per point: (dv, F, J, current columns (∂N/∂x, ∂N/∂y, N/r), reference position, response).
-    #[allow(clippy::type_complexity)]
-    fn finite_points(
+    /// Per point: F, J, the reference position, the small strain and the response; the forces
+    /// and the tangent.
+    fn nonlinear_points(
         &self,
         m: &MaterialModel,
         u: &[[f64; 2]],
+        kin: Kinematics,
         committed: &[f64],
         updated: &mut [f64],
         tangent: bool,
-    ) -> Result<(Finite, Vec<(Mat3, f64, [f64; 3], Response)>)> {
-        self.check_finite(u)?;
+    ) -> Result<(Finite, Vec<PlanePoint>)> {
+        self.check_nonlinear(u, kin)?;
         let n = self.ty.nodes;
         let rule = self.rule().points();
         let len = check_history(m, rule.len(), committed, updated)?;
@@ -316,6 +416,7 @@ impl Plane<'_> {
         let mut force = vec![0.0; ndof];
         let mut k = tangent.then(|| Matrix::zeros(ndof));
         let axisymmetric = self.theory == Theory::Axisymmetric;
+        let small = kin == Kinematics::Small;
         let mut out = Vec::with_capacity(rule.len());
         for (i, p) in rule.iter().enumerate() {
             let at = self.at(p)?;
@@ -331,16 +432,23 @@ impl Plane<'_> {
                     f[2][2] += ua[0] * g[2];
                 }
             }
-            let det2 = f[0][0] * f[1][1] - f[0][1] * f[1][0];
-            let j = det2 * f[2][2];
-            if !(j > 0.0 && j.is_finite() && f[2][2] > 0.0) {
-                return Err(inverted(&format!("2D {}", self.ty.kind.name()), j));
-            }
-            let inv = [[f[1][1] / det2, -f[0][1] / det2], [-f[1][0] / det2, f[0][0] / det2]];
-            let current: Vec<[f64; 3]> =
-                cols.iter().map(|g| [g[0] * inv[0][0] + g[1] * inv[1][0], g[0] * inv[0][1] + g[1] * inv[1][1], g[2] / f[2][2]]).collect();
             let updated_point = updated.get_mut(i * len..(i + 1) * len).unwrap_or_default();
-            let r = m.response(&f, point_history(committed, len, i), updated_point)?;
+            let (j, strain, current, r) = if small {
+                let strain = small_strain(&f);
+                let r = m.small_strain(&strain, None, point_history(committed, len, i), updated_point, true)?;
+                (1.0 + strain[0] + strain[1] + strain[2], strain, cols.to_vec(), r)
+            } else {
+                let det2 = f[0][0] * f[1][1] - f[0][1] * f[1][0];
+                let j = det2 * f[2][2];
+                if !(j > 0.0 && j.is_finite() && f[2][2] > 0.0) {
+                    return Err(inverted(&format!("2D {}", self.ty.kind.name()), j));
+                }
+                let inv = [[f[1][1] / det2, -f[0][1] / det2], [-f[1][0] / det2, f[0][0] / det2]];
+                let current: Vec<[f64; 3]> =
+                    cols.iter().map(|g| [g[0] * inv[0][0] + g[1] * inv[1][0], g[0] * inv[0][1] + g[1] * inv[1][1], g[2] / f[2][2]]).collect();
+                let r = m.response(&f, None, point_history(committed, len, i), updated_point)?;
+                (j, [0.0; 6], current, r)
+            };
             let tau = [r.stress[0], r.stress[1], r.stress[2], r.stress[3]];
             for (a, g) in current.iter().enumerate() {
                 if let Some(x) = force.get_mut(2 * a) {
@@ -353,15 +461,17 @@ impl Plane<'_> {
             if let Some(k) = k.as_mut() {
                 let c4: [[f64; 4]; 4] = std::array::from_fn(|a| std::array::from_fn(|b| r.tangent[a][b]));
                 crate::plane::accumulate(&mut k.data, &current, &c4, at.dv);
-                for (a, ga) in current.iter().enumerate() {
-                    for (b, gb) in current.iter().enumerate() {
-                        let s = ga[0] * gb[0] * tau[0] + ga[1] * gb[1] * tau[1] + (ga[0] * gb[1] + ga[1] * gb[0]) * tau[3];
-                        k.add(2 * a, 2 * b, at.dv * (s + tau[2] * ga[2] * gb[2]));
-                        k.add(2 * a + 1, 2 * b + 1, at.dv * s);
+                if !small {
+                    for (a, ga) in current.iter().enumerate() {
+                        for (b, gb) in current.iter().enumerate() {
+                            let s = ga[0] * gb[0] * tau[0] + ga[1] * gb[1] * tau[1] + (ga[0] * gb[1] + ga[1] * gb[0]) * tau[3];
+                            k.add(2 * a, 2 * b, at.dv * (s + tau[2] * ga[2] * gb[2]));
+                            k.add(2 * a + 1, 2 * b + 1, at.dv * s);
+                        }
                     }
                 }
             }
-            out.push((f, j, at.x, r));
+            out.push(PlanePoint { f, j, x: at.x, strain, response: r });
         }
         Ok((Finite { force, tangent: k }, out))
     }
@@ -369,19 +479,35 @@ impl Plane<'_> {
     /// Internal forces (u, v per node) and (with `tangent`) the tangent stiffness at nodal
     /// displacements `u`, as for [`Solid::finite_strain`].
     pub fn finite_strain(&self, m: &MaterialModel, u: &[[f64; 2]], committed: &[f64], updated: &mut [f64], tangent: bool) -> Result<Finite> {
-        Ok(self.finite_points(m, u, committed, updated, tangent)?.0)
+        Ok(self.nonlinear_points(m, u, Kinematics::Finite, committed, updated, tangent)?.0)
     }
 
-    /// The Cauchy stress and logarithmic strain at every integration point, as the solid
-    /// elements' 6-vectors (axisymmetric: rr, axial, hoop, r-axial).
-    pub fn finite_strain_stresses(&self, m: &MaterialModel, u: &[[f64; 2]], committed: &[f64]) -> Result<Vec<PointStress>> {
+    /// As [`Solid::nonlinear`], on the 2D element.
+    pub fn nonlinear(
+        &self,
+        m: &MaterialModel,
+        u: &[[f64; 2]],
+        kin: Kinematics,
+        committed: &[f64],
+        updated: &mut [f64],
+        tangent: bool,
+    ) -> Result<Finite> {
+        Ok(self.nonlinear_points(m, u, kin, committed, updated, tangent)?.0)
+    }
+
+    /// The Cauchy stress and the strain at every integration point, as the solid elements'
+    /// 6-vectors (axisymmetric: rr, axial, hoop, r-axial).
+    pub fn nonlinear_stresses(&self, m: &MaterialModel, u: &[[f64; 2]], kin: Kinematics, committed: &[f64]) -> Result<Vec<PointStress>> {
         let mut scratch = committed.to_vec();
-        let (_, points) = self.finite_points(m, u, committed, &mut scratch, false)?;
+        let (_, points) = self.nonlinear_points(m, u, kin, committed, &mut scratch, false)?;
         points
             .iter()
-            .map(|(f, j, x, r)| {
-                let strain = logarithmic_strain(f).ok_or_else(|| inverted(&format!("2D {}", self.ty.kind.name()), *j))?;
-                Ok(PointStress { xyz: *x, strain, stress: r.stress.map(|v| v / j) })
+            .map(|p| {
+                if kin == Kinematics::Small {
+                    return Ok(PointStress { xyz: p.x, strain: p.strain, stress: p.response.stress });
+                }
+                let strain = logarithmic_strain(&p.f).ok_or_else(|| inverted(&format!("2D {}", self.ty.kind.name()), p.j))?;
+                Ok(PointStress { xyz: p.x, strain, stress: p.response.stress.map(|v| v / p.j) })
             })
             .collect()
     }
@@ -420,20 +546,26 @@ mod tests {
 
     /// Central differences of the internal forces against the tangent.
     fn check_tangent(solid: &Solid, m: &MaterialModel, u: &[[f64; 3]]) {
-        let k = solid.finite_strain(m, u, &[], &mut [], true).unwrap().tangent.unwrap();
+        check_tangent_from(solid, m, u, Kinematics::Finite, &[], 1e-6);
+    }
+
+    /// The same from a committed history, with either kinematics.
+    fn check_tangent_from(solid: &Solid, m: &MaterialModel, u: &[[f64; 3]], kin: Kinematics, committed: &[f64], tol: f64) {
+        let mut scratch = committed.to_vec();
+        let k = solid.nonlinear_points(m, u, kin, committed, &mut scratch, true).unwrap().0.tangent.unwrap();
         let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs()));
         let n = 3 * u.len();
-        let h = 1e-6;
+        let h = if kin == Kinematics::Small { 1e-9 } else { 1e-6 };
         for j in 0..n {
             let shifted = |s: f64| {
                 let mut v = u.to_vec();
                 v[j / 3][j % 3] += s;
-                solid.finite_strain(m, &v, &[], &mut [], false).unwrap().force
+                solid.nonlinear_points(m, &v, kin, committed, &mut committed.to_vec(), false).unwrap().0.force
             };
             let (p, q) = (shifted(h), shifted(-h));
             for i in 0..n {
                 let numeric = (p[i] - q[i]) / (2.0 * h);
-                assert!((numeric - k.get(i, j)).abs() < 1e-6 * scale, "{:?} {i} {j}: {numeric} vs {}", solid.ty.kind, k.get(i, j));
+                assert!((numeric - k.get(i, j)).abs() < tol * scale, "{:?} {kin:?} {i} {j}: {numeric} vs {}", solid.ty.kind, k.get(i, j));
             }
         }
         for i in 0..n {
@@ -493,6 +625,108 @@ mod tests {
         }
     }
 
+    /// J2 plasticity in full and hybrid elements, at finite and small strain: the tangent is
+    /// the derivative of the forces from a plastic history too, and the plastic points are
+    /// really plastic.
+    #[test]
+    fn plastic_tangents_are_the_derivative_of_the_forces() {
+        use crate::elastoplastic::{HISTORY, PLASTIC_STRAIN};
+        use fieldcraft_model::{Elastic, Hardening, Plastic};
+        let j2 = crate::plasticity::J2::from_model(
+            &Elastic::Isotropic { young: 1000.0, poisson: 0.3 },
+            Some(&Plastic::Bilinear { yield_stress: 5.0, hardening_modulus: 50.0, hardening: Hardening::Isotropic }),
+            None,
+        )
+        .unwrap();
+        let m = MaterialModel::Plastic(std::sync::Arc::new(j2));
+        let (x20, x10) = (hex20(), tet10());
+        let x8: Vec<[f64; 3]> = x20[..8].to_vec();
+        for (kind, x, f) in [
+            (ElementKind::Hex8, &x8, Formulation::Full),
+            (ElementKind::Tet10, &x10, Formulation::Full),
+            (ElementKind::Hex20, &x20, Formulation::Reduced),
+            (ElementKind::Hex8, &x8, Formulation::Hybrid { pressures: 1 }),
+            (ElementKind::Hex20, &x20, Formulation::Hybrid { pressures: 4 }),
+        ] {
+            let solid = Solid::new(kind, x, f).unwrap();
+            let points = solid.rule().points().len();
+            for (kin, scale) in [(Kinematics::Finite, 0.05), (Kinematics::Small, 0.05)] {
+                // A first increment from the virgin state, committed, then a second from it.
+                let mut committed = vec![0.0; points * HISTORY];
+                let first = displacement(x.len(), 0.5 * scale);
+                solid.nonlinear_points(&m, &first, kin, &vec![0.0; points * HISTORY], &mut committed, false).unwrap();
+                assert!((0..points).all(|p| committed[p * HISTORY + PLASTIC_STRAIN] > 0.0), "{kind:?} {f:?} {kin:?} yields everywhere");
+                check_tangent_from(&solid, &m, &displacement(x.len(), scale), kin, &committed, 1e-6);
+            }
+        }
+        // Incompatible modes have no inelastic form.
+        let im = Solid::new(ElementKind::Hex8, &x8, Formulation::IncompatibleModes).unwrap();
+        assert!(matches!(
+            im.nonlinear(&m, &displacement(8, 0.0), Kinematics::Small, &[0.0; 8 * HISTORY], &mut [0.0; 8 * HISTORY], true),
+            Err(Error::NotSupported(_))
+        ));
+    }
+
+    #[test]
+    fn plastic_plane_strain_and_axisymmetric_tangents_are_the_derivative_of_the_forces() {
+        use crate::elastoplastic::{HISTORY, PLASTIC_STRAIN};
+        use fieldcraft_model::{Elastic, Hardening, Plastic};
+        let j2 = crate::plasticity::J2::from_model(
+            &Elastic::Isotropic { young: 1000.0, poisson: 0.3 },
+            Some(&Plastic::Bilinear { yield_stress: 5.0, hardening_modulus: 50.0, hardening: Hardening::Isotropic }),
+            None,
+        )
+        .unwrap();
+        let m = MaterialModel::Plastic(std::sync::Arc::new(j2));
+        let quad8 = [
+            [1.0, 0.0, 0.0],
+            [2.0, 0.1, 0.0],
+            [2.2, 1.0, 0.0],
+            [0.9, 1.1, 0.0],
+            [1.5, 0.0, 0.0],
+            [2.1, 0.55, 0.0],
+            [1.5, 1.1, 0.0],
+            [0.95, 0.5, 0.0],
+        ];
+        for theory in [Theory::PlaneStrain { thickness: 2.0 }, Theory::Axisymmetric] {
+            for (kind, x) in [(ElementKind::Quad8, &quad8[..]), (ElementKind::Tri3, &quad8[..3])] {
+                let plane = Plane::new(kind, x, Formulation::Full, theory).unwrap();
+                let points = plane.rule().points().len();
+                for kin in [Kinematics::Finite, Kinematics::Small] {
+                    let u = |scale: f64| -> Vec<[f64; 2]> { displacement(x.len(), scale).iter().map(|v| [v[0], v[1]]).collect() };
+                    let mut committed = vec![0.0; points * HISTORY];
+                    plane.nonlinear_points(&m, &u(0.03), kin, &vec![0.0; points * HISTORY], &mut committed, false).unwrap();
+                    assert!((0..points).all(|p| committed[p * HISTORY + PLASTIC_STRAIN] > 0.0), "{theory:?} {kind:?} {kin:?}");
+                    let at = u(0.06);
+                    let k = plane.nonlinear_points(&m, &at, kin, &committed, &mut committed.clone(), true).unwrap().0.tangent.unwrap();
+                    let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs()));
+                    let h = 1e-8;
+                    for j in 0..2 * x.len() {
+                        let shifted = |s: f64| {
+                            let mut v = at.clone();
+                            v[j / 2][j % 2] += s;
+                            plane.nonlinear_points(&m, &v, kin, &committed, &mut committed.clone(), false).unwrap().0.force
+                        };
+                        let (p, q) = (shifted(h), shifted(-h));
+                        for i in 0..2 * x.len() {
+                            let numeric = (p[i] - q[i]) / (2.0 * h);
+                            assert!(
+                                (numeric - k.get(i, j)).abs() < 1e-6 * scale,
+                                "{theory:?} {kind:?} {kin:?} {i} {j}: {numeric} vs {}",
+                                k.get(i, j)
+                            );
+                        }
+                    }
+                }
+            }
+        }
+        let stress = Plane::new(ElementKind::Quad4, &quad8[..4], Formulation::Full, Theory::PlaneStress { thickness: 1.0 }).unwrap();
+        assert!(matches!(
+            stress.nonlinear(&m, &[[0.0; 2]; 4], Kinematics::Small, &[0.0; 4 * HISTORY], &mut [0.0; 4 * HISTORY], true),
+            Err(Error::NotSupported(_))
+        ));
+    }
+
     #[test]
     fn rigid_rotation_strains_and_stresses_nothing() {
         let x20 = hex20();
@@ -501,7 +735,7 @@ mod tests {
         let neo = MaterialModel::NeoHookean { shear: 300.0, bulk: 3000.0 };
         for f in [Formulation::Full, Formulation::Hybrid { pressures: 4 }] {
             let solid = Solid::new(ElementKind::Hex20, &x20, f).unwrap();
-            let s = solid.finite_strain_stresses(&neo, &u, &[]).unwrap();
+            let s = solid.nonlinear_stresses(&neo, &u, Kinematics::Finite, &[]).unwrap();
             assert!(s.iter().all(|p| p.stress.iter().chain(&p.strain).all(|v| v.abs() < 1e-12 * 3000.0)), "{f:?}: {:?}", s[0]);
             let r = solid.finite_strain(&neo, &u, &[], &mut [], false).unwrap();
             assert!(r.force.iter().all(|v| v.abs() < 1e-11 * 3000.0));
diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs
index c2d15da..d561a1d 100644
--- a/crates/fem/src/lib.rs
+++ b/crates/fem/src/lib.rs
@@ -14,12 +14,13 @@
 //!   the solid kernels evaluate;
 //! - [`material`]: isotropic and orthotropic elasticity matrices, von Mises stress;
 //! - [`constitutive`]: the material point of finite-strain analysis (Kirchhoff stress and
-//!   spatial tangent from the deformation gradient, with history): St Venant–Kirchhoff and
-//!   neo-Hookean;
+//!   spatial tangent from the deformation gradient, with history): St Venant–Kirchhoff,
+//!   neo-Hookean and J2 plasticity;
 //! - [`element`]: the element table, stiffness (full, reduced, incompatible modes), consistent
 //!   and lumped mass, body forces, integration-point stresses and their extrapolation to nodes;
 //! - [`finite`]: total Lagrangian solid and 2D elements (internal forces, the consistent tangent
-//!   with its geometric part, mixed u–p at finite strain, Cauchy stress and logarithmic strain);
+//!   with its geometric part, mixed u–p at finite strain, Cauchy stress and logarithmic strain),
+//!   and the same with small strain for inelastic materials;
 //! - [`plane`]: 2D elements (plane stress, plane strain, axisymmetric) on tri3, tri6, quad4 and
 //!   quad8, with two displacement components per node;
 //! - [`shell`]: flat MITC4 quadrilateral and MITC3 triangle shells with a drilling rotation, six
@@ -33,6 +34,9 @@
 //!   shape vectors and element frequency bounds for explicit dynamics;
 //! - [`rotation`]: incrementally objective stress updates (Hughes–Winget Jaumann, Flanagan–Taylor
 //!   Green–Naghdi);
+//! - [`elastoplastic`]: J2 plasticity at a material point of implicit analysis, at finite strain
+//!   (Simo's exponential map in the principal logarithmic strains) and at small strain, with the
+//!   consistent tangents, on the radial return of [`plasticity`];
 //! - [`plasticity`]: J2 radial return with isotropic, kinematic and Johnson–Cook hardening,
 //!   Johnson–Cook damage;
 //! - [`truss`]: two- and three-node trusses (straight or curved), their stiffness, mass, body
@@ -52,6 +56,7 @@ pub mod beam;
 pub mod constitutive;
 pub mod dense;
 pub mod dof;
+pub mod elastoplastic;
 pub mod element;
 pub mod face;
 pub mod finite;
diff --git a/crates/fem/src/plasticity.rs b/crates/fem/src/plasticity.rs
index 8a36692..7f13484 100644
--- a/crates/fem/src/plasticity.rs
+++ b/crates/fem/src/plasticity.rs
@@ -30,6 +30,11 @@
 //! [`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.
 //!
+//! **One return map for both solvers.** [`J2::radial_return`] is the return above; explicit
+//! dynamics calls it through [`J2::update`] on its hypoelastic stress, the implicit material
+//! point ([`crate::elastoplastic`]) on its trial elastic strain, small or logarithmic, and adds
+//! the consistent tangent [`J2::tangent`].
+//!
 //! The strain increment the update takes is the **mechanical** one: a caller with initial or
 //! thermal strains subtracts their increment first.
 //!
@@ -98,6 +103,21 @@ impl PointState {
     }
 }
 
+/// What a radial return did: what the consistent tangent needs.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Return {
+    /// The equivalent plastic strain increment Δε_p.
+    pub plastic_strain: f64,
+    /// The trial von Mises stress of ξ* = s* − α, q* = √(3/2)‖ξ*‖.
+    pub trial: f64,
+    /// The flow direction ξ*/‖ξ*‖ (tensor components).
+    pub normal: [f64; 6],
+    /// The kinematic hardening modulus H_k.
+    pub kinematic: f64,
+    /// The slope of the flow stress against the plastic strain at the end of the increment.
+    pub slope: f64,
+}
+
 /// What one update did.
 #[derive(Debug, Clone, Copy, Default, PartialEq)]
 pub struct Increment {
@@ -232,27 +252,11 @@ impl J2 {
         }
         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;
-                }
+        if let Some(r) = self.radial_return(state, &mut trial, dt) {
+            let q_new = r.trial - (3.0 * self.shear + r.kinematic) * r.plastic_strain;
+            out = Increment { plastic_strain: r.plastic_strain, plastic_work: q_new * r.plastic_strain };
+            if self.heating > 0.0 && density > 0.0 {
+                state.temperature += self.heating * out.plastic_work / density;
             }
         }
         for k in 0..3 {
@@ -273,6 +277,68 @@ impl J2 {
         out
     }
 
+    /// The radial return of a trial deviator `trial` (tensor components) from `state`: when the
+    /// trial von Mises stress exceeds the yield stress, `trial` is moved back to the yield
+    /// surface, the back stress and the equivalent plastic strain of `state` are advanced, and
+    /// what the consistent tangent needs is returned; `None` when the point stays elastic.
+    pub fn radial_return(&self, state: &mut PointState, trial: &mut [f64; 6], dt: f64) -> Option<Return> {
+        let flow = self.flow.as_ref()?;
+        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) {
+            return None;
+        }
+        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 = 2.0 * self.shear * (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];
+        }
+        let rate = if dt > 0.0 { de / dt } else { 0.0 };
+        let (_, slope) = flow.stress(state.plastic_strain + de, rate, state.temperature, dt);
+        state.plastic_strain += de;
+        Some(Return { plastic_strain: de, trial: q, normal: n, kinematic: hk, slope })
+    }
+
+    /// The consistent (algorithmic) tangent of the return map, dσ/dε of the stress after the
+    /// update against the total strain of the increment (engineering-shear columns): the
+    /// elasticity K I ⊗ I + 2G ℙ while elastic, and after a return
+    /// K I ⊗ I + 2G θ ℙ − 2G θ̄ n ⊗ n with θ = 1 − 3G Δε_p / q*, θ̄ = 1/(1 + H/3G) − (1 − θ) and
+    /// H the slope of the flow stress plus the kinematic modulus at the end of the increment
+    /// (Simo & Hughes 1998, Box 3.2, eq. 3.3.10, in the equivalent plastic strain). `bulk`
+    /// leaves out the volumetric part K I ⊗ I when false (the isochoric part a mixed element
+    /// takes). `None` when 3G + H is not positive: the point softens faster than it unloads
+    /// elastically and has no tangent.
+    pub fn tangent(&self, r: Option<&Return>, bulk: bool) -> Option<[[f64; 6]; 6]> {
+        let g = self.shear;
+        let (theta, theta_bar, n) = match r {
+            None => (1.0, 0.0, [0.0; 6]),
+            Some(r) => {
+                let stiff = 3.0 * g + r.slope + r.kinematic;
+                if !(stiff > 0.0 && stiff.is_finite() && r.trial > 0.0) {
+                    return None;
+                }
+                let theta = 1.0 - 3.0 * g * r.plastic_strain / r.trial;
+                (theta, 3.0 * g / stiff - (1.0 - theta), r.normal)
+            }
+        };
+        let k = if bulk { self.bulk } else { 0.0 };
+        let identity = crate::constitutive::identity4();
+        let mut d = [[0.0; 6]; 6];
+        for a in 0..6 {
+            for b in 0..6 {
+                let (ia, ib) = (if a < 3 { 1.0 } else { 0.0 }, if b < 3 { 1.0 } else { 0.0 });
+                d[a][b] = k * ia * ib + 2.0 * g * theta * (identity[a][b] - ia * ib / 3.0) - 2.0 * g * theta_bar * n[a] * n[b];
+            }
+        }
+        Some(d)
+    }
+
     /// 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;
diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs
index 7cad70a..690a281 100644
--- a/crates/solver/src/nonlinear.rs
+++ b/crates/solver/src/nonlinear.rs
@@ -62,6 +62,7 @@
 use fieldcraft_fem::constitutive::MaterialModel;
 use fieldcraft_fem::dense::Matrix;
 use fieldcraft_fem::element::{Extrapolation, extrapolation};
+use fieldcraft_fem::finite::Kinematics;
 use fieldcraft_fem::quadrature::Rule;
 use fieldcraft_fem::{DofKind, DofSet, ElementType, PointStress, face, von_mises};
 use fieldcraft_model::{Extent, Load, Model, NonlinearStatic, Procedure, Step};
@@ -250,7 +251,7 @@ impl Problem<'_> {
         }
         let material = self.materials.get(e).ok_or_else(|| Error::Internal("element material".into()))?;
         let mut updated = committed.to_vec();
-        let r = kernel.finite_strain(material, &ue, committed, &mut updated, tangent).map_err(at)?;
+        let r = kernel.nonlinear(material, &ue, Kinematics::Finite, committed, &mut updated, tangent).map_err(at)?;
         Ok((r.force, r.tangent, updated))
     }
 
@@ -833,7 +834,7 @@ fn finish(
         let kernel = el.kernel()?;
         let points = if problem.finite {
             let material = problem.materials.get(i).ok_or_else(|| Error::Internal("element material".into()))?;
-            kernel.finite_strain_stresses(material, &ue, committed.get(i).map(Vec::as_slice).unwrap_or_default())
+            kernel.nonlinear_stresses(material, &ue, Kinematics::Finite, committed.get(i).map(Vec::as_slice).unwrap_or_default())
         } else {
             let per_node: Vec<[f64; 3]> = el.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect();
             kernel.stresses_with(&el.d, &per_node, None)
diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs
index 103c795..7d624eb 100644
--- a/crates/solver/src/system.rs
+++ b/crates/solver/src/system.rs
@@ -11,6 +11,7 @@ use fieldcraft_fem::beam::{Beam, BeamSection};
 use fieldcraft_fem::constitutive::MaterialModel;
 use fieldcraft_fem::dense::Matrix;
 use fieldcraft_fem::element::{Formulation, Solid, ThermalStrain};
+use fieldcraft_fem::finite::Kinematics;
 use fieldcraft_fem::plane::{Plane, Theory};
 use fieldcraft_fem::shell::{Shell, ShellSection};
 use fieldcraft_fem::truss::{Truss, young_modulus};
@@ -113,27 +114,30 @@ impl Kernel<'_> {
         }
     }
 
-    /// Finite-strain internal forces and tangent ([`fieldcraft_fem::finite`]) for the element's
-    /// displacements `u` (its degrees of freedom in order: 3 per node, or 2 for 2D elements).
-    pub fn finite_strain(
+    /// Nonlinear internal forces and tangent ([`fieldcraft_fem::finite`]) for the element's
+    /// displacements `u` (its degrees of freedom in order: 3 per node, or 2 for 2D elements),
+    /// total Lagrangian or with small strain.
+    pub fn nonlinear(
         &self,
         m: &MaterialModel,
         u: &[f64],
+        kin: Kinematics,
         committed: &[f64],
         updated: &mut [f64],
         tangent: bool,
     ) -> fieldcraft_fem::Result<fieldcraft_fem::finite::Finite> {
         match self {
-            Kernel::Solid(s) => s.finite_strain(m, u.as_chunks::<3>().0, committed, updated, tangent),
-            Kernel::Plane(p) => p.finite_strain(m, u.as_chunks::<2>().0, committed, updated, tangent),
+            Kernel::Solid(s) => s.nonlinear(m, u.as_chunks::<3>().0, kin, committed, updated, tangent),
+            Kernel::Plane(p) => p.nonlinear(m, u.as_chunks::<2>().0, kin, committed, updated, tangent),
         }
     }
 
-    /// Cauchy stress and logarithmic strain at the integration points, as [`Self::finite_strain`].
-    pub fn finite_strain_stresses(&self, m: &MaterialModel, u: &[f64], committed: &[f64]) -> fieldcraft_fem::Result<Vec<PointStress>> {
+    /// Cauchy stress and strain (logarithmic at finite strain) at the integration points, as
+    /// [`Self::nonlinear`].
+    pub fn nonlinear_stresses(&self, m: &MaterialModel, u: &[f64], kin: Kinematics, committed: &[f64]) -> fieldcraft_fem::Result<Vec<PointStress>> {
         match self {
-            Kernel::Solid(s) => s.finite_strain_stresses(m, u.as_chunks::<3>().0, committed),
-            Kernel::Plane(p) => p.finite_strain_stresses(m, u.as_chunks::<2>().0, committed),
+            Kernel::Solid(s) => s.nonlinear_stresses(m, u.as_chunks::<3>().0, kin, committed),
+            Kernel::Plane(p) => p.nonlinear_stresses(m, u.as_chunks::<2>().0, kin, committed),
         }
     }
 
diff --git a/crates/solver/tests/nonlinear.rs b/crates/solver/tests/nonlinear.rs
index 99a3578..69d5468 100644
--- a/crates/solver/tests/nonlinear.rs
+++ b/crates/solver/tests/nonlinear.rs
@@ -130,7 +130,7 @@ fn nonlinear_patch_test() {
     };
     let (e, nu) = (100.0, 0.3);
     let neo = fieldcraft_fem::constitutive::MaterialModel::from_model(&iso(e, nu), Some(&Hyperelastic::NeoHookean {})).unwrap();
-    let tau = neo.response(&f, &[], &mut []).unwrap().stress;
+    let tau = neo.response(&f, None, &[], &mut []).unwrap().stress;
     let j = fieldcraft_fem::dense::det3(&f);
     let hybrid = Formulation { hybrid: true, ..Formulation::default() };
     let mut cases: Vec<(ElementKind, Formulation, &str)> = variants().into_iter().filter(|v| v.1 != incompatible()).collect();