Constitutive laws in Python (PYEXT)

A constitutive law can be written in Python and driven by the C++ material-point solver exactly like a built-in kernel. The C++ side keeps the Newton loop, the step cutting, the start-of-increment rollback and the finite-strain kinematics; the Python side only integrates one increment at one material point. The law is served under the UMAT name PYEXT and reaches the C++ code through a process-wide callback registered by the bindings.

Typical uses: a research model prototyped in numpy, a machine-learning model (PyTorch, JAX), or any law for which a C++ port is not worth the effort yet.

Quick start

import numpy as np
import simcoon as sim
from simcoon.solver import StepMeca, solve

class LinearElastic(sim.PythonUMAT):
    nstatev = 0                      # no internal variable

    def __init__(self, E, nu):
        self.L = sim.L_iso([E, nu], "Enu")

    def integrate(self, *, Etot, DEtot, sigma, statev, Wm, **kw):
        stress = self.L @ (Etot + DEtot)
        Wm = Wm.copy()
        Wm[0] += 0.5 * (sigma + stress) @ DEtot     # cumulative work
        return stress, self.L, statev, Wm            # sigma, Lt, statev, Wm

step = StepMeca(control=["strain"] + ["stress"] * 5, value=[0.01, 0, 0, 0, 0, 0], ninc=50)
res = solve(step, LinearElastic(70000., 0.3))       # the law object replaces the name
res["Stress"][0, -1]                                 # 700.0

solve registers the object under PYEXT for the duration of the call; props and nstatev are taken from the object (props may stay empty).

The contract

simcoon.PythonUMAT mirrors the C++ small-strain UMAT convention. The integrate method receives keyword arguments and returns a tuple:

def integrate(self, *, Etot, DEtot, sigma, DR, props, statev, T, DT, Time, DTime,
              Wm, ndi, nshr, start, tangent_mode, **_):
    ...
    return sigma, Lt, statev, Wm            # optionally a 5th value: L (elastic operator)

Voigt convention: order 11 22 33 12 13 23, engineering shear strains, quantities in the material (local) frame. Etot is the strain at the beginning of the increment, DEtot the increment, sigma the stress at the beginning of the increment; DR is the rotation increment (identity in small strain).
State: the law must be a pure function of its arguments except through statev and Wm. The solver calls the same increment several times (Newton iterations, step cuts) with statev, sigma and Wm reset to their start-of-increment values — nothing may be cached in the Python object between calls. A recurrent model keeps its hidden state inside statev.
start is True on the first call of a block (Time == 0); statev arrives zero-filled — initialise it there if needed. The solver primes the tangent at the start of every block with a zero-increment call (DTime == 0, DEtot == 0), and finite-element couplers do the same at initialisation; simcoon.PythonUMAT treats that probe as a pure tangent query: stress, statev and Wm are returned exactly as received and only the tangent comes from the evaluation, so a history-dependent law does not count it as a loading step (a bare callable registered instead of a PythonUMAT must do it itself).
ndi follows the classical convention: 3 = 3D (plane strain is 3D with a null out-of-plane strain), 2 = plane stress (the law must condense, cf. el_pred), 1 = uniaxial. nshr is the number of shear components.
tangent_mode: 0 = return the elastic operator as Lt, 1 = continuum tangent, 2 = algorithmic (consistent) tangent (default of the solver).
Wm is the accumulated (Wm, Wm_r, Wm_ir, Wm_d) at the beginning of the increment and must be returned accumulated.
Return: sigma (6,), Lt (6,6), statev (nstatev,), Wm (4,) and optionally L (6,6) (defaults to Lt). Lists, float32 arrays and CPU torch tensors are converted; shapes are checked and a non-finite value is an error.
Finite strain: under NLGEOM the caller feeds the corotational logarithmic strain and expects the Kirchhoff stress — the same convention as ELISO / EPICP (PYEXT belongs to the Kirchhoff-box set). Internal tensorial history is not rotated by the solver; rotate it with DR in the law if needed.

Step cuts and errors

Raise simcoon.StepCut to ask the solver for a smaller increment (the trial is discarded and the increment retried; the effective factor is the solver’s div_tnew_dt). Any other exception aborts the solve and is re-raised unchanged in Python (type and traceback preserved), including KeyboardInterrupt. With inforce=0 the solver aborts (status = 1, a RuntimeError unless raise_on_abort=False) when the increment falls below Dn_mini; with the default inforce=1 it forces the minimal increment.

Batch entry point and explicit registration

sim.umat("PYEXT", ...) (the per-Gauss-point batch call used by finite-element couplers) works with a registered law:

with sim.registered(law):
    stress, statev, Wm, Lt = sim.umat("PYEXT", etot, Detot, F0, F1, sigma, DR,
                                      props, statev, time, dtime, Wm, ndi=3)

The points are integrated serially on the calling thread (the callback re-enters the interpreter, so it never runs inside the parallel region; n_threads is ignored). simcoon.registered() restores the previously registered law on exit; simcoon.pyumat.register() / simcoon.pyumat.unregister() are the explicit forms.

Performance

Per call, the bridge acquires the GIL (released by the solver around the C++ loop), builds small numpy copies of the inputs and copies the outputs back — about 5–10 µs, negligible against any non-trivial law. A numpy J2 law costs a few tens of µs per call, a small LSTM step with its autograd tangent about a millisecond. For finite-element scale, batch the evaluation on the Python side (a batched step of its own) rather than calling sim.umat point by point.

Examples

J2 plasticity written in Python (PYEXT) — J2 plasticity in numpy compared with EPICP.

API reference

class simcoon.PythonUMAT

Point-wise small-strain constitutive law implemented in Python.

Subclasses set nstatev (number of internal variables), optionally props (material properties forwarded to the C++ material record; may stay empty since the object holds its own parameters) and implement integrate().

Contract of integrate() (identical to the C++ kernels):

stateless between calls except through statev and Wm: the solver re-calls the same increment several times (Newton iterations, step cuts) with statev/sigma/Wm reset to their start-of-increment values. A recurrent model must keep its hidden state inside statev.
start is True on the first call of a block (Time == 0): initialise statev/Wm there. statev arrives zero-filled.
the solver primes the tangent at the start of every block with a zero-increment call (DTime, DT and DEtot all zero), and finite-element couplers do the same at initialisation. That probe is not a step of the loading history: PythonUMAT answers it as a pure tangent query (see __call__()), so integrate need not care. A bare callable registered instead of a PythonUMAT must handle it itself.
Voigt order 11 22 33 12 13 23, engineering shear strains, quantities in the material (local) frame. Etot is the strain at the beginning of the increment, DEtot the increment, sigma the stress at the beginning of the increment. Under finite strain the caller feeds the logarithmic strain and expects the Kirchhoff stress (same convention as ELISO/EPICP).
ndi follows the classical convention: 3 = 3D (plane strain is 3D with a null out-of-plane strain), 2 = plane stress, 1 = uniaxial. nshr = number of shear components.
tangent_mode: 0 = return the elastic operator as Lt, 1 = continuum tangent, 2 = algorithmic (consistent) tangent.
Wm is the accumulated (Wm, Wm_r, Wm_ir, Wm_d) at the beginning of the increment and must be returned accumulated.
returns (sigma (6,), Lt (6,6), statev (nstatev,), Wm (4,)[, L (6,6)]) as float arrays (lists / float32 / CPU torch tensors are converted). L defaults to Lt. Raise StepCut to request a smaller increment; any other exception aborts the solve and is re-raised unchanged.

abstractmethod integrate(*, Etot: ndarray, DEtot: ndarray, sigma: ndarray, DR: ndarray, props: ndarray, statev: ndarray, T: float, DT: float, Time: float, DTime: float, Wm: ndarray, ndi: int, nshr: int, start: bool, tangent_mode: int, **_): Integrate one increment at one material point (see class docstring).

nstatev: int = 0: Number of internal state variables.

props: ndarray = array([], dtype=float64): Material properties forwarded to the C++ material record (floats). May be empty. The default is a shared read-only array: a law with properties assigns its own (self.props = np.array([...])) rather than mutating this one in place.

class simcoon.StepCut(ratio: float = 0.5, msg: str = 'step cut requested')

Raise from PythonUMAT.integrate() to ask the solver for a smaller increment.

The current trial is discarded (the solver restores the start-of-increment state) and the increment is retried with a smaller time step.

Parameters:

ratio (float) – Suggested reduction factor in (0, 1), forwarded as tnew_dt. The material-point solver applies it directly when the cut interrupts a strain-driven increment and replaces it by its own div_tnew_dt when a Newton loop had to be abandoned.
msg (str) – Message attached to the exception.

simcoon.registered(umat)

Context manager: register umat for the block, restore the previous law on exit.

Parameters:: umat (PythonUMAT or callable) – The law to serve under the PYEXT name.

simcoon.pyumat.register(umat) → None

Register umat as the process-wide PYEXT law (prefer registered()).

The slot is process-wide: registering from a second thread while another thread’s law is in place raises instead of silently rebinding that thread’s running solve.

simcoon.pyumat.unregister() → None: Remove the process-wide PYEXT law (only the registering thread may do so while it is alive: a solve running on that thread must not lose its law).

simcoon.pyumat.current(): Return the currently registered Python law (or None).