Add Drucker-Prager plasticity with apex return

[petlenz]

Summary

Adds Drucker-Prager (DP) plasticity with non-associative flow and apex return
to the existing J2 framework. The same small_strain_plasticity template now
handles both J2 and DP via a yield-function policy. Also fixes a class of
tangent-consistency bugs and adds extensive documentation.

What's new

Drucker-Prager yield function

• drucker_prager_yield_function<T, Dim> policy with friction η, dilatancy
  β, bulk modulus K_bulk
• Non-associative flow rule (η ≠ β) — yield normal M = dF/dσ and flow
  normal N = dG/dσ are distinct
• Stateful policy passed via parameter handler at construction
• Type alias: drucker_prager_plasticity<Traits> mirrors j2_plasticity<Traits>

Apex return

• Detects when standard return would push √J₂ below zero (pressure-dominated
  loading)
• Pre-checks the apex condition from trial quantities to avoid wasted smooth
  Newton iterations
• Dedicated apex residual: r = η·p_trial - K·η·β·Δκ - k - H = 0
• Rank-1 volumetric apex tangent: K·H' / (K·η·β + H') · I⊗I
• Compile-time dispatch via has_apex_return concept — J2 path is unaffected

Tangent consistency fixes

• Original DP tangent had a constant ~1.8% error from a pressure-term
  mismatch between modified_equivalent_stress and yield_normal
• Root cause: √J₂ + η·I₁ paired with s/(2q) + η/3·I (factor-of-3 mismatch
  in the volumetric coupling)
• Fixed by unifying on √J₂ + η·p convention everywhere
• Tangent now uses trial-state quantities (matches the algorithm: ε_p update
  uses N_trial)
• DP tangent error now ~10⁻¹⁰ (machine precision); J2 tangent improved from
  ~10⁻⁴ to ~10⁻¹⁰ as a side effect

Refactoring

• Extracted shared helpers in plasticity_utils.h: evaluate_at_state,
  compute_trial, compute_tangent, make_IIdev
• small_strain_plasticity::compute() is now a flat dispatch: elastic /
  apex / smooth — each path is a named private member function
• solve_scalar_return(phi, G_eff, kappa_n) is shared between smooth and
  apex Newton (matches the unified residual structure)
• Renamed hardening variable α → κ (avoids clash with DP friction
  coefficient; matches the document notation)
• DP friction coefficient renamed α → η in code

Documentation

• New docs/small_strain_plasticity.md: full derivation of yield functions,
  return-mapping algorithm, consistent tangent (smooth and apex), policy
  interface, consistency requirements, and common failure modes
• Includes the J2-to-DP normalization relationship and stress sign convention

Tests and tooling

• test_drucker_prager.cpp: 6 tests covering elastic-before-yield, yielding,
  volumetric plastic strain, pressure sensitivity, consistent tangent, and
  tangent convergence (max error 5.2e-10)
• scripts/plot_plasticity.py + tests/plot_data.cpp: stress-strain and
  tangent-accuracy plots for three loading paths (ε₁₁, ε₂₂, ε₁₂)
• tests/debug_apex.cpp: standalone diagnostic for the apex-overshoot
  condition

Bug fixes (also applied to existing code)

• rk_plasticity.h: jacobian(m_G, ...) should be
  jacobian(effective_modulus(m_G), ...) — pre-existing scaling bug
• backward_euler::solve(): now reports convergence via converged() and
  clamps result to ≥ 0 (no backward plastic flow)
• small_strain_plasticity::compute(): throws on non-convergence (was silent
  wrong-result bug); falls back to apex if smooth Newton fails and apex
  support exists
• apex_tangent: scale-relative zero-denominator threshold
  (eps · (|Kηβ| + |H'|)) instead of numeric_limits::min()
• flow_normal_stress_derivative now takes sig_dev directly — avoids
  cancellation error from reconstructing s from non-associative N
• evaluate_at_state threshold scaled to σ_0 — prevents subnormal-stress
  divisions in downstream 1/J₂ terms

Test plan

• All 38 existing tests pass (5 suites: J2, DP, RK, materials, damage)
• DP consistent-tangent test passes at machine precision (~10⁻¹⁰)
• DP convergence test: max tangent error 5.2e-10 (was 1.84e-02 before
  fixes)
• J2 tests still pass — no regression from the unified policy interface
• Apex return verified on uniaxial-strain loading (oscillating-stress
  bug fixed)
• Three loading directions (ε₁₁, ε₂₂, ε₁₂) tested and plotted

[petlenz]

Review from the numsim-codegen side

Reviewing with the codegen hat on: numsim-codegen's NumSimMaterialTarget generates constitutive materials that target exactly this runtime's solver/property contract (<function>::rate/rate_derivative and residual/jacobian over Local edges, driven by update_source()). So this solver layer is shared infrastructure and changes ripple into generated code. Overall this is careful, well-structured plasticity — findings are mostly about the generic solver picking up domain-specific behavior, plus one tangent-consistency question.

Done well

• Scalar return mapping in full-tensor representation. The local unknown is the scalar multiplier Δλ; the tensor update is closed-form (ε_p += Δλ·N). That's why a scalar solver suffices and there is no Voigt/Mandel flattening. (So the Mandel + vector-solver gaps I filed — #11/#12 — are for a different class, coupled tensor-valued local systems, not this. This PR correctly sidesteps them.)
• Convergence-failure handling is robust in small_strain_plasticity::compute() — checks converged(), falls back to apex, and throws if both fail rather than using a bad iterate. Better rigor than codegen's current in-function Newton (numsim-codegen#85).
• The conservative apex pre-check (dλ_max = F_trial/G_eff) to skip a doomed smooth Newton is a nice optimization.

Findings

HIGH — backward_euler is a generic name with domain-specific behavior baked in. (solvers/backward_euler.h)

• update() ends with m_delta = std::abs(m_delta); — "curing degree can only increase" (:80). A curing-specific monotonicity assumption hardcoded into a generic solver.
• solve() clamps every result with std::max(x, value_type{0}) — "negative plastic-multiplier unphysical" (:93,94,98). A plasticity-specific assumption, also in the generic solver.

So "backward_euler" is really two domain-specific solvers wearing a generic name. This blocks codegen reuse: a generated material whose state can legitimately decrease or go negative would be silently corrupted by abs()/max(·,0). Suggest hoisting the sign/clamp policy out of the solver (a clamp_nonnegative parameter defaulting off, or keep the clamp in the plasticity caller).

MEDIUM — update() and solve() disagree on convergence reporting. solve() sets m_converged; update() never does — on max_iter exhaustion it silently proceeds with the last iterate. The graph-driven path (curing, and any codegen rate material wired through it) has no non-convergence signal. Mirror converged() into update().

MEDIUM — Drucker-Prager algorithmic tangent uses the trial flow direction. materials/small_strain_plasticity.h:168-169. For DP the converged flow direction differs from trial, so N_trial gives a tangent that is not the consistent algorithmic tangent — degrading host-FE Newton from quadratic toward linear near yielding. The material advertises a tangent output; is the trial-direction tangent an intentional documented approximation, or should DP assemble at the converged direction?

LOW — magic constants in the solver. The m_delta = 5e-12 seed (:59), the 1e-30 singular-Jacobian threshold (:94), and the 5-iteration NaN-damping cap should be named/parameterized.

Cross-repo note

The residual/jacobian (and rate/rate_derivative) Local-edge contract itself is clean and is what codegen emits against. The one thing currently preventing a generated material from reusing backward_euler as-is is the HIGH finding (the baked-in sign policy). Hoist that out and this solver becomes directly consumable by codegen output — the goal of the graph-coupled architecture.

[petlenz]

HIGH (codegen-relevant): std::max(x, 0) bakes a plasticity assumption (Δλ ≥ 0) into a generically-named solver. A codegen-generated material solving a general residual whose root is legitimately negative would be silently clamped to 0. Suggest moving the clamp into the plasticity caller, or gating it behind a clamp_nonnegative parameter (default off).

[petlenz]

The trial-direction tangent is not the consistent algorithmic tangent for DP (N differs trial vs converged), which degrades host-FE Newton convergence near yielding. Intentional documented approximation, or should the DP tangent be assembled at the converged flow direction? The material advertises a tangent output, so worth pinning down.