Mandel transform utilities (sym tensor ↔ vector / rank-4 ↔ 6×6) for vector Newton over tensor-valued state

[petlenz]

Context

Under the numsim-codegen graph-coupled architecture, codegen emits constitutive materials that expose rate/residual + jacobian (and tangent blocks) as properties, and the numsim-materials solver layer drives the Newton iteration. The existing solvers (solvers/backward_euler.h, solvers/rk_integrator.h) are scalar-only — rk_integrator's Eigen system is the RK stage coupling for a single scalar ODE (s×s, m_J = I − h·diag(df)·A), not coupling across physical state variables.

A generated material whose local unknown is a symmetric second-rank tensor (e.g. a tensor-valued plastic strain, or coupled tensor internals) needs a vector Newton solver that flattens the symmetric tensor ↔ a vector and the rank-4 Jacobian ↔ a dense matrix. That flattening must use Mandel (not engineering Voigt) notation, and no such utility exists.

Confirmed absent across all branches (main, feature/{j2-plasticity, rk-integrators, drucker-prager, damage-materials, numerical-diff-checker}): no mandel/voigt symbol anywhere.

Why Mandel, not Voigt

The dense Newton solve J·Δx = −R must be numerically equivalent to the tensor system ∂R/∂x : Δx = −R. Mandel's √2 scaling on the off-diagonal/shear components makes the map an isometry:

• A : B == mandel(A) · mandel(B) (double contraction → vector dot product)
• (C :: X) (rank-4 on rank-2) == mandel4(C) · mandel2(X) (matrix–vector)

so inv/solve and norms (the convergence test) are preserved. Engineering Voigt does not preserve these without asymmetric factors, which corrupts the Jacobian solve and the consistent tangent.

Proposal

A header-only Mandel transform utility (3-D; extendable):

• mandel2(sym tensor<3,2>) → Eigen::Vector<double,6> and inverse unmandel2
• mandel4(sym tensor<3,4>) → Eigen::Matrix<double,6,6> and inverse unmandel4
• √2 on the 3 shear/off-diagonal slots (rank-2), 2 on the corresponding rank-4 cross blocks
• Works on tmech tensors on the codegen side and Eigen on the solve side (the two type systems meet here)

Tests

• Round-trip: unmandel2(mandel2(A)) == A, unmandel4(mandel4(C)) == C (to roundoff).
• Inner-product isometry: A:B == mandel2(A)·mandel2(B) over random symmetric A,B.
• Operator equivalence: mandel2(C :: X) == mandel4(C) · mandel2(X).
• Solve equivalence: for a symmetric rank-4 C and rank-2 R, unmandel2(mandel4(C).lu().solve(mandel2(R))) equals the tensor solution of C :: X = R.

Related / scope

• Prerequisite for a generic vector/tensor-valued local Newton solver — the existing backward_euler (scalar residual/jacobian) and rk_integrator (scalar rate/rate_derivative; Eigen is RK-stage, not state coupling) are both scalar-only. That solver is a separate, larger piece; this issue is the self-contained transform it depends on.
• Driven by the numsim-codegen graph-coupled roadmap (codegen emits the tensor-valued constitutive properties full-tensor; the solver does the Mandel flattening — codegen stays Voigt/Mandel-free).