verification: implicit plasticity oracle cases and documentation

Four cases, each against its closed form: uniaxial tension of a J2 cube to a
logarithmic strain of 0.2 (Kirchhoff stress exact; 0.13 % from the explicit
solver's hypoelastic Cauchy stress, the elastic volume change between the two
laws), the Bauschinger effect of linear kinematic hardening on a load path
back to zero strain (exact), Hill's elastic-perfectly plastic thick cylinder
(bore displacement at a plastic radius of 1.5a within 0.009 %, collapse
pressure within 2e-8; the pending case is now active) and a cantilever's
plastic collapse, converged 3.1 % above beam theory's hinge, the bound a
continuum exceeds. docs/fem.md, the literature, ROADMAP (M15), parity and the
oracle page follow.

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

docs/fem.md

@@ -172,6 +172,22 @@ converges, so a cut-back increment leaves no trace. The materials:
 | Elastic (isotropic or orthotropic), no option | St Venant–Kirchhoff: S = D E, E = ½(FᵀF − I). Right for large rotations with small strains; it softens without bound in large compression. |
 | Isotropic elastic with `"hyperelastic": {"type": "neo_hookean"}` | Compressible neo-Hookean with the isochoric–volumetric split, W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)², μ and κ from E and ν (C₁₀ = μ/2, D₁ = 2/κ). Its small-strain limit is the linear material, which linear steps use. |
 
+| Isotropic elastic with `"plastic"` (bilinear isotropic or kinematic, table, Johnson–Cook) | J2 plasticity (`fieldcraft_fem::elastoplastic`): at finite strain Simo's multiplicative model (CMAME 99, 1992) with a Hencky elastic energy, τ = K ln J I + 2G dev ½ ln bₑ, the history Cₚ⁻¹ and the equivalent plastic strain, integrated by the exponential return map on the principal logarithmic strains; at small strain the additive radial return with the back stress. Both call `J2::radial_return`, the update of explicit dynamics, and its consistent tangent `J2::tangent`. |
+
+**Plasticity's tangent.** With Cₚ⁻¹ frozen over an iteration, L_v τ = 𝔻ᵉᵖ : ½ (∂ ln b/∂b) : (d b + b d) − (d τ + τ d), 𝔻ᵉᵖ the algorithmic tangent of the return
+(K I⊗I + 2Gθ ℙ − 2Gθ̄ n⊗n, Simo & Hughes 1998 Box 3.2) and the derivative of the logarithm the
+Daleckii–Krein divided differences θ_AB = artanh(r)/r, r = (λ²_A − λ²_B)/(λ²_A + λ²_B), in the
+principal frame of the trial bₑ: smooth through equal stretches, so uniaxial and undeformed
+states need no special case. The tangent is symmetric (isotropic hardening) and checked against
+central differences of the stress and of every element's forces from a plastic history. Hybrid
+elements take the deviatoric return and the volume K/2 (ln J)². Refused by name in nonlinear
+steps: kinematic hardening at finite strain (a back stress is not coaxial with bₑ, which the
+principal-frame return needs; it works with small strain), the Johnson–Cook strain-rate term (a
+static step has no time; C = 0 is the rate-independent law), and incompatible-mode elements of a
+plastic material. Damage, erosion and adiabatic heating act in explicit dynamics only (a warning
+says so). A hardening slope below −3G leaves the point without a tangent: the increment cuts
+back and the step stops at the minimum increment, naming the element.
+
 "Isotropic elastic" keeps meaning the same thing everywhere: the neo-Hookean model is chosen
 only by the explicit option, and explicit dynamics refuses a hyperelastic material by name.
 
@@ -190,9 +206,11 @@ differences of its forces.
 | plane stress | refused (σ_zz = 0 needs a local iteration for the thickness) |
 | trusses, shells, beams | refused (large rotations of structural elements come later in M15) |
 
-Without nonlinear geometry (`"nonlinear_geometry": false`) the elements are the linear ones:
-every material is elastic today, so such a step is the linear solution reached by Newton (the
-path the next materials, plasticity first, take with small strains).
+Without nonlinear geometry (`"nonlinear_geometry": false`) an elastic element is the linear one
+(K u), and an element of a plastic material takes the same loops with small strain
+(`Kinematics::Small`: ε = sym ∇u, the reference gradients, f = ∫ Bᵀσ dV, K = ∫ Bᵀ𝔻B dV with the
+consistent tangent, no geometric stiffness; hybrid elements project 1 + tr ε and K tr ε), so an
+elastic model is still the linear solution reached by Newton.
 
 **Follower pressure** (`fieldcraft_fem::face::follower_pressure`, `follower_edge_pressure`):
 fᵢ = −p ∫ Nᵢ (x,ξ × x,η) dξ dη on the current face, with its load stiffness
@@ -207,6 +225,23 @@ the rate of convergence of open loaded surfaces; it is revisited with the symmet
 (LDLᵀ) factorization that arc-length and Taylor–Hood need. Without nonlinear geometry a pressure
 is a dead load.
 
+**Load path.** A step's `load_path`, `[[fraction, factor], …]` from `[0, 0]` to fraction 1, makes
+the load factor piecewise linear in the step fraction (without one it equals the fraction), so
+loads and prescribed displacements can come down or reverse within a step: cyclic loading.
+Increments end at every point of the path, and the extrapolated predictor scales with the change
+of the load factor, so it turns with the path. LS-DYNA writes the path as the step's curve;
+Abaqus refuses it by name for now (it would be an amplitude on every load and boundary
+condition).
+
+**Predictors.** An increment starts from the linear extrapolation of the last one; with
+path-dependent materials an increment without one (the first, after a turn of the path to zero
+change) starts from the tangent predictor: one solve with the converged state's tangent for the
+change of the loads and of the prescribed displacements, so a prescribed motion is spread into the
+interior before the material sees it (started from the boundary's jump alone, the elements next
+to it would flow plastically far beyond the solution). Its work is then the energy criterion's
+reference. A first tangent that is not positive definite is a cutback when the undeformed
+model's tangent is (the strain of the increment, not a model that can move freely).
+
 **Newton–Raphson** with the consistent tangent refactored every iteration, an optional line
 search (off by default: on large elastic rotations the full step converges in fewer
 iterations), and convergence when the force residual (∞-norm over the largest force seen), the
@@ -219,8 +254,10 @@ tangent is not positive definite, or the iterations run out; the step stops with
 and the reason at the minimum increment, the increment cap, or `Limits::max_newton_iterations`.
 
 **Results**: the final displacements and reactions, the **Cauchy stress** and the
-**logarithmic strain** (ln V, engineering shears) at the integration points and averaged at the
-nodes, every increment's load factor, iterations and residuals, and a displacement frame per
+**logarithmic strain** (ln V, engineering shears; the small strain without nonlinear geometry) at
+the integration points and averaged at the nodes, the **equivalent plastic strain** at the
+integration points of plastic elements (`solve.probe` with `elements`, the largest in the step
+summary), every increment's load factor, iterations and residuals, and a displacement frame per
 converged increment (`results.frame`, the slider in the results panel).
 
 **Couplings** stay the linear constraints of the linear solver (small rotations of the
@@ -229,9 +266,11 @@ centrifugal loads are refused by name in nonlinear steps.
 
 **Decks**: Abaqus `*STEP, NLGEOM=YES|NO, INC=` with `*STATIC` (initial, period 1, minimum,
 maximum) and `*CONTROLS, PARAMETERS=LINE SEARCH`; the neo-Hookean option as `*HYPERELASTIC, NEO
-HOOKE`. LS-DYNA: the ramped implicit run with full Newton (`ilimit` 1, `maxref`) and
+HOOKE`; plasticity as for explicit dynamics (`*PLASTIC`); a load path, and a plastic deck of
+several general steps (Abaqus carries the plastic state from one to the next, every model step
+starts undeformed), are refused. LS-DYNA: the ramped implicit run with full Newton (`ilimit` 1, `maxref`) and
 `*CONTROL_IMPLICIT_AUTO`; hyperelastic materials in nonlinear steps and steps without nonlinear
-geometry are refused. Nastran refuses nonlinear steps (SOL 400) by name. Iteration caps and
+geometry are refused; a load path is the step's curve. Nastran refuses nonlinear steps (SOL 400) by name. Iteration caps and
 tolerances are left to each solver's own controls, whose norms differ from ours.
 
 **Verification** (`crates/solver/tests/nonlinear.rs`, `crates/fem/src/{constitutive,finite,face}.rs`,
@@ -242,4 +281,16 @@ test on a distorted mesh (all solids, hybrid ones included); a 90° rigid rotati
 the elastica (Bisshopp–Drucker) at P L²/(E I) = 1, 2, 5, 10 with hex20 (0.2–0.5 %); the inflation of a thick neo-Hookean tube by a follower pressure (Rivlin's closed form);
 Newton's quadratic convergence; the linear limit; a finite-strain Cook's membrane is pending
 its transcribed reference; and the failures (an element inverted by
-crushing stops at the minimum increment, naming it; the caps; what is refused).
+crushing stops at the minimum increment, naming it; the caps; what is refused). Plasticity
+(`crates/solver/tests/plasticity.rs`, `crates/fem/src/elastoplastic.rs`, `tests/oracle/`
+`nl-j2-uniaxial`, `j2-bauschinger`, `plastic-thick-cylinder`, `plastic-hinge-cantilever`): the
+consistent tangents against central differences from a plastic history (material point and
+every element, finite and small strain, full and hybrid); uniaxial tension to a logarithmic
+strain of 0.2 in every element kind, hybrid ones included (Kirchhoff stress exact to 10⁻⁷, and
+0.13 % from explicit dynamics' hypoelastic Cauchy stress, the elastic volume change between the
+two laws); the Bauschinger effect on a load path back to zero strain (exact); Hill's thick
+cylinder, the bore displacement at a plastic radius of 1.5a (0.009 %) and the collapse pressure
+(2·10⁻⁸); a cantilever's collapse 3.1 % above beam theory's hinge, converged under refinement;
+Newton's order with the consistent tangent 1.9 (small strain) and 2.3 (finite strain) in a
+plastic cantilever; and the refusals and the softening material that stops at the minimum
+increment.

docs/oracle.md

@@ -2,7 +2,7 @@
 
 Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the policy is [verification.md](verification.md). Every case is our own model, built with FieldCraft commands from the problem's published description, and compared with the reference value its `expected.json` cites, within the stated tolerance and for the stated reason. Accepted results are frozen in each case's `fixture.json`: a later change that moves one beyond its band (10⁻⁹ of its size) fails CI even while the reference still passes. No reference value was produced with a commercial product.
 
-**142 of 148 quantities pass** (6 known failure(s)) in 40 active cases: 31 theory, 9 published benchmark, 0 experiment; 10 pending cases wait for capabilities still being built.
+**149 of 155 quantities pass** (6 known failure(s)) in 44 active cases: 35 theory, 9 published benchmark, 0 experiment; 9 pending cases wait for capabilities still being built.
 
 ## Theory: closed-form solutions, patch tests and mesh invariants
 
@@ -31,6 +31,8 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 | [free-beam-frequencies](../tests/oracle/free-beam-frequencies/) | sixth frequency (the last rigid-body mode) | 0 | ≈ 0 | ≈ 0 | ± 8.1264e-5 | pass | matches |
 | [free-beam-frequencies](../tests/oracle/free-beam-frequencies/) | first free–free bending frequency (mode 7) | 0.081264 | 0.080618 | -0.795% | 1.5% | pass | matches |
 | [free-beam-frequencies](../tests/oracle/free-beam-frequencies/) | total mass | 20 | 20 | < 1e-7% | 1e-12 rel. | pass | matches |
+| [j2-bauschinger](../tests/oracle/j2-bauschinger/) | axial stress back at zero strain after a pull to ε_a = 0.01 | -247.52 | -247.52 | < 1e-7% | 1e-6 rel. | pass | blessed |
+| [j2-bauschinger](../tests/oracle/j2-bauschinger/) | accumulated equivalent plastic strain | 0.016089 | 0.016089 | < 1e-7% | 1e-6 rel. | pass | blessed |
 | [j2-uniaxial-explicit](../tests/oracle/j2-uniaxial-explicit/) | axial Cauchy stress at a logarithmic strain of 0.2 | 643.56 | 643.57 | 0.000174% | 1% | pass | matches |
 | [j2-uniaxial-explicit](../tests/oracle/j2-uniaxial-explicit/) | equivalent plastic strain | 0.19678 | 0.19678 | -0.00000113% | 1% | pass | matches |
 | [kirsch-plate](../tests/oracle/kirsch-plate/) | σ_xx at the top of the hole (0, 1, 0) over the remote stress | 3 | 3.0121 | 0.402% | 2% | pass | matches |
@@ -57,6 +59,8 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip shortening δ_h/L at P L²/(E I) = 5 | 0.38763 | 0.38955 | 0.495% | 2% | pass | matches |
 | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip deflection δ_v/L at P L²/(E I) = 10 | 0.81061 | 0.81421 | 0.444% | 1% | pass | matches |
 | [nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/) | tip shortening δ_h/L at P L²/(E I) = 10 | 0.555 | 0.558 | 0.541% | 2% | pass | matches |
+| [nl-j2-uniaxial](../tests/oracle/nl-j2-uniaxial/) | axial Cauchy stress at a logarithmic strain of 0.2 | 642.74 | 642.74 | < 1e-7% | 1e-6 rel. | pass | blessed |
+| [nl-j2-uniaxial](../tests/oracle/nl-j2-uniaxial/) | equivalent plastic strain | 0.19678 | 0.19678 | < 1e-7% | 1e-6 rel. | pass | blessed |
 | [nl-neo-hookean-uniaxial](../tests/oracle/nl-neo-hookean-uniaxial/) | axial Cauchy stress at the corner (1, 1, 1), stretch λ = 2 | 8.0601 | 8.0601 | < 1e-7% | 1e-6 rel. | pass | matches |
 | [nl-neo-hookean-uniaxial](../tests/oracle/nl-neo-hookean-uniaxial/) | lateral displacement u_y at the corner (1, 1, 1), l − 1 | -0.18686 | -0.18686 | < 1e-7% | 1e-6 rel. | pass | matches |
 | [nl-patch-hex8](../tests/oracle/nl-patch-hex8/) | u_x at inner node 22 | 0.13583 | 0.13583 | < 1e-7% | 1e-8 rel. | pass | matches |
@@ -121,6 +125,9 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 | [patch-wedge6](../tests/oracle/patch-wedge6/) | σ_yy at the inner node (0.46, 0.55, 0.53) | 2 | 2 | < 1e-7% | 1e-9 rel. | pass | matches |
 | [patch-wedge6](../tests/oracle/patch-wedge6/) | σ_zz at the inner node (0.46, 0.55, 0.53) | 3 | 3 | < 1e-7% | 1e-9 rel. | pass | matches |
 | [patch-wedge6](../tests/oracle/patch-wedge6/) | σ_xy at the inner node | 0 | ≈ 0 | ≈ 0 | ± 1.0000e-9 | pass | matches |
+| [plastic-hinge-cantilever](../tests/oracle/plastic-hinge-cantilever/) | tip force at a tip deflection of 15 (4.5 times the first yield's) | -312.5 | -322.61 | -3.23% | 5% | pass | blessed |
+| [plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/) | bore displacement at the pressure that takes the plastic zone to c = 15 (a = 10, b = 20, σ_y = 250) | 0.024355 | 0.024358 | 0.0089% | 0.2% | pass | blessed |
+| [plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/) | radial force on the bore (per radian times 2π, height 1) at a bore displacement of 0.3, 24 times the first yield's: the collapse pressure 2k ln 2 = 200.094 times 2πa | 12572 | 12572 | 0.00000178% | 0.1% | pass | blessed |
 | [pure-bending-coupling](../tests/oracle/pure-bending-coupling/) | bending stress σ_xx at mid-span, top fibre (10, 0.5, 1) | 15 | 15 | -0.0000103% | 0.5% | pass | matches |
 | [pure-bending-coupling](../tests/oracle/pure-bending-coupling/) | rotation of the reference point about y | 0.3 | 0.30355 | 1.18% | 1% | known failure: 1.2% high: the distributing coupling weights every node of the end face equally, corner and mid-side nodes of hex20 alike, which is not the consistent linear traction; the end face deforms locally and the fitted rotation includes it. Area-consistent weights would remove it | — |
 | [rotating-cylinder-axisymmetric](../tests/oracle/rotating-cylinder-axisymmetric/) | radial stress σ_r on the axis (r = 0) | 0.38571 | 0.38608 | 0.0958% | 1% | pass | matches |
@@ -173,6 +180,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 - **[cdt-cylinder](../tests/oracle/cdt-cylinder/)** (theory): Constrained Delaunay mesh of a cylinder: volume within the chordal bound, valid elements. Source: mesh invariants: an inscribed polygon of sagitta s loses at most 2s/R of a circle's area (elementary geometry); constrained Delaunay tetrahedralization after Diazzi et al., ACM TOG 42(6) 2023. Tolerances — volume: the tessellation inscribes the cylinder within the chordal tolerance s = 0.05: the polygon section loses between 0 and 2s/R = 0.5% of the area; the reference is the middle of that band; invalid: every element of a valid mesh has positive volume
 - **[cdt-stepped-block](../tests/oracle/cdt-stepped-block/)** (theory): Constrained Delaunay mesh of a stepped block: exact volume, valid elements. Source: mesh invariants: a tetrahedralization of a polyhedron fills its volume exactly (Si, TetGen paper ACM TOMS 41(2) 2015 for the definition of a CDT; Diazzi et al., ACM TOG 42(6) 2023); the volume is 60·40·10 + 40·40·10 + 20·40·10 − 10·10·10. Tolerances — volume: a planar boundary is tessellated exactly and the constrained Delaunay mesh fills it exactly: only rounding; invalid: every element of a valid mesh has positive volume; undetermined: a linear tetrahedron's sign is exact
 - **[free-beam-frequencies](../tests/oracle/free-beam-frequencies/)** (theory): Natural frequencies of a free beam: rigid-body modes and the first bending mode. Source: Euler–Bernoulli free–free beam, β₁L = 4.7300; Timoshenko, Young and Weaver, Vibration Problems in Engineering, 4th ed., Wiley 1974, §5.3. Tolerances — rigid_6: a free body has six zero frequencies; numerical rigid-body frequencies below 10⁻³ of the first elastic one (eigenvalues below 10⁻⁶ of it) are rounding of the shifted solve; f7: Euler–Bernoulli neglects shear and rotary inertia, about 0.8% at L/h = 20 for the free–free mode; mass: the mass matrix integrates the density exactly
+- **[j2-bauschinger](../tests/oracle/j2-bauschinger/)** (theory): Bauschinger effect of linear kinematic hardening: tension, then back to zero strain. Source: J2 plasticity with Prager's linear kinematic hardening in uniaxial stress (W. Prager, Ing.-Arch. 1955; Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3): from σ_a = σ_y + E_t (ε_a − σ_y/E) the material unloads elastically by 2σ_y and then hardens with E_t, so at zero strain σ = σ_a − 2σ_y − E_t (ε_a − 2σ_y/E). Tolerances — residual_stress: a homogeneous small-strain state the element represents exactly, with increments that end at the load path's turn; isotropic hardening would give −281.835 (it yields again only at −σ_a); plastic_strain: (σ_a − σ_y)/H forward plus E/(E + H) (ε_a − 2σ_y/E) in reverse
 - **[j2-uniaxial-explicit](../tests/oracle/j2-uniaxial-explicit/)** (theory): Uniaxial tension of a J2 bilinear-hardening cube to large strain (explicit). Source: von Mises plasticity with linear isotropic hardening under uniaxial stress (Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3); logarithmic strain for a hypoelastic–plastic law under coaxial stretching. Tolerances — sigma: uniaxial J2 with linear hardening: σ = σ_y + E_t (ε − σ_y/E) in logarithmic strain; the stretch takes 10⁴ wave transits, so inertia and the ringing of the start (ρcv ≈ 0.1% of σ) stay well under 1%; plastic_strain: ε_p = ε − σ/E, as for the stress
 - **[kirsch-plate](../tests/oracle/kirsch-plate/)** (theory): Plate with a circular hole in tension (Kirsch). Source: Kirsch, Zeitschrift des VDI 42 (1898) 797–807; Timoshenko and Goodier, Theory of Elasticity, 3rd ed., 1970, §35. Tolerances — kt: the reference is the infinite plate; the finite plate (W = 60a, Heywood's width correction about +0.1%, and the loaded ends 30a away) raises the peak by a few tenths of a percent, and the stress extrapolated to the hole surface adds up to about 1%; side: infinite-plate value −S; finite width and extrapolation as for Kt, on a smaller value
 - **[kirsch-plate-plane-stress](../tests/oracle/kirsch-plate-plane-stress/)** (theory): Plate with a circular hole in tension, 2D plane stress (Kirsch). Source: Kirsch, Zeitschrift des VDI 42 (1898) 797–807; Timoshenko and Goodier, Theory of Elasticity, 3rd ed., 1970, §35. Tolerances — kt: the reference is the infinite plate; the finite plate (W = 60a, loaded ends 30a away) raises the peak by a few tenths of a percent, and the stress extrapolated to the hole surface adds up to about 1%; side: infinite-plate value −S; finite width and extrapolation as for Kt, on a smaller value
@@ -189,6 +197,7 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 - **[nafems-le10](../tests/oracle/nafems-le10/)** (benchmark): NAFEMS LE10 thick plate under pressure (hex20). Source: NAFEMS, The Standard NAFEMS Benchmarks, TNSB Rev. 3, October 1990, test LE10 (target σ_yy = −5.38 MPa at D); as reproduced in the Abaqus 2017 Benchmarks Guide, 'LE10: Thick plate under pressure'. Tolerances — syy_D: the target is published to three figures; at D the stress is extrapolated to a corner where the hole edge meets the loaded face, and published converged solutions scatter by 1–2% (Abaqus Benchmarks Guide fine meshes: −5.51 to −5.89)
 - **[nl-cylinder-inflation](../tests/oracle/nl-cylinder-inflation/)** (theory): Inflation of a thick neo-Hookean cylinder by a follower pressure (plane strain). Source: Finite inflation of an incompressible neo-Hookean thick-walled tube in plane strain: r² = R² + a² − A², σ_θ − σ_r = μ (λ² − λ⁻²) with λ = r/R, integrated through the wall: p = (μ/2)[(ln λ_a² − λ_a⁻²) − (ln λ_b² − λ_b⁻²)] (R. S. Rivlin, Large elastic deformations of isotropic materials VI, Phil. Trans. R. Soc. A 242:173–195, 1949; R. W. Ogden, Non-Linear Elastic Deformations, Dover 1997, §5.3.2). Tolerances — ua: the closed form is incompressible and our material nearly so (κ/μ = 1000, a volume change of order p/κ ≈ 4·10⁻⁴); hybrid hex20, 6 × 16, follow the wall's strain gradient well inside; ub: as at the bore
 - **[nl-elastica-cantilever](../tests/oracle/nl-elastica-cantilever/)** (theory): Large-deflection cantilever under a dead tip load against the elastica (solids). Source: The elastica of a cantilever under a vertical end load: K. E. Bisshopp, D. C. Drucker, Large deflection of cantilever beams, Quarterly of Applied Mathematics 3(3):272–275, 1945 (elliptic-integral solution); values computed by integrating θ″ = −(P/EI) cos θ by shooting (our own RK4, 20000 steps), agreeing with the table of K. Mattiasson, Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals, IJNME 17(1):145–153, 1981, to all six digits. Tolerances — v1: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u1: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v2: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u2: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v5: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u5: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's; v10: hex20, 40 × 2 through the depth, L/h = 20: shear deformation (0.6 (h/L)² ≈ 0.15 % of the deflection) and the clamp's local restraint, both absent from the elastica; u10: as δ_v; the shortening is the smaller, second-order quantity, so its relative error is about twice the deflection's
+- **[nl-j2-uniaxial](../tests/oracle/nl-j2-uniaxial/)** (theory): Uniaxial tension of a J2 bilinear-hardening cube to large strain (implicit, finite strain). Source: J2 plasticity with linear isotropic hardening in uniaxial stress (Simo and Hughes, Computational Inelasticity, Springer 1998, §2.3) in the logarithmic strain, which Simo's exponential return map (CMAME 99:61–112, 1992) integrates exactly on a proportional path: Kirchhoff stress τ = σ_y + E_t (ε − σ_y/E), Cauchy stress σ = τ/J with ln J = τ/(3K) (Hencky's volume law). Tolerances — sigma: a homogeneous deformation the element represents exactly, reached on a proportional path: the tolerance only admits the Newton tolerance and rounding; plastic_strain: ε_p = (τ − σ_y)/H, as for the stress
 - **[nl-neo-hookean-uniaxial](../tests/oracle/nl-neo-hookean-uniaxial/)** (theory): Uniaxial stretch of a neo-Hookean block to twice its length. Source: Closed form of the compressible neo-Hookean model W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)² in uniaxial stress (G. A. Holzapfel, Nonlinear Solid Mechanics, Wiley 2000, §6.5, eq. 6.148 with the volumetric term of §6.4): σ_yy = μ J^(−5/3)(l² − Ī/3) + κ(J − 1) = 0 fixes the lateral stretch l, then σ_xx = μ J^(−5/3)(λ² − l²). Tolerances — sigma: a homogeneous deformation every element represents exactly: the tolerance only admits rounding and the Newton tolerance; lateral: a homogeneous deformation every element represents exactly: the tolerance only admits rounding and the Newton tolerance
 - **[nl-patch-hex8](../tests/oracle/nl-patch-hex8/)** (theory): Nonlinear patch test, hex8: a homogeneous finite deformation of a distorted mesh. Source: The patch test (B. M. Irons, A. Razzaque, in A. K. Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, 557–587) at finite strain: a boundary moved by x = F X makes F the solution everywhere (J. Bonet, R. D. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis, 2nd ed., Cambridge 2008, §8); the neo-Hookean Cauchy stress σ = μ J^(−5/3) dev b + κ (J − 1) I of that F. Tolerances — ux: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; uy: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; uz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; sxx: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; syy: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; szz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; sxy: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; sxz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance; syz: a patch test is passed exactly: the tolerance admits rounding and the Newton tolerance
 - **[nl-rigid-rotation](../tests/oracle/nl-rigid-rotation/)** (theory): A quarter turn strains and stresses nothing (objectivity). Source: Objectivity of the total Lagrangian formulation and of the hyperelastic stress: a rigid motion x = R X + c gives C = I and σ = 0 whatever the rotation (Bonet & Wood 2008, §3.4 and §5.2; T. Belytschko, W. K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, 2nd ed., Wiley 2014, §5.5). Tolerances — ux: the inner nodes follow the rigid rotation exactly; uy: the inner nodes follow the rigid rotation exactly; uz: the inner nodes follow the rigid rotation exactly; von_mises: zero up to rounding of stresses of order E = 1000: a 90° turn leaves no strain
@@ -199,6 +208,8 @@ Generated by `cargo xtask oracle` from [`tests/oracle/`](../tests/oracle/); the
 - **[patch-tet10](../tests/oracle/patch-tet10/)** (theory): Constant-stress patch test, tet10. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1)
 - **[patch-tet4](../tests/oracle/patch-tet4/)** (theory): Constant-stress patch test, tet4. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1)
 - **[patch-wedge6](../tests/oracle/patch-wedge6/)** (theory): Constant-stress patch test, wedge6. Source: Irons and Razzaque, 'Experience with the patch test for convergence of finite elements', in Aziz (ed.), The Mathematical Foundations of the Finite Element Method, Academic Press 1972, pp. 557–587; MacNeal and Harder, Finite Elements in Analysis and Design 1 (1985) 3–20, constant-stress loading, p. 14 Table 7. Tolerances — ux: a patch test is passed exactly: the tolerance only admits rounding; uy: a patch test is passed exactly: the tolerance only admits rounding; uz: a patch test is passed exactly: the tolerance only admits rounding; sxx: a patch test is passed exactly: the tolerance only admits rounding; syy: a patch test is passed exactly: the tolerance only admits rounding; szz: a patch test is passed exactly: the tolerance only admits rounding; sxy: a patch test is passed exactly: the tolerance only admits rounding (stresses of order 1)
+- **[plastic-hinge-cantilever](../tests/oracle/plastic-hinge-cantilever/)** (theory): Plastic collapse of a cantilever: the fully plastic moment at the root. Source: Limit analysis (B. G. Neal, The Plastic Methods of Structural Analysis, 3rd ed., Chapman & Hall 1977, ch. 1–2): beam theory's plastic hinge at the root, the fully plastic moment M_p = σ_y b h²/4 of a rectangular section, collapses a cantilever of length L at P_c = M_p/L = 312.5 (b = h = 10, L = 200, σ_y = 250); a bound for the solid, whose collapse lies above it (W. Prager, P. G. Hodge, Theory of Perfectly Plastic Solids, Wiley 1951: a continuum mechanism allows no normal velocity jump, which the beam's hinge rotation is). Tolerances — collapse_force: the solid's collapse lies above beam theory's hinge: refined through the depth (4 to 8 elements) and along the length (40 to 160) the force settles 3.1 % above M_p/L (322.3), so the band is that bound's gap; the root is a symmetry plane (no rotation, free lateral contraction) held on average through a coupling, the tip pushed through another; reduced hex20, 20 × 2 × 4
+- **[plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/)** (theory): Elastic–perfectly plastic thick cylinder: plastic zone and collapse pressure. Source: Hill, The Mathematical Theory of Plasticity, Oxford 1950, ch. V §2: a von Mises cylinder in plane strain, elastically incompressible, with the plastic zone a ≤ r ≤ c under the pressure p = k (2 ln(c/a) + 1 − c²/b²), k = σ_y/√3, and the bore displacement u_a = k c²/(2 G a); collapse when c = b at p = 2k ln(b/a). Tolerances — u_bore: 20 reduced quad8 elements through the wall resolve the logarithmic plastic field and ν = 0.4999 stands in for Hill's incompressible elasticity (0.009 % apart); u_a ∝ c², so 0.2 % on u_a is 0.1 % on the plastic radius; collapse_force: once the wall is plastic through, perfect plasticity holds the force at the collapse value (σ_θ − σ_r = 2k everywhere, σ_z the mean of the others as plastic flow dominates); the band admits the quadrature of the logarithmic field
 - **[pure-bending-coupling](../tests/oracle/pure-bending-coupling/)** (theory): Pure bending through a distributing coupling. Source: Euler–Bernoulli pure bending, σ = M_y z/I, θ = ML/(EI) (exact in elasticity for ν = 0); Timoshenko and Goodier, Theory of Elasticity, 3rd ed., 1970, §102. Tolerances — sigma_top: with ν = 0 pure bending is a quadratic displacement field that hex20 represents exactly; the coupling's nodal force pattern differs from the linear traction only within a depth or two of the end (Saint-Venant), five depths from the probe; rotation: the reference point follows the best-fit rigid motion of the end section, which the end zone's local deformation disturbs by well under 1%
 - **[rotating-cylinder-axisymmetric](../tests/oracle/rotating-cylinder-axisymmetric/)** (theory): Long rotating solid cylinder, axisymmetric. Source: Timoshenko and Goodier, Theory of Elasticity, 3rd ed., McGraw-Hill 1970, §32 (long rotating cylinder, ε_z = 0). Tolerances — radial_axis: stress extrapolated to the axis, where σ_r = σ_θ: under 1% with sixteen quadratic elements across; hoop_rim: stress extrapolated to the surface: under 1%
 - **[rotating-disc-plane-stress](../tests/oracle/rotating-disc-plane-stress/)** (theory): Rotating disc with a central hole, plane stress. Source: Timoshenko and Goodier, Theory of Elasticity, 3rd ed., McGraw-Hill 1970, §32 (rotating disc of uniform thickness with a hole). Tolerances — hoop_bore: stress extrapolated to the bore surface, where the a²b²/r² term is steepest: 16 radial quadratic elements graded towards the bore leave about 0.2% (8 uniform: −2.4%; 32 graded: −0.05%); radial_max: the nearest node to √(ab) is within half an element of it, where σ_r is stationary (an error of order (h/r)², under 0.5%), and nodal averaging adds a few tenths of a percent; hoop_rim: as at the bore
@@ -218,6 +229,5 @@ Reference values already cited; each case gets its recipe when the capability it
 | [nl-cook-membrane](../tests/oracle/nl-cook-membrane/) | benchmark | Cook's membrane, nearly incompressible neo-Hookean, large load (plane strain) | the published finite-strain reference value, transcribed from the paper with its load and material constants | vertical displacement of the top corner C (48, 60) = to be transcribed | R. D. Cook, Improved two-dimensional finite element, J. Struct. Div. ASCE 100(ST9):1851–1863, 1974 (the tapered panel (0, 0), (48, 44), (48, 60), (0, 44), clamped on x = 0, shear on x = 48); finite-strain versions with a neo-Hookean material are published by J. C. Simo, F. Armero, IJNME 33:1413–1449, 1992; the values to compare with are to be transcribed from the paper |
 | [pinched-cylinder](../tests/oracle/pinched-cylinder/) | benchmark | Pinched cylinder with rigid diaphragms | shell elements (M8) | radial displacement under the load = -1.8248e-5 | Belytschko, Stolarski, Liu, Carpenter and Ong, 'Stress projection for membrane and shear locking in shell finite elements', Computer Methods in Applied Mechanics and Engineering 51 (1985) 221–258 (the shell obstacle course); reference 1.8248 × 10⁻⁵ under the load |
 | [pinched-hemisphere](../tests/oracle/pinched-hemisphere/) | benchmark | Hemispherical shell with an 18° hole under alternating point loads | shell elements (M8), or mapped blocks that follow doubly curved faces | displacement under the load along the load = 0.094 | MacNeal and Harder, 'A proposed standard set of problems to test finite element accuracy', Finite Elements in Analysis and Design 1 (1985) 3–20, Fig. 9 (p. 9), Table 5(b) (p. 11): 0.0940 (lower bound 0.0924 from their [7]) |
-| [plastic-thick-cylinder](../tests/oracle/plastic-thick-cylinder/) | theory | Elastic–perfectly plastic thick cylinder: collapse pressure | implicit (static) J2 plasticity: explicit J2 is in, a collapse load needs the static solver | collapse pressure for a = 10, b = 20, σ_y = 250 = 200.09 | Hill, The Mathematical Theory of Plasticity, Oxford 1950, ch. V §2: plane strain von Mises collapse pressure p = (2/√3)σ_y ln(b/a) |
 | [scordelis-lo-roof-shell](../tests/oracle/scordelis-lo-roof-shell/) | benchmark | Scordelis–Lo roof with shell elements | shell elements (MITC4, M8) | vertical displacement at the midpoint of the free edge = -0.3024 | MacNeal and Harder, 'A proposed standard set of problems to test finite element accuracy', Finite Elements in Analysis and Design 1 (1985) 3–20, Fig. 8 (p. 8), Table 5(a) (p. 11); Belytschko, Stolarski, Liu, Carpenter and Ong, 'Stress projection for membrane and shear locking in shell finite elements', Computer Methods in Applied Mechanics and Engineering 51 (1985) 221–258 (the shell obstacle course) |
 | [taylor-bar](../tests/oracle/taylor-bar/) | experiment | Taylor bar impact (copper) | the measured final shapes, to be transcribed from the paper (explicit Johnson–Cook plasticity is in) | final over initial length = to be transcribed | Wilkins and Guinan, 'Impact of cylinders on a rigid boundary', Journal of Applied Physics 44 (1973) 1200–1206 (final length and mushroom diameter vs impact speed) |

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 | | started (von Mises plasticity: bilinear, table, Johnson–Cook; Johnson–Cook damage, erosion; written to Abaqus `*PLASTIC`, Nastran `MATS1`, LS-DYNA plastic materials where expressible; compressible neo-Hookean for nonlinear steps, Abaqus `*HYPERELASTIC, NEO HOOKE`) |
+| Model | Materials: plasticity, hyperelasticity, damage, thermal | M8 | | started (von Mises plasticity: bilinear, table, Johnson–Cook; Johnson–Cook damage, erosion; plasticity in explicit and nonlinear static steps; written to Abaqus `*PLASTIC`, Nastran `MATS1`, LS-DYNA plastic materials where expressible; compressible neo-Hookean for nonlinear steps, Abaqus `*HYPERELASTIC, NEO HOOKE`) |
 | 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 |
@@ -53,7 +53,7 @@ products are run, no documentation text is copied). Status: **done**, **started*
 | Model | Loads: force and moment, pressure, gravity, body force | M3 | | done |
 | Model | Reference points, kinematic and distributing coupling | M3 | | done |
 | Model | Steps: linear static, frequency; output presets | M3 | | done |
-| Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8, M15 | | started (nonlinear static with nonlinear geometry, increments and tolerances, written to Abaqus `*STEP, NLGEOM` and LS-DYNA implicit, refused by Nastran; explicit dynamics; heat transfer; amplitudes planned) |
+| Model | Steps: nonlinear static, explicit dynamics, heat transfer; amplitudes | M7–M8, M15 | | started (nonlinear static with or without nonlinear geometry, increments, tolerances and a load path (cyclic loading), written to Abaqus `*STEP, NLGEOM` and LS-DYNA implicit (the load path as LS-DYNA's curve; refused by Abaqus for now), refused by Nastran; explicit dynamics; heat transfer; amplitudes per load planned) |
 | Interaction | Tie between non-matching meshes | M7 | | planned |
 | Interaction | Contact (general and pairs), friction, automatic pair detection | M7 | | planned |
 | Interaction | Springs and dashpots (ground or two-point; DOF, vector or line), point masses with rotary inertia | M4 | | done (implicit; explicit translational; all three decks with what each cannot take reported) |
@@ -63,7 +63,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); quad4/tri3 MITC shells and Timoshenko beams with profiles, mixed with solids; reference points, kinematic and distributing couplings, moments; centrifugal loads) |
-| Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff and neo-Hookean; verified by the elastica, tube inflation, patch and objectivity tests; plasticity, shells, contact, arc-length next) |
+| Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff, neo-Hookean and J2 plasticity (Simo's exponential return map at finite strain, radial return at small strain, consistent tangents, the return shared with explicit dynamics); verified by the elastica, tube inflation, patch and objectivity tests, uniaxial tension to large strain, Hill's thick cylinder, a cantilever's collapse and the Bauschinger effect; Mooney–Rivlin, thermal strain, shells, contact, arc-length next) |
 | 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 |

docs/solvers/literature.md

@@ -1186,3 +1186,28 @@ update are written from (cited again in the module docs):
   transcribed finite-strain values), R. D. Cook, J. Struct. Div. ASCE 100(ST9):1851–1863, 1974,
   and J. C. Simo, F. Armero, IJNME 33:1413–1449, 1992; tube inflation, R. S. Rivlin, Phil. Trans. R. Soc. A
   242:173–195, 1949, and R. W. Ogden, *Non-Linear Elastic Deformations*, Dover 1997.
+
+## Part N: Implicit plasticity (finite and small strain), as implemented
+
+- Multiplicative J2 plasticity with a Hencky elastic energy and the exponential return map in
+  the principal logarithmic strains, which keeps the classical radial return of the
+  infinitesimal theory: J. C. Simo, *Algorithms for static and dynamic multiplicative plasticity
+  that preserve the classical return mapping schemes of the infinitesimal theory*, CMAME
+  99(1):61–112, 1992, DOI [10.1016/0045-7825(92)90123-2](https://doi.org/10.1016/0045-7825(92)90123-2);
+  J. C. Simo, T. J. R. Hughes, *Computational Inelasticity*, Springer 1998, ch. 9.
+- The radial return (shared with explicit dynamics) and its consistent tangent
+  K I⊗I + 2Gθ ℙ − 2Gθ̄ n⊗n: Simo & Hughes 1998, §3.3, Box 3.2 (A4.1); J. C. Simo, R. L. Taylor,
+  *Consistent tangent operators for rate-independent elastoplasticity*, CMAME 48(1):101–118,
+  1985, DOI [10.1016/0045-7825(85)90070-2](https://doi.org/10.1016/0045-7825(85)90070-2).
+- The derivative of the matrix logarithm by divided differences (Daleckii–Krein), smooth through
+  equal stretches: R. Bhatia, *Matrix Analysis*, Springer 1997, §V.3.
+- Kinematic hardening at finite strain (refused for now; the route when it comes): additive
+  plasticity in the Lagrangian logarithmic strain space, C. Miehe, N. Apel, M. Lambrecht, CMAME
+  191(47–48):5383–5425, 2002, DOI [10.1016/S0045-7825(02)00438-3](https://doi.org/10.1016/S0045-7825(02)00438-3).
+- Prager's linear kinematic hardening: W. Prager, *The theory of plasticity: a survey of recent
+  achievements*, Proc. Inst. Mech. Eng. 169:41–57, 1955.
+- Verification references: the thick cylinder's plastic zone and collapse pressure, R. Hill,
+  *The Mathematical Theory of Plasticity*, Oxford 1950, ch. V; the cantilever's plastic hinge,
+  B. G. Neal, *The Plastic Methods of Structural Analysis*, 3rd ed., Chapman & Hall 1977, and the
+  continuum's velocity fields, W. Prager, P. G. Hodge, *Theory of Perfectly Plastic Solids*,
+  Wiley 1951.