Hands-on: writing a UMAT with simcoon

This tutorial builds a complete constitutive subroutine (UMAT) from scratch: J2 plasticity with linear isotropic hardening, written with the typed Tensor2/Tensor4 API so the code reads like the equations. The model is chosen because everything has a closed form — the return mapping needs no local iteration and the continuum tangent can be derived by hand, so every line of the implementation can be checked analytically.

The complete companion code lives in examples/umat_tutorial/ (umat_tutorial_J2.hpp/.cpp) and is continuously validated by the Ttutorial_umat test against the in-tree EPICP reference, the shared tangent assembly, and a finite-difference Jacobian.

The model

Additive small-strain decomposition, linear isotropic elasticity, von Mises yield with a linear hardening law:

\[\boldsymbol{\sigma} = \mathbf{L} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{th} - \boldsymbol{\varepsilon}^p), \qquad f = \sigma_{eq} - (\sigma_Y + H\,p) \le 0, \qquad \sigma_{eq} = \sqrt{\tfrac{3}{2}\,\mathbf{s}:\mathbf{s}},\]

with associated flow \(\dot{\boldsymbol{\varepsilon}}^p = \dot{p}\,\boldsymbol{\Lambda}\), \(\boldsymbol{\Lambda} = \tfrac{3}{2}\,\mathbf{s}/\sigma_{eq}\).

Props (5): E, nu, alpha, sigma_Y, H. State (8): T_init, p, EP(6) — deliberately the EPICP layout, so the two are directly comparable (EPICP with m = 1 and k = H is the same model).

Step 1 — anatomy of a UMAT

Every simcoon UMAT has the same signature (see any kernel or the tutorial header): strain state in (Etot, DEtot, Voigt order [11, 22, 33, 12, 13, 23], MPa), stress and tangent out (sigma, Lt), material properties (props), persistent state (statev), the rotation increment DR for objectivity, cumulative work outputs (Wm, Wm_r, Wm_ir, Wm_d) and the tangent_mode selector.

Three contract points that every UMAT must honor:

start: on the first increment, initialize the state and RESET the cumulative work accumulators.
objectivity: stored tensorial state co-rotates with DR. With the typed API this is one line — EP.rotate(...) picks the strain rotation kernel (factor-2 shear convention) from the tensor’s own type.
cumulative work: Wm += increment — the solver reports path integrals, never overwrite them.

Step 2 — trial state, typed

    // ------------------------------------------------------------------ 4.
    // Elastic prediction (total form): the trial state freezes plastic flow.
    //   eps_el^trial = Etot + DEtot - eps_th - EP
    //   sigma^trial  = L : eps_el^trial
    const tensor2 eps_th =
        strain(vec(alpha * (T + DT - T_init) * Ith()));
    const tensor2 eps_el_trial =
        strain(vec(Etot + DEtot)) - eps_th - EP;
    const tensor2 sig_trial = L_el.contract(eps_el_trial);  // stress-typed

Because L_el is a stiffness-typed tensor4 and the operand is a strain-typed tensor2, the contraction produces a stress-typed result — the shear-convention bookkeeping (engineering factor 2) is carried by the types, not by the programmer.

Step 3 — radial return, closed form

The trial direction is preserved by the return (the correction is radial in deviatoric space), and with linear hardening the consistency condition is linear in \(\Delta p\):

\[f(\Delta p) = f^{trial} - (3\mu + H)\,\Delta p = 0 \quad\Longrightarrow\quad \Delta p = \frac{f^{trial}}{3\mu + H}.\]

    if (f_trial > 0.) {
        // -------------------------------------------------------------- 6.
        // RADIAL RETURN, closed form. The flow direction
        //   Lambda = (3/2) dev(sigma)/sigma_eq
        // is preserved by the return (the correction is radial in deviatoric
        // space), and with LINEAR hardening the consistency condition
        //   f(Dp) = f^trial - (3 mu + H) Dp = 0
        // is linear in Dp — no local Newton iteration is needed:
        Lam = flow_normal(sig_trial);            // strain-typed
        Dp = f_trial / (3. * mu + H);

        p += Dp;
        EP += Dp * Lam;
    }

    // Converged stress from the total form (elastic when Dp = 0).
    const tensor2 sig = L_el.contract(eps_el_trial - Dp * Lam);
    sigma = sig.to_arma_voigt();

Step 4 — the continuum tangent, derived and built

Differentiating the converged stress w.r.t. the total strain with the flow direction frozen gives the classical rank-one update

\[\mathbf{L}_t = \mathbf{L} - \frac{(\mathbf{L}:\boldsymbol{\Lambda}) \otimes (\mathbf{L}:\boldsymbol{\Lambda})} {\boldsymbol{\Lambda}:\mathbf{L}:\boldsymbol{\Lambda} + H}.\]

For isotropic \(\mathbf{L}\) and the deviatoric von Mises direction the pieces are analytic — \(\mathbf{L}:\boldsymbol{\Lambda} = 2\mu\boldsymbol{\Lambda}\) and \(\boldsymbol{\Lambda}:\mathbf{L}: \boldsymbol{\Lambda} = 3\mu\) — hence the closed form

\[\mathbf{L}_t = \mathbf{L} - \frac{4\mu^2}{3\mu + H}\, \boldsymbol{\Lambda}\otimes\boldsymbol{\Lambda}.\]

    if (Dp <= simcoon::iota || tangent_mode == tangent_none) {
        Lt = L;  // elastic step, or explicit integration
    } else if (tangent_mode == tangent_continuum) {
        const tensor2 kappa = L_el.contract(Lam);            // stress-typed
        const double Lam_L_Lam = accu(Lam.mat() % kappa.mat());  // = 3 mu
        const tensor4 Lt_t =
            L_el - auto_dyadic(kappa) * (1. / (Lam_L_Lam + H));
        Lt = mat(Lt_t.mat());
    } else {
        // Algorithmic (Simo-Hughes) and beyond: hand the physical inputs to
        // the shared dispatch — Bhat = Lambda:kappa + H, the flow Hessian is
        // dLambda/dsigma = deta_stress. This is what the production kernels do.
        const tensor2 kappa = L_el.contract(Lam);
        const double Bhat = accu(Lam.mat() % kappa.mat()) + H;
        const vec sig_v = sigma;
        const ContinuumTangent ct = compute_tangent_operator(
            tangent_mode, Bhat, kappa.to_arma_voigt(), Lam.to_arma_voigt(),
            Dp, L, [&sig_v]() { return deta_stress(sig_v); });
        Lt = ct.Lt;
    }

Warning

The classic Voigt trap. The closed form above holds for the tensorial \(\boldsymbol{\Lambda}\). In Voigt components the outer product must be taken with the stress-Voigt normal (shear not doubled): \(\boldsymbol{\kappa} = \mathbf{L}:\boldsymbol{\Lambda}\) lands in stress Voigt automatically, but eta_stress() returns the engineering (doubled-shear) normal — using it directly in \(\boldsymbol{\Lambda}\boldsymbol{\Lambda}^T\) overweights the shear block of the correction by up to a factor 4. The typed auto_dyadic route sidesteps the trap (the type carries the convention); note also that sym_dyadic is the \((ik)(jl)\)-symmetrized product — a different operator from the plain dyadic needed here.

For the algorithmic (Simo–Hughes) operator the tutorial does what the production kernels do — hand the physical inputs to the shared dispatch compute_tangent_operator with a lazy flow-Hessian provider; for J2 the resulting operator is the exact Jacobian of the discrete radial-return map.

Step 5 — validation

test/Libraries/Umat/Ttutorial_umat.cpp asserts, on every test run:

Physics: the tutorial matches the EPICP reference to the CCP tolerance (< 1e-8 relative) over a cyclic tension/shear path.
Tangent: the hand-built continuum operator equals both the independent closed form (with the stress-Voigt normal!) and assemble_continuum_tangent to machine precision.
Consistency: the algorithmic operator equals the central finite-difference Jacobian of the discrete stress update (< 1e-5).

Where to go next

Swap the yield criterion: replace flow_normal by a Hill or DFA gradient (criteria.hpp) — the generic \(\kappa\)-form of the tangent in the code keeps working; only the closed-form shortcut is J2-specific.
Add kinematic hardening: shift the criterion to \(\boldsymbol{\sigma} - \mathbf{X}\) and consult the kept EPCHA kernel for the classic Armstrong–Frederick treatment.
Or skip hand-written kernels entirely: the same model is one line of simcoon.modular.ModularMaterial — see Constitutive model (UMAT) catalog.