OneCAD-Simulation/docs/fem.md

# FEM foundation: degrees of freedom and shape functions

What every element in `crates/fem` and every solver builds on. Code: `crates/fem/src/dof.rs`,
`crates/fem/src/basis.rs`, `crates/fem/src/jet.rs`.

## Degrees of freedom

Every unknown is a pair (node, `DofKind`):

| `DofKind` | Meaning | Model / deck number |
|---|---|---|
| `Ux`, `Uy`, `Uz` | translations along X, Y, Z | 1, 2, 3 |
| `Rx`, `Ry`, `Rz` | rotations about X, Y, Z | 4, 5, 6 |
| `P` | pressure (mixed u–p elements) | 8 |
| `T` | temperature | 11 |

The enum order (Ux, Uy, Uz, Rx, Ry, Rz, T, P) is the numbering order within a node.
`DofKind::from_number` / `number` convert to the 1-based numbers boundary conditions use;
`from_direction` / `direction` to the structural directions 0–5.

**Elements declare their unknowns** with `ElementDofs { node: DofSet, internal: Internal }`:
the kinds every node of the element carries (`DofSet::TRANSLATIONS` for solids, `PLANE` for
2D elements, `STRUCTURAL` for shells and beams, `TEMPERATURE` for heat transfer, translations
plus `PRESSURE` for continuous-pressure mixed elements) and the internal unknowns the element
condenses out itself (`incompatible_modes`, discontinuous `pressure`), which never reach the
global system. An element's matrices are ordered node by node, its `node` kinds in `DofKind`
order within a node (a shell's 24 × 24 stiffness: node 1 Ux Uy Uz Rx Ry Rz, node 2 …). Solid
elements: `Formulation::dofs()` / `Solid::dofs()`.

**Numbering** (`DofLayout`): collect what every node needs with `DofLayout::builder()`
(`element(labels, dofs)`, `node(label, kinds)` for constraints and boundary conditions that
need a kind no element gives, `point(label, kinds)` for reference points), then `build()`. A
node carries the union of what was asked for. Unknowns are numbered mesh nodes in label order,
each node's kinds in `DofKind` order, then the points in the order added; a slot's unknowns are
contiguous, so a node pair is one dense block in the assembled matrix. The same model always
gives the same numbering.

Lookups: `node(label)` / `slot(label)` (a slot is a position: mesh nodes `0..nodes()`, points
`point_slot(k)`), `dof(slot, kind)`, `range(slot)`, `kinds(slot)`, `locate(dof) -> (slot,
kind)`, `element_dofs(slots, kinds)`, `describe(dof)`, `stride()` (uniform layouts).

Trusses declare `fieldcraft_fem::truss::DOFS` (translations, on `Basis::Line2`/`Line3`).
Point masses, springs and dashpots are not mesh elements: the implicit solver asks the layout
for what each needs at its nodes (`fieldcraft_solver::discrete::needs`: translations for a mass
or a translational spring, rotations for an inertia or a torsional spring; a reference point
they use is a point with six) and finds its unknowns with `discrete::dof_of`.

Solvers: `fieldcraft_solver::system::Mesh::layout` numbers the implicit system (mesh nodes,
then six per coupled reference point); assembly (`sparse::add_block(k, layout, slots, kinds,
data)`, `add_scalar`, `add_diagonal`) works for any mix of kinds per node, so a solid sharing
nodes with a shell adds its translation block into the shell nodes' six-unknown blocks.
Boundary conditions and loads address (node, kind): a boundary condition holds the kinds a node
carries, and a moment needs rotations. The explicit solver has solid nodes only and checks that
its layout is the uniform node-major one (`NODE_DOFS = 3`).

## Shape functions

`Basis` is every family's interpolation:

| Family | Reference element | Nodes |
|---|---|---|
| `Line2`, `Line3` | ξ ∈ [−1, 1] | ends −1, 1; line3 adds the middle |
| `Tri3`, `Tri6` | ξ, η ≥ 0, ξ + η ≤ 1 | corners (0,0), (1,0), (0,1); mid-sides of 1-2, 2-3, 3-1 |
| `Quad4`, `Quad8`, `Quad9` | [−1, 1]² | corners counter-clockwise from (−1,−1); mid-sides of 1-2, 2-3, 3-4, 4-1; centre |
| `Tet4`, `Tet10` | unit tetrahedron | corners (0,0,0), (1,0,0), (0,1,0), (0,0,1); mid-sides of 1-2, 2-3, 3-1, 1-4, 2-4, 3-4 |
| `Wedge6`, `Wedge15` | triangle × ζ ∈ [−1, 1] | bottom 1-3, top 4-6; mid-sides of the bottom edges, top edges, vertical edges |
| `Pyramid5` | base [−1, 1]² at ζ = 0, apex (0, 0, 1) | base corners as quad4, then the apex |
| `Hex8`, `Hex20`, `Hex27` | [−1, 1]³ | corners 1-4 at ζ = −1, 5-8 above; hex20 mid-sides (bottom, top, vertical edges); hex27 face centres in `ElementKind::Hex8.faces()` order, then the centre |

Node orders are the model's (`ElementKind::edges()` / `faces()` are the single source):
corners, mid-side nodes in edge order, face nodes, body nodes. `Basis::of(ElementKind)` maps a
model element; `nodes()` gives natural node coordinates; `dim()`, `count()`, `corners()`,
`degree()` (the complete polynomial degree reproduced: 1 for linear families including the
rational pyramid, 2 for the quadratic ones).

Evaluation at a natural point `xi: [f64; 3]` (unused coordinates are zero):

- `eval(xi) -> Shape`: values `n` and first derivatives `dn[i] = ∂Nᵢ/∂(ξ, η, ζ)`. tet4, tet10,
  wedge6, hex8, hex20 use the hand-written functions of `shape.rs` (the solid kernels' path).
- `eval2(xi) -> Shape2`: adds `d2n[i][a][b] = ∂²Nᵢ/∂ξₐ∂ξ_b`, from second-order forward-mode
  automatic differentiation of the same formulas (`jet::Jet`): exact to round-off.
- `gradients(&shape, xyz) -> Gradients`: physical `dx[i] = ∇ₓNᵢ = J⁻¹ ∇_ξNᵢ` and `det J`, with
  J_ab = ∂x_b/∂ξₐ; an inverted or collapsed element is `Error::Degenerate`.
- `hessians(&shape2, xyz) -> Hessians`: physical `d2x[i]` =
  J⁻¹ (H_ξ − Σ_b (∂N/∂x_b) X_b″) J⁻ᵀ with X_b″ = Σᵢ x_{i,b} H_ξ(Nᵢ), the map's second
  derivative (T. J. R. Hughes, *The Finite Element Method*, 2000, §3.9, §4.2), and
  `laplacian(i)`. Without the X″ term a field linear in x would have a nonzero Hessian on a
  curved element.

A family of dimension d maps onto the first d physical coordinates (plane elements in x–y,
lines along x). Elements embedded in a higher dimension (shells, beams, trusses in 3D) pass
node coordinates in their own local frame.

## Verification (`basis::tests`, `dof::tests`, `jet::tests`)

For every family: Kronecker delta at the nodes; partition of unity (Σ N = 1, Σ ∇N = 0,
Σ H = 0); exact reproduction of every monomial up to `degree()` with its gradient and Hessian
in natural coordinates; first and second derivatives against central differences;
on an affine map, physical gradients, Hessians and Laplacians exact for a complete polynomial
of the family's degree; on a curved map, linear fields exact with zero Hessian, and the
physical Hessian of a quadratic interpolant consistent with finite differences of its physical
gradient. The solver and explicit migration to `DofLayout` was checked bit for bit: every
result of the solver, explicit, engine and io test suites hashes identically before and after.

## Heat transfer and thermal strain

A heat conduction element carries one temperature per node (`heat::HEAT_DOFS`, kind `T`);
`heat::Conduction` gives the conductivity ∫ ∇Nᵀ k ∇N dV (full rule), the capacity ∫ ρc N Nᵀ dV
and a source ∫ N Q dV (mass rule), and the heat flux −k∇T at the integration points;
`heat::face_weights` and `face_film` integrate surface flux and convection on tri3/tri6/quad4/quad8
faces. The solver (`fieldcraft_solver::heat`) numbers the temperatures with
`DofLayout::uniform(…, DofSet::TEMPERATURE)` and assembles with the same `sparse::add_block`.

Thermal strain ε_th = α (T − T_ref) along the material axes enters a solid as the load
`Solid::thermal_load(d, &ThermalStrain { alpha, delta })` = ∫ Bᵀ D ε_th dV and its stress as
`Solid::stresses_with(d, u, Some(&thermal))` = D (ε − ε_th); `delta` is the nodal T − T_ref.
2D elements have the same pair (`Plane::thermal_load`, `Plane::stresses_with`): plane stress
leaves εzz free (σzz = 0), plane strain holds it (σzz = −D₃·ε_th), axisymmetric elements strain
the hoop direction too. Trusses (rod2, rod3) take the axial part Σ αₖ tₖ² ΔT along their tangent (`Truss::thermal_load`,
`Truss::stresses_with`). Shells and beams will provide the same two functions for their strains;
springs, dashpots and point masses have no thermal strain.

Explicit dynamics takes a step's temperature loads as a thermal strain rate per element,
α ΔT / end time (the change reached linearly over the step, like a prescribed displacement, so
there is no thermal shock), subtracted from each increment's strain increment before the
constitutive update, so it feeds J2 and Johnson–Cook plasticity alike; isotropic expansion only
(the rate is the same in any frame). Explicit trusses take no thermal strain yet (a warning
says so).

Model side: `Material::thermal` (`conductivity`, `specific_heat`, `expansion`, each isotropic
or per axis, and `reference_temperature`), `BoundaryCondition::Temperature`, loads
`SurfaceFlux`, `Convection`, `HeatSource` (heat transfer) and `Temperature` (a temperature field
in a static step), procedures `HeatTransfer { transient }` and `CoupledTemperatureDisplacement`.
Commands: `solve.heat`, and `solve.static` for coupled steps. The deck writers refuse heat
transfer by name for now.
## Hybrid (mixed u–p) solids

A section with `"formulation": {"hybrid": true}` gives `Formulation::Hybrid { pressures }`
(`crate::hybrid`): the isotropic elasticity splits as D = D_dev + K m mᵀ, the mean stress is an
independent field discontinuous between elements (hex20: 1, x, y, z; hex8 and tet10: constant),
and it is condensed out element by element, p̂ = K M⁻¹ Gᵀ u, leaving the stiffness
K_dev + K G M⁻¹ Gᵀ. `ElementDofs.internal.pressure` counts the condensed unknowns; nothing
global changes, and Cholesky still applies. Stresses are D_dev ε + p m with the element's
pressure field. The choice of pressure space follows the inf-sup analysis (docs in
`crates/fem/src/hybrid.rs`): Q2–P1disc and P2–P0 are stable; Q1–P0 does not lock but may show
checkerboard pressures; linear tetrahedra and wedges lock even so and are refused. A thermal
strain splits the same way: ∫ Bᵀ D_dev ε_th on the nodes and, through the pressure,
K G M⁻¹ ∫ N_p tr ε_th, so a heated nearly incompressible solid takes its bulk stress −3KαΔT
when confined. A continuous
nodal pressure (`DofKind::P`, Taylor–Hood) needs a symmetric indefinite factorization and is
planned.

## Symmetric indefinite factorization (LDLᵀ)

`fieldcraft_solver::ldlt` factors P A Pᵀ = L D Lᵀ (L unit lower triangular, D with 1 × 1 and
2 × 2 pivots) on the Cholesky factor's symbolic factorization: the same nested dissection
ordering, elimination tree, supernodes, fill, factor-size cap and parallel subtree tasks, so it
is deterministic on any number of threads. Pivots are chosen by the bounded Bunch–Kaufman (rook)
search of Ashcraft, Grimes & Lewis (1998) among the columns of the supernode being factored
(Schenk & Gärtner's supernode Bunch–Kaufman, 2006): an interchange inside a supernode changes no
row structure, so nothing is delayed and the static structure holds. Its limits, stated in the
module: L's entries below the supernode are not bounded by the search (`solve_refined` adds
iterative refinement), and a column whose supernode offers no pivot above round-off is a located
`Singular` error (a matrix singular at that column, or a zero diagonal coupled only to later
supernodes, which no FieldCraft system has: hybrid pressures are condensed per element).

The factor reports its **inertia** (positive and negative eigenvalue counts, Sylvester's law):
a 2 × 2 pivot of the search has a negative determinant (one of each). `Symmetric::new(a, cap,
kind)` picks Cholesky (`Definite`, the default everywhere a matrix is positive definite, so
linear results are bit for bit unchanged), LDLᵀ (`Indefinite`), or Cholesky with LDLᵀ when the
matrix turns out not to be positive definite (`Auto`, one symbolic factorization for both).
On a hex8 stiffness pattern of 140 000 equations (1.5 · 10⁸ entries) the LDLᵀ factor costs 1.1
times the Cholesky factor, and 1.2 times on a shift with 10 600 negative eigenvalues and 7 700
2 × 2 pivots (`ldlt::tests::speed_against_cholesky`).

## Linear buckling

A step `{"type": "buckle", "modes": 4}` (`Procedure::Buckle`, solved by `solve.buckle`,
`fieldcraft_solver::buckle`) solves its loads and prescribed motions as a linear static preload,
forms the initial-stress stiffness K_σ = ∫ (∇N_a · σ ∇N_b) I dV of every solid and 2D element
from the preload's stresses at the stiffness points (`Solid::geometric_stiffness`,
`Plane::geometric_stiffness`, with the axisymmetric hoop term), and finds the load factors λ
with (K + λ K_σ) φ = 0 by block Lanczos in the K inner product (Grimes, Lewis & Simon 1994):
without a shift on K⁻¹(−K_σ) (θ = 1/λ, Cholesky of K), the lowest positive λ; with `"shift": σ`
on the buckling transformation (K + σ K_σ)⁻¹ K (θ = λ/(λ − σ), LDLᵀ), the λ above σ, and the
LDLᵀ inertia counts the load factors between 0 and σ, reported. Negative load factors (the
reversed load buckles) are reported, not returned. Modes are scaled to a largest component of 1
and display and probe like frequency modes (their value is λ; `FrequencyResult::buckling`).
Boundary conditions hold their directions at zero in the eigenproblem; pressures are dead loads
(no load stiffness); trusses, shells and beams are refused (no K_σ yet). Decks: Abaqus
`*BUCKLE, EIGENSOLVER=LANCZOS` (refused after a general step, whose loaded state Abaqus would
buckle from, and with a shift); Nastran SOL 105 and LS-DYNA refuse it by name. Verified by Euler
columns of hex20 (pinned and clamped–free, 0.18 % from Engesser's shear-corrected load) and the
shifted form against the plain one.

**Imperfections.** A nonlinear static step's `"imperfection": {"step": "Buckle", "mode": 1,
"scale": 0.01}` solves that buckling step first and adds its mode, scaled to a largest component
of `scale`, to the node coordinates; the step is solved and reported on that geometry, with a
warning naming the mode. Abaqus and LS-DYNA decks refuse it by name (their imperfection reads a
results file). Verified: a clamped–free column with its first mode as imperfection under half its
buckling load deflects by δ₀ α/(1 − α) (2.8 %).

## Nonlinear static analysis (finite strain)

A step `{"type": "nonlinear_static", …}` (`Procedure::NonlinearStatic`, solved by
`solve.nonlinear`, `fieldcraft_solver::nonlinear`) follows equilibrium in load increments with
Newton–Raphson iterations. The step's load factor λ runs from 0 to 1: dead loads are λ times
their values, a pressure follows the deformed faces with λ times its magnitude, and prescribed
displacements are λ times theirs (a linear ramp).

**Material point** (`fieldcraft_fem::constitutive`). Given the deformation gradient F and the
history committed at the end of the last converged increment, a material returns the
**Kirchhoff stress** τ = J σ and the **spatial tangent** c of its Lie derivative (L_v τ = c : d),
in Voigt order with engineering-shear columns, so c = D at F = I. This is the push-forward of
the second Piola–Kirchhoff stress and its tangent: τ = F S Fᵀ, c_ijkl = F_iI F_jJ F_kK F_lL
C_IJKL. History variables (`MaterialModel::history_len` per integration point) are read from the
committed slice and written to a separate one, which the solver commits when the increment
converges, so a cut-back increment leaves no trace. The materials:

| Material in a nonlinear step | Model |
|---|---|
| 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 `{"type": "mooney_rivlin", "c01": …}` or `{"type": "yeoh", "c20": …, "c30": …}` | Decoupled polynomials W̄ = C₁₀ (Ī₁ − 3) + C₀₁ (Ī₂ − 3) + C₂₀ (Ī₁ − 3)² + C₃₀ (Ī₁ − 3)³ with the same volume κ/2 (J − 1)²: Mooney–Rivlin with C₁₀ = μ/2 − C₀₁ (0 ≤ C₀₁ ≤ μ/2), Yeoh with C₁₀ = μ/2, so 2(C₁₀ + C₀₁) is the elastic μ and the small-strain material stays the elastic one. Stress and tangent from the fictitious τ̄ and c̄ (Holzapfel 2000, eqs. 6.168, 6.193). Abaqus `*HYPERELASTIC, MOONEY-RIVLIN` / `YEOH`; LS-DYNA refuses them by name. |

| 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.

**Elements** (`fieldcraft_fem::finite`), total Lagrangian over the reference element with its
own quadrature: f_a = ∫ Bₓ,aᵀ τ dV₀ and K_ab = ∫ (Bₓ,aᵀ c Bₓ,b + ∇ₓN_a·τ∇ₓN_b I) dV₀ with the
current gradients ∇ₓN = F⁻ᵀ∇₀N (material plus geometric stiffness), identical term by term to
∫ B₀ᵀ S dV₀ and its linearisation. Every element's tangent is checked against central
differences of its forces.

| Element | Nonlinear steps |
|---|---|
| tet4, tet10, wedge6, hex8, hex20 (full), hex20 (reduced) | total Lagrangian |
| hybrid hex8, hex20, tet10 | Simo–Taylor–Pister three-field form on the linear hybrid elements' pressure fields (θ the projection of J, p the projection of U′(θ)); needs the neo-Hookean option (St Venant–Kirchhoff has no volumetric split) |
| hex8 and quad4 incompatible modes | refused: their finite-strain form (enhanced assumed strain, Simo–Armero 1992) needs the modes as unknowns carried between iterations and hourglasses in large compression; use full integration, hybrid or quadratic elements |
| plane strain, axisymmetric (tri3, tri6, quad4, quad8) | total Lagrangian with F₃₃ = 1 or r/R (hoop column N/r, geometric term τ_θθ N_a N_b/r²) |
| 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`) 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
K = p ∫ N_a (N_b,η [x,ξ]× − N_b,ξ [x,η]×). The load stiffness is not symmetric. The solver keeps
its symmetric part and leaves the skew part out, so the tangent stays symmetric for the
supernodal Cholesky factorization: the residual is exact (the converged answer does not
change), and over a closed surface or one whose boundary is held the skew parts of neighbouring
faces cancel (Schweizerhof & Ramm 1984), so nothing is lost there; elsewhere Newton may converge
linearly instead of quadratically, and the result says how large the left-out part is. An LU
factorization would double the factor's memory and lose the symmetric solver for a gain only in
the rate of convergence of open loaded surfaces; it is revisited with the symmetric indefinite
(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
displacement correction (over the increment's largest displacement) and the energy (over the
first iteration's) are all below their tolerances, or the residual is round-off. Each increment
starts from the linear extrapolation of the last one. Automatic incrementation grows an
increment by half after two fast ones (≤ 4 iterations), shrinks it after a slow one (> 8), and
cuts it to a quarter when an element inverts, the residual is not finite or grows twice, the
tangent is singular, or the iterations run out; the step stops with the load factor
and the reason at the minimum increment, the increment cap, or `Limits::max_newton_iterations`.

**Indefinite tangents.** A tangent that is not positive definite (a strained nearly
incompressible mixed element, an unstable material, a state past a limit or bifurcation point)
is factored by LDLᵀ and Newton goes on; each iteration records the tangent's negative
eigenvalues, and the summary warns when a converged state keeps some (an unstable equilibrium
under load control). This lifts the old limit on hybrid elements near ν = ½ under prescribed
displacement (at ν = 0.4999 every increment used to be cut back to the cap). Only the first
tangent of a model that can move freely is an error. The tangent predictor of path-dependent
materials keeps Cholesky and is skipped when the tangent is indefinite.

**Arc length.** `"arc_length": {"constraint": "cylindrical"|"spherical", "max_load_factor": 1,
"displacement_limit": {"node", "direction", "displacement"}}` makes the load factor an unknown
(Riks 1979, Crisfield 1981): each increment keeps ‖Δq‖² + w Δλ² = Δl² (w = 0 cylindrical, w the
first tangent displacement's ‖q_f‖² for the spherical arc), predicts along the tangent with the
sign of its angle to the last increment (Feng, Perić & Owen 1996), and corrects with
δq = δq_r + δλ δq_g, δλ the root of the constraint's quadratic whose increment turns least
(Crisfield's angle criterion). The first increment's load factor change is `initial_increment`,
which sets the arc; arcs grow, shrink and cut back as increments do, within the minimum and
maximum increments' arcs. Prescribed motions scale with λ (∂r/∂λ includes −K T_p u_p) and
follower pressures with their current faces. The tangent at every converged state is factored
for the next predictor (LDLᵀ when indefinite); its negative eigenvalues are recorded per
increment (`negative_pivots`), and where their count changes the step reports a **critical
point**: a limit point where the load factor turns, a bifurcation otherwise
(`History::critical_points`, in the `solve.nonlinear` summary). The step ends at the maximum load
factor, at the displacement limit, or at its increment cap (a warning). Refused: a load path,
temperatures. Decks: Abaqus `*STATIC, RIKS` (cylindrical only); LS-DYNA refuses arc-length steps
by name. Verified by the snap-through of a shallow sinusoidal arch of hex20 solids
(`crates/solver/tests/arc_length.rs`, `tests/oracle/arc-shallow-arch`): the path within 1.9 % of
the one-mode closed form, both limit points, and the inertia 0 → 1 → 0 exactly between them.

**Results**: the final displacements and reactions, the **Cauchy stress** and the
**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
reference points), with a warning; springs act along their fixed directions. Centrifugal loads
are refused by name in nonlinear steps.

**Thermal strain.** A nonlinear step's temperatures are its temperature loads, or the steady
conduction of its temperature boundary conditions and heat loads solved first (one-way, as the
linear coupled temperature–displacement step); they ramp with the load factor from each
material's reference temperature. With nonlinear geometry F = F_m F_θ, F_θ = diag(1 + αᵢ ΔT) at
each integration point (V. A. Lubarda, Appl. Mech. Rev. 57:95–108, 2004): the material sees F_m,
the loops run on the expanded volume J_θ dV₀ (mixed elements project det F_m), and the reported
strain is the total logarithmic one; the linear stretch 1 + α ΔT keeps α the engineering
coefficient, so a free body grows by exactly α ΔT. Without nonlinear geometry the split is
additive (ε_m = ε − α ΔT), and elastic elements on their linear stiffness take the equivalent
nodal forces as a dead load. Johnson–Cook plasticity softens at the point's temperature. The
expansion is constant in temperature (the model has no temperature-dependent coefficient yet).

**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`; 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; 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`,
and `tests/oracle/nl-*`): consistent tangents of both materials and of every element (central
differences); uniaxial stretch to λ = 2 of a neo-Hookean block in every element kind, plane
strain and axisymmetric (closed form); St Venant–Kirchhoff uniaxial stretch; the nonlinear patch
test on a distorted mesh (all solids, hybrid ones included); a 90° rigid rotation (zero stress);
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). 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. Mooney–Rivlin and Yeoh (`nl-mooney-rivlin-modes`, `nl-yeoh-uniaxial`,
`nl-mooney-rivlin-tube`): tangents by central differences (material point and hybrid
elements), Treloar's uniaxial, equibiaxial and pure-shear modes against the compressible closed
form (exact), and a tube inflated by a follower pressure while stretched axially to λ_z = 1.5
(Rivlin's incompressible solution, 0.6 % at ν = 0.4995). Thermal strain (`nl-thermal-expansion`,
`j2-heated-bar`, `crates/solver/tests/nonlinear_thermal.rs`): free expansion and a fully held
block at a thermal stretch of 1.2 (exact, full and hybrid), a bar held at its ends heated past
yield (exact), the linear limit against the linear and coupled solvers, and tangents with a
temperature field.

crushing stops at the minimum increment, naming it; the caps; what is refused).

## Contact and friction (explicit dynamics)

**Model** (`fieldcraft_model::contact`): a contact names its surfaces as sets of element faces —
a **pair** (`main`, `secondary`) or **general** contact on one surface (self contact) or on every
exterior face of the solids — with Coulomb **friction** μ = μ_k + (μ_s − μ_k) e^{−d γ̇} (static,
kinetic, decay per slip rate), a **penalty scale** s ≤ 1 and what happens to nodes that start
inside the other surface (**adjust**: moved onto it strain-free and reported; **refuse**: the
analysis stops naming the node, the face and the depth). Planar rigid walls take the same friction.
Commands `interaction.contact.pair|general|edit|delete|find`; Abaqus writes `*SURFACE
INTERACTION` with `*FRICTION` (exponential decay), `*CONTACT PAIR` and general `*CONTACT`, LS-DYNA
`*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE` and `_SINGLE_SURFACE` (`fs`, `fd`, `dc`, `sfsa`,
`sfsb`); Nastran and the linear implicit analyses refuse contact by name; nonlinear static steps
take contact pairs (next section).

**Geometry** (`fieldcraft_fem::contact`, shared with the implicit solver to come): the closest
point of a triangle (exact region tests) or of a bilinear quadrilateral (projected Newton with the
exact Hessian) to a point, with the signed gap along the face's normal and the lateral offset;
a bucket grid of face boxes whose cell is twice the mean box, grown until the grid and its entries
fit their caps, queried by point or by box with the candidates in face order.

**Kinematics** (`fieldcraft_explicit::contact`): segment to segment. Each secondary face and
each main face facing it (normals within 84° of opposite) are projected onto the secondary face's
plane, split into triangles and clipped against each other; on each overlap triangle a
three-point rule takes the gap g = (x_s − x_m)·n_m between the faces' closest points. Where
g < 0 (no deeper than half the main face's element), a penalty pressure p = ε(−g) acts on the
secondary face's nodes with weights N_i^s along n_m and on the main face's nodes with −N_j^m. For
a uniform pressure across flat faces both meshes receive exactly the consistent nodal forces,
so the contact patch test passes on non-matching meshes. Every contact runs both ways (a pair
with each surface as the secondary, single-surface contact with every face against every other
face that shares no node with it) at half the penalty each.

**Penalty**: ε = s γ min_i(m_i/A_i)/Δt² over the contact's nodes (m_i lumped mass, A_i the node's
area facing one way, Δt the stable increment, γ = 0.3). By Gershgorin's and Weyl's inequalities
the contact adds at most 2 s γ/Δt² to the square of any frequency, which the default increment
scale 0.9 leaves room for: contact never shortens the increment.

**Friction**: at the secondary nodes, by a regularised return map against the main face that
pressed each node most: the elastic predictor t* = P t_n − k_t Δu_T (k_t = ε A_i/2), stick when
|t*| ≤ μ f_n, otherwise slip back to the cone, dissipating μ f_n (|t*| − μ f_n)/k_t. On rigid
walls, a node's tangential velocity changes by at most μ times its normal velocity change.

**Energy**: the contact forces' work (penalty energy plus friction dissipated) is integrated with
the trapezoidal rule and enters the balance; friction dissipation (contact and walls) is
reported apart. **Search** every increment, in parallel by faces, forces added in face order: the
results are the same on any number of threads.

**Verification** (`crates/explicit/tests/contact.rs`, `tests/oracle/explicit-*`): two bars
colliding (contact for 2L/c, interface force ρcAv₀/2, velocities exchanged, kinetic energy kept to
0.01%, the contact energy back to 10⁻⁴ of it), a block sliding to a stop under constant and
decaying friction on a rigid wall and on a deformable plate (distance and time to 0.1–0.5%), the
contact patch test on non-matching meshes (σ_zz uniform to 0.1%), two spheres colliding against
Hertz (peak force +2.3%, duration −0.4% on a graded quarter-sphere mesh), a free body untouched
by general contact on all its faces, initial penetrations adjusted or refused, the face cap, and
bit-for-bit results on one and four threads.

## Contact in nonlinear static steps (mortar, augmented Lagrangian)

**Discretisation** (`fieldcraft_solver::contact`): mortar contact, one pass, the pair's
`secondary` surface the non-mortar side. Each secondary node A carries a pressure p_A and the
weighted gap g̃_A = Σ ∫ Φ_A (x_m − x_s)·n dA over the overlaps of its faces with the main faces
facing them, Φ_A the secondary face's own shape functions; ĝ_A = g̃_A / A_A with A_A = ∫ Φ_A over
the node's faces. The pressure field is interpolated like the displacements, so a uniform
pressure crosses a non-matching interface as the consistent nodal forces of both meshes: the
contact patch test. Linear faces (tri3, quad4: tet4, wedge6, hex8); general contact, quadratic
faces, 2D models and arc-length steps are refused by name.

**Segmentation** (`fieldcraft_fem::mortar`): both faces are projected along the secondary
face's unit normal onto its plane, the main polygon is clipped by the secondary one
(Sutherland–Hodgman; a projected quadrilateral that is not convex goes as its two triangles), and
the overlap's fan of triangles is integrated with Radon's seven-point degree-5 rule. At each point
both faces' coordinates come from inverting their projected maps (linear, or bilinear by the root
of its quadratic), so the main point is on the main surface along the normal, warped or not. A
vertex within 10⁻⁸ of a clip edge's length of its line is on it and never cut: the intersection
of two nearly coincident edges is ill-conditioned and its derivative swamped by round-off. The
kernel is written over a generic scalar: `f64` for values, forward-mode `Dual` for gradients and
`Hyper` (gradient and packed Hessian in the 24 corner coordinates) for the consistent tangent,
the exact linearisation of the segmentation, normal and projection (`fieldcraft_fem::dual`).

**Enforcement**: augmented Lagrangian. With the multipliers λ_A fixed, p_A = ⟨λ_A − ε_A ĝ_A⟩,
the derivative of Π = Σ A_A/(2ε_A) ⟨λ_A − ε_A ĝ_A⟩²: internal forces −p_A ∂g̃_A/∂x and the tangent
Σ_closed (ε_A/A_A) ∂g̃_A ⊗ ∂g̃_A − p_A ∂²g̃_A (the second term the contact's geometric stiffness).
Once Newton converges an increment, λ_A ← p_A (Uzawa) and the increment is solved again from its
solution, until every closed node's gap is within 10⁻⁵ of the thickness of the element under it
and no open node keeps a multiplier (at most 25 augmentations, then a warning). The penalty per
length is ε_A = s · 10 · 3Ḡ/h (Ḡ the mean shear modulus, h = V/A of the element under the face,
s the contact's `penalty_scale`): about ten times the local element stiffness, and 3Ḡ rather than
the constrained modulus keeps it so for nearly incompressible rubber.

**Active set and Newton**: the max() decides at every evaluation which nodes are closed
(semismooth Newton); an increment whose active set alternates between two states three times is
cut back (chattering). The tangent's pattern takes the blocks of the closed nodes' face pairs on
top of the mesh's. With contact, the force scale of the convergence test is that of the
converged states, the loads and the current state rather than the largest force of every trial
(a trial can press a body deep into another before the active set settles). The augmentations'
re-solutions are measured against the increment's start and its first correction; the next
increment's size follows the increment's own iterations.

**Small and finite sliding**: without nonlinear geometry the face pairs (those within a main
element's thickness) and their mortar integrals D_ai = ∫ Φ_a N_i, M_aj = ∫ Φ_a N_j are fixed in the
reference configuration, so g̃_a = n · (Σ M_aj x_j − Σ D_ai x_i) is linear: exact and quadratic-free.
With nonlinear geometry the pairs are searched (bucket grid of the main faces' boxes grown by
their depth limit) and integrated on the current configuration at every evaluation; a pair deeper
than half the main element's thickness has passed through and is left out.

**Initial overclosure**: nodes whose gap starts negative are reported; `adjust` takes the overlap
out of their weighted gap (an offset: strain-free, nodes are not moved), `refuse` stops naming
the node and the depth.

**Output**: per secondary node the pressure (its normal contact force over the area it carries,
∫ Φ_A over the overlaps, which is smoother than the multipliers at the contact's edge), the
multiplier, the gap ĝ and the status (open, closed); per contact the force on the secondary
surface, the normal force Σ p_A ∫Φ_A, the deepest penetration and the augmentations.
`solve.nonlinear` summarises each contact, `solve.probe` gives a secondary node's `contact`
entries, and the results display draws `contact_pressure`, `contact_gap` and `contact_status`.

**Verification** (`crates/solver/tests/contact.rs`, `tests/oracle/nl-contact-patch`,
`tests/oracle/nl-hertz-cylinder`): the contact patch test on non-matching hex8 and tet4 meshes
(uniform σ_zz, pressure equal to it, no lateral stress: to 10⁻⁹ with small sliding, 10⁻⁴ with
finite sliding, where the gaps the augmentation leaves enter through their geometric terms);
Hertz's line contact, plane strain (peak pressure within 0.3%, the pressure at two nodes within
3% and 6%, the vertical stress under the contact within 2%, the half-width bracketed); Hertz's
sphere on a flat, a quarter model (peak pressure within 4%, contact radius bracketed; heavy);
Newton's quadratic rate near the solution with finite sliding; initial overclosure adjusted and
refused; bit-for-bit results on one and four threads; the refusals.

## Hyperelastic materials in explicit dynamics

Hyperelastic materials (neo-Hookean, Mooney–Rivlin, Yeoh) use the shared finite-strain material
point (`fieldcraft_fem::constitutive`) in explicit dynamics too, in total form, which is objective
by construction (no stress rate): the one-point hexahedron takes F = I + ∂u/∂X, its mean over the
undeformed element, and the Cauchy stress of the nodal forces τ/J; the tetrahedron takes the
isochoric stress of its own F and the **average nodal pressure** of a decoupled material (each
node's U′(J_a) from its tetrahedra's volume ratio, the element's the mean of its nodes': J. Bonet,
A. J. Burton, Commun. Numer. Methods Eng. 14:437–449, 1998, at finite strain), which does not lock
when the material is nearly incompressible. (Averaging J into an F-bar instead was unstable on
the Kuhn mesh of the tests.) The stable increment follows the tangent: each element's frequency
bound is scaled by the ratio of its tangent's largest normal stiffness, plus the stress, to the
small-strain one. Logarithmic strain is the strain output. Hyperelastic elements need nonlinear
geometry and take no thermal strain yet (both refused by name). Verified by uniaxial strain to
λ = 1.5 of all three laws on hex8 and tet4 against the material point (`crates/explicit/tests/hyperelastic.rs`)
and by `tests/oracle/explicit-neo-hookean-implicit`, a quasi-static explicit pull against the
implicit nonlinear solver.

## Mass scaling (explicit dynamics)

Opt-in per explicit step (`mass_scaling: {target_increment, set}`): every solid element (of the
set, or of the model) whose stable increment times the step's scale is below the target has its
density raised by (target / increment)², which brings its increment to the target (its lumped
masses, frequency bound and hourglass and bulk-viscosity scales follow; adiabatic heating keeps
the physical density). The run reports the elements scaled, the added mass in all and per
section, its fraction of the model's mass, a warning, and the kinetic energy of the added mass
(`energy.added_kinetic`, part of the kinetic energy, so the balance holds). Verified on a uniform
bar scaled to twice its increment: added mass three times the bar's, the increment doubled, the
added mass's share of the kinetic energy three quarters (`crates/explicit/tests/one_d.rs`).