From 85eaad027377bc6c1e082c7fd79720a995df12f9 Mon Sep 17 00:00:00 2001 From: petlenz Date: Mon, 14 Sep 2026 13:22:37 +0200 Subject: [PATCH 1/3] umat: add a CalculiX external-behaviour interface for the native umat_user hook CalculiX has two external hooks and they do not share a convention. call_external_umat is reached from umat_abaqus.f after that wrapper has already converted to the Abaqus UMAT convention, so NUMSIM_MATERIALS_DEFINE_UMAT already serves it. call_external_umat_user is the native hook, and this adapter translates its conventions onto umat_dispatch: - kode carries the constant count as -100 - nconst; - emec/emec0 are TENSORIAL, so the shear is doubled to the engineering form the dispatch buffer expects, exactly as umat_abaqus.f:280-283 does; - xstateini/xstate are the FULL (nstate_, mi(1), #elements) arrays, not a slice, so state is indexed at nstatv*((iint-1) + mi1*(iel-1)); - TIME is rebased onto the start of the increment per umat_abaqus.f:187-188; - the 6x6 tangent is symmetrized into stiff(21) as umat_abaqus.f:335-355 does, reading the column-major buffer narrow_matrix writes. A local orientation and a nonzero initial stress are refused rather than silently ignored: the native hook makes frame rotation the material's own responsibility, so returning unrotated results would be wrong with no symptom. Scope is geometrically linear; the PK2/Green-Lagrange pairing is exact only for a linear C:E law. REVIEW.md records the three-lens review this came out of, including the three silent-wrong-answer bugs it found. Each has a regression test verified to fail against the pre-fix code. --- REVIEW.md | 228 ++++++ .../umat/calculix_interface.h | 321 ++++++++ tests/CMakeLists.txt | 1 + tests/test_calculix_interface.cpp | 735 ++++++++++++++++++ 4 files changed, 1285 insertions(+) create mode 100644 REVIEW.md create mode 100644 include/numsim-materials/umat/calculix_interface.h create mode 100644 tests/test_calculix_interface.cpp diff --git a/REVIEW.md b/REVIEW.md new file mode 100644 index 0000000..ab65a84 --- /dev/null +++ b/REVIEW.md @@ -0,0 +1,228 @@ +# Review — CalculiX external-behaviour adapter + +> **STATUS: all findings addressed.** See "Resolution" at the end. The three +> bugs that could produce a silently wrong result (C1, H1, H2) each now have a +> regression test that was *verified to fail* against the old code before the +> fix was restored. + + +Scope: `include/numsim-materials/umat/calculix_interface.h` and +`tests/test_calculix_interface.cpp` (the native `call_external_umat_user` hook). + +Method: 3-lens parallel review (architecture / C++-ABI / code-quality & +coverage). Every Critical and High finding below was independently +re-verified by me against the CalculiX 2.22 source before landing here — one +agent's "packing is correct" was **wrong** (see H1), so the severities are mine, +not a straight merge of the agents'. + +The convention translation is correct **for the one configuration the single +test exercises** (`nstatv==0`, `iel=iint=1`, `iorien=0`, `time=0`, `emec0=0`, +one increment, symmetric tangent). Every finding below is a case that leaves +that configuration — and all are currently masked by it. + +--- + +## Critical + +### C1 — STATEV is addressed at point (1,1) for *every* integration point +`calculix_interface.h:90-91,103` (macro drops `IEL`/`IINT`/`MI`, `:149,156`) + +CalculiX passes the **whole** state arrays' base to the external hook, not a +per-point slice. Verified: +- `umat_main.f:40` — `real*8 xstate(nstate_,mi(1),*), xstateini(nstate_,mi(1),*)` +- `umat_main.f:206` (native `umat_user`) and the `@` branch at `:233` both pass + bare `xstateini,xstate` — **no** `(1,iint,iel)` slice, unlike + `umat_abaqus.f:295` which slices `xstate(1,iint,iel)` before its call. + +So the external `.so` must index state itself. The adapter uses offset 0 +(`std::copy_n(statev_old, …, statev_new)` and `a.statev = statev_new`), so every +Gauss point of every element reads and writes the **same** state slot; all other +points are never seeded or updated. Silent wrong physics for any model with +`nstatv > 0`. Invisible today only because linear-elastic has `nstatv == 0` +(the copy is a no-op) and the test uses `iel=iint=1` (offset 0). + +`emec`/`emec0` (6) and `stre`/`stiff` (6/21) are per-point vectors in +`umat_main.f:36`, **not** 3-D arrays — so only STATEV needs slicing. + +**Fix:** wire `IEL`, `IINT`, `MI` through the macro; compute +`offset = nstatv_v * ((iint-1) + (*mi)*(iel-1))` and use `statev_old+offset` / +`statev_new+offset` for both the seed copy and `a.statev`. Add a stateful test +at a non-`(1,1)` point so a regression cannot hide at offset 0. + +--- + +## High + +### H1 — Tangent is transposed for major-asymmetric materials +`calculix_interface.h:131` + `tensor_conversion.h:146`, `material_point_evaluator.h:82,181` + +`umat_dispatch` writes `a.ddsdde` (my `ddsdde36`) **column-major**: +`material_point_evaluator.h:82` documents it, `narrow_matrix` +(`tensor_conversion.h:146`) does `host[a + b*n] = C(a,b)`. So +`ddsdde36[i + j*6] == C(i,j)`, and my packing read `ddsdde36[i*6 + j]` is +`C(j,i)` — the transpose. Harmless for a major-symmetric tangent (linear +elastic), silently wrong for non-associative plasticity / damage. The comment +"row- vs column-major indexing … is moot" asserts an unenforced assumption. + +The cross-check test cannot catch it: `ref_ddsdde` is filled by the same +column-major path and read with the same `[i*6+j]`, so both sides carry the same +transpose; the shear test probes only the diagonal `(3,3)`, index-invariant. + +**Fix:** symmetrize exactly as the reference `umat_abaqus.f:335-355` does +(`stiff = (ddsdde(i,j)+ddsdde(j,i))/2`): +`stiff[i+j*(j+1)/2] = 0.5*(ddsdde36[i + j*6] + ddsdde36[j + i*6])`. +For symmetric tangents this is identity; for asymmetric it matches ccx (which +keeps only the symmetric part in `stiff(21)`). Add a test with a known +asymmetric `C` asserting an off-diagonal lands at a specific `stiff` index. + +### H2 — TIME(1)/TIME(2) mapping is wrong (and the comment claims it is right) +`calculix_interface.h:97-99` + +Verified against `umat_abaqus.f:187-188`: +`abqtime(1) = time - dtime`, `abqtime(2) = ttime + time - dtime`, where ccx +`time` = step time at the **end** of the increment, `ttime` = total time at the +**start of the step**. The adapter passes `{time, ttime}`, so TIME(1) is off by +`dtime` and TIME(2) is short by `(time - dtime)`. Correct only on the first +increment of the first step. Any time/rate-dependent model (creep, +viscoelasticity, the `external_scalar_source` time consumers) gets the wrong +absolute time on every later increment. The comment on `:98` is a false claim. + +**Fix:** `dt = dtime?*dtime:0; t = time?*time:0; tt = ttime?*ttime:0;` +`time2 = { t - dt, tt + t - dt }`. Correct the comment. + +### H3 — `iorien != 0` silently ignored; the native hook makes rotation the user's job +`calculix_interface.h:44-49` (macro drops `IORIEN`/`PGAUSS`/`ORAB`, `:156`) + +For `umat_abaqus.f` ccx rotates strain in / stress+stiffness out around the +call, so a UMAT never sees orientation. The native `umat_user` has **no** such +wrapper: `umat_user.f:86-104` requires results in the material frame and tells +the user to call `transformatrix(orab(1,iorien),…)`; `umat_main.f` does no +rotation around the external call. So with `*ORIENTATION` in the deck, +`iorien != 0` reaches the adapter and results come back in the wrong frame — +silently wrong. The header reuses the Abaqus-entry justification, which does not +transfer. + +**Fix:** wire `IORIEN`; `throw fatal_error` when `*iorien != 0` until real +rotation exists. Converts a silent wrong answer into a hard stop; cheap. + +### H4 — Stress/strain measure is small-strain-only, unguarded +`calculix_interface.h:31,44-49` + +`emec` is Green-Lagrange, `stre` is PK2; the models consume/emit them as small +strain. For a `C:E` elastic law this coincidentally *is* St-Venant–Kirchhoff, so +it is correct. Under `NLGEOM` ccx passes large Green-Lagrange strains; the first +inelastic model then feeds finite strain into a small/log-strain return map and +labels the output PK2 — a wrong converged stress, not merely a slow tangent. +The single test (`emec0=0`, one increment) cannot see this. + +**Fix:** document the scope as *geometrically-linear / small-strain only* far +more forcefully than "targets first"; consider a magnitude guard. + +--- + +## Medium + +- **M1 — No real STATEV-flow test.** All tests are single-call linear-elastic, + so the `xstateini→xstate` seeding is a no-op and untested. Add a + multi-increment J2/hardening driver committing `statev_new→statev_old` between + steps, cross-checked vs the direct evaluator. (Also guards C1.) + `test_calculix_interface.cpp` +- **M2 — Nonzero `emec0` never tested.** Every case has `emec0=0`, so + `stran = 2·emec0 = 0` and the stran/dstran split is unverified; only the sum + is exercised. Add a case with `emec0` nonzero in all six slots. +- **M3 — 16-arg positional raw-pointer API.** `calculix_dispatch` abandons the + misuse-resistant `dispatch_args` aggregate the Abaqus side uses; seven + `const double*` in a row make a `time`/`ttime` or `emec`/`emec0` swap compile + cleanly — exactly the class of bug H2 is. Consider a `calculix_args` + aggregate or do the translation inside the macro. `calculix_interface.h:57-64` +- **M4 — `beta` and `ielas` dropped without a guard.** `beta` = + `*INITIAL CONDITIONS,TYPE=STRESS` (`umat_user.f:52`) — a preloaded model is + wrong from step 1; `ielas==1` requests an elastic response ccx will later + mis-use. Wire and guard (fatal if `beta` nonzero) rather than leave un-named. +- **M5 — No error-path coverage.** The sibling suite's `FatalProbe` + + `set_fatal_handler` (`test_umat_interface.cpp:322-335`) is available. Cover: + unknown model → fatal + zeroed outputs; too-few constants → `require_props` + fatal (the `nconst=-kode-100` decode is otherwise happy-path only); a + convergence failure → `pnewdt` cutback propagated back through the adapter. +- **M6 — `emec`/`emec0` inconsistently null-guarded** vs every other pointer; + a null there segfaults rather than degrading. Guard for consistency. + `calculix_interface.h:78-84` + +## Low + +- **L1 — `nstatv_v` not clamped ≥ 0** before flowing to + `umat_interface.h:440`'s `static_cast`, where a negative wraps to a + huge span. Shared with the pre-existing `umat_` path; unreachable via a valid + deck. Clamp `*nstatv > 0 ? *nstatv : 0`. +- **L2 — `-(*kode)-100` is UB if `*kode==INT_MIN`.** Theoretical; widen to + `long long` if desired. `:71` +- **L3 — Fortran hidden-length hardcoded `int`** vs the configurable + `NUMSIM_MATERIALS_FORTRAN_STRLEN` the Abaqus macro exposes. Benign (last arg, + `amat` unused) but inconsistent with the stated lesson. +- **L4 — Magic `6`/`36`/`21`.** A `constexpr` would self-document the packing. +- **L5 — Doc overclaim** `:35` "verified bit-identical against ccx built-in + *ELASTIC" — that verification was the compiled-in `umat_user_` target + (S11=0.027), not this external `.so` path. Soften or cite. +- **L6 — Abaqus concepts behind a CalculiX-named file.** The `umat_dispatch` + reuse is the right seam, but the header should name the inherited + error/cutback semantics (`pnewdt=0.25`, zero-outputs) so a reader doesn't + assume this file owns them. + +--- + +## Suggested sequencing (one PR) + +1. **C1 + H1 + H2 + H3** — the code-correctness fixes; all small, all currently + masked. C1 changes the macro signature (wire `IEL`/`IINT`/`MI`), so do it + first. +2. **M1 + M2 + M5** — the tests that would have caught C1/H1/H2 and prevent + regressions (stateful multi-increment at a non-(1,1) point; nonzero `emec0`; + asymmetric-tangent packing; error paths). +3. **H4 + M3 + M4 + L*** — hardening and documentation; M3 (aggregate) is the + structural change that closes the transposition class M-wide. + +--- + +## Resolution + +All findings applied in `calculix_interface.h` / `test_calculix_interface.cpp`. +Full suite: **293 tests, 0 failures** (282 before, 11 new). + +| # | Fix | +|---|---| +| C1 | `IEL`/`IINT`/`MI` wired through the macro; state sliced at `nstatv·((iint-1) + mi1·(iel-1))`, with the index triple validated (`iint ≤ mi1`, all ≥ 1) and a fatal on inconsistency. | +| H1 | `stiff(21)` now carries the **symmetrized** tangent, `0.5·(D[i+j·6] + D[j+i·6])`, exactly as `umat_abaqus.f:335-355`. Reads the buffer column-major, matching `narrow_matrix`. | +| H2 | TIME rebased onto the start of the increment: `{time-dtime, ttime+time-dtime}`, mirroring `umat_abaqus.f:187-188`. Comment corrected. | +| H3 | `iorien != 0` now zeroes outputs and reports a fatal instead of returning wrong-frame results. | +| H4 | Header carries an explicit **SCOPE — GEOMETRICALLY LINEAR (SMALL STRAIN) ONLY** section stating the PK2/Green-Lagrange pairing is exact only for `C:E`, and that an inelastic model under NLGEOM would be silently wrong. | +| M1 | `IndexesStateByElementAndIntegrationPoint`: 40-increment J2 driven at `iel=2, iint=3`, cross-checked against the Abaqus path, asserting state landed in the right block, that it is non-trivial, and that **no other block was touched**. | +| M2 | The cross-check now uses a nonzero `emec0` in all six slots, so the `stran`/`dstran` split is observable. | +| M3 | Introduced the `calculix_args` aggregate; the 16-arg positional list is gone, closing the transposition class that produced H2. | +| M4 | `beta` wired and guarded (nonzero ⇒ fatal), with a test that zero `beta` still passes. `ielas` wired and documented as ignored. | +| M5 | Error paths covered: unknown model, `kode`-decoded constant shortfall, orientation, initial stress — all via `FatalProbe`. | +| M6 | `emec`/`emec0` null-guarded consistently with the rest; pointer contract documented. | +| L1 | `nstatv` clamped to ≥ 0 before it reaches the `size_t` cast. | +| L2 | `-kode-100` widened to `long long` before negation; count clamped ≥ 0. | +| L3 | Resolved *opposite* to the suggestion: the trailing length is `int` because `call_external_umat_user.c`'s C typedef fixes it — unlike Abaqus, where the Fortran compiler chooses. Documented; deliberately **not** made configurable. | +| L4 | `calculix_ntens` / `calculix_nstiff` replace the magic 6 / 21 / 36. | +| L5 | The "verified bit-identical against ccx built-in *ELASTIC" claim removed — that run was the compiled-in `umat_user_` target, not this `.so` path. | +| L6 | Header documents the inherited `umat_dispatch` error semantics (fatal-zeroes-and-terminates, `PNEWDT = 0.25` cutback, zeroed 6×6 packed through). | + +### Regression tests verified to fail against the old code + +Each was re-run with the fix reverted, to prove it is not vacuous: + +- **C1** → `state written outside this point's block at 4` +- **H1** → `stiff(0,1)`, `stiff(0,2)`, `stiff(1,2)` mismatched +- **H2** → `passing ttime straight through would give 5.0 here` + +### Still open + +- No test forces a genuine convergence failure to assert the `PNEWDT` cutback + round-trips through the adapter. The path is inherited unchanged from + `umat_dispatch` (which the sibling suite covers), so this is a gap in + *adapter-level* coverage only. +- `ielas == 1` (elastic-iteration request) is accepted and ignored. Harmless for + an elastic model; an inelastic one will need to honour it. +- End-to-end validation against a real ccx run through the `@LIB,FUNC` deck path + has not been done for this adapter. diff --git a/include/numsim-materials/umat/calculix_interface.h b/include/numsim-materials/umat/calculix_interface.h new file mode 100644 index 0000000..743622e --- /dev/null +++ b/include/numsim-materials/umat/calculix_interface.h @@ -0,0 +1,321 @@ +#ifndef CALCULIX_INTERFACE_H +#define CALCULIX_INTERFACE_H + +#include +#include +#include +#include + +#include "numsim-materials/umat/umat_interface.h" + +/// The CalculiX external-behaviour entry point (`call_external_umat_user`). +/// +/// CalculiX has two external hooks, and they do NOT share a convention: +/// +/// * `call_external_umat` is reached from umat_abaqus.f AFTER that wrapper has +/// already converted to the Abaqus UMAT convention (STRAN = strain at the +/// start with engineering shear, DSTRAN the increment, DDSDDE a full 6x6, +/// STATEV already sliced to the integration point). A material built with +/// NUMSIM_MATERIALS_DEFINE_UMAT already serves it — no adapter is needed. +/// +/// * `call_external_umat_user` is the NATIVE hook (umat_user.f). It speaks +/// CalculiX's own convention, which this adapter translates to and from: +/// +/// - kode, not nprops: the 4th argument is `kode = -100 - #constants`, so +/// the constant count is `-kode - 100` (umat_user.f:37-43). +/// - emec / emec0 are the Lagrange mechanical strain at the END and START +/// of the increment, component order {11,22,33,12,13,23}, TENSORIAL +/// (no engineering doubling — umat_abaqus.f:280-283 doubles the shear +/// itself before handing strain to an Abaqus UMAT, which is exactly the +/// conversion this adapter performs). +/// - stre is the second Piola-Kirchhoff stress, in on entry / out on exit, +/// plain components with no scaling. +/// - stiff(21) is the upper triangle of the symmetric 6x6 tangent, packed +/// column-major, carrying no engineering factors. +/// - xstateini / xstate are the FULL three-dimensional state arrays +/// `(nstate_, mi(1), #elements)`, NOT a slice for this point. The callee +/// must index them itself — which is why iel, iint and mi are in the +/// signature. See the STATE VARIABLES note below. +/// - TIME must be rebased onto the START of the increment; ccx passes the +/// step time at its END. See the TIME note below. +/// +/// One emitted symbol maps to one registered model (the CalculiX deck selects a +/// specific `@LIB,FUNC` per material), so the model name is baked into the macro +/// rather than taken from `amat`. That also means one FUNC implies one constant +/// set: the registry's props-invariance check is what enforces it, and its +/// message is phrased for Abaqus ("use distinct *MATERIAL names"), which in ccx +/// means a distinct FUNC. +/// +/// SCOPE — GEOMETRICALLY LINEAR (SMALL STRAIN) ONLY. `emec` is the Green-Lagrange +/// strain and `stre` is PK2, while the models below consume and produce these as +/// small-strain quantities. For a linear `C:E` law the two coincide (that pairing +/// IS St-Venant-Kirchhoff), so the elastic case is exactly right. It is NOT +/// merely a slower tangent for anything else: under NLGEOM ccx passes large +/// Green-Lagrange strains, and feeding those to a small- or logarithmic-strain +/// return map produces a wrong CONVERGED stress, silently. Do not use this hook +/// for an inelastic model under finite deformation without first choosing a +/// conjugate stress-strain pair and converting here. +/// +/// INHERITED ERROR SEMANTICS. Fault handling belongs to umat_dispatch, not to +/// this file: a setup fault zeroes the outputs and terminates via the fatal +/// handler, and any other exception zeroes the outputs and requests a cutback +/// with PNEWDT = 0.25 — which is a valid ccx pnewdt (0 < pnewdt < 1), so the +/// Abaqus-side convention transfers unchanged. The zeroed 6x6 is packed into +/// stiff(21) like any other result. +/// +/// NOT handled, and rejected rather than ignored: a local material orientation +/// (`iorien != 0`) and a nonzero initial stress (`beta`). Both would otherwise be +/// silently wrong — see the guards in calculix_dispatch. Deliberately IGNORED: +/// the deformation gradients and Jacobians (small-strain scope, above), +/// temperature (constants are read once when the graph is built), and `ielas`, +/// whose elastic-iteration request an inelastic model would need to honour. +namespace numsim::materials::umat { + +/// The native hook is always full 3D: six components, and the symmetric 6x6 +/// tangent packed as 21 upper-triangular entries. +inline constexpr std::size_t calculix_ntens = 6; +inline constexpr std::size_t calculix_nstiff = + calculix_ntens * (calculix_ntens + 1) / 2; + +/// Every argument the CalculiX shim forwards, named. +/// +/// The Abaqus side funnels through dispatch_args for the same reason: a long +/// positional list of same-typed pointers lets a transposition compile cleanly, +/// and `time`/`ttime` or `emec`/`emec0` are exactly the pairs that would be +/// swapped. Scalars are held by value; CalculiX passes them by reference as +/// Fortran always does, and the macro dereferences them, so a null there would +/// be a CalculiX defect rather than a case to handle. +struct calculix_args { + const char* amat{nullptr}; + int iel{0}; ///< element number, 1-based + int iint{0}; ///< integration point number, 1-based + int kode{0}; ///< -100 - #constants + const double* mprops{nullptr}; + const double* emec{nullptr}; ///< tensorial strain at the END + const double* emec0{nullptr}; ///< tensorial strain at the START + const double* beta{nullptr}; ///< initial (residual) stress + double dtime{0}; + double time{0}; ///< STEP time at the END of the increment + double ttime{0}; ///< TOTAL time at the START of the step + int icmd{0}; ///< 3 => stress only, leave stiff alone + int ielas{0}; + int mi1{0}; ///< mi(1): max integration points per element + int nstatv{0}; ///< nstate_: state variables PER integration point + const double* statev_old{nullptr}; ///< xstateini, FULL array + double* statev_new{nullptr}; ///< xstate, FULL array + double* stress{nullptr}; + double* stiff{nullptr}; ///< 21 packed entries + int iorien{0}; + double* pnewdt{nullptr}; +}; + +/// Zero what CalculiX reads back, in CalculiX's own shapes. +inline void calculix_zero_outputs(double* stress, double* stiff) noexcept { + if (stress) + for (std::size_t i = 0; i < calculix_ntens; ++i) stress[i] = 0.0; + if (stiff) + for (std::size_t i = 0; i < calculix_nstiff; ++i) stiff[i] = 0.0; +} + +/// Shared implementation behind the CalculiX external symbol. +/// +/// Everything here is fixed-size stack work with no allocation, so it cannot +/// throw before reaching umat_dispatch, which owns the try/catch that keeps a +/// C++ exception from unwinding into CalculiX's Fortran. The guards below +/// report through the same fatal handler rather than throwing, for that reason. +template +void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept { + using T = typename Traits::value_type; + static_assert(std::is_same_v, + "the CalculiX external interface is double-precision (ccxreal)"); + + // A local orientation is NOT applied around this hook. umat_abaqus.f rotates + // strain in and stress/stiffness out around its call (lines 202-276, 305-330), + // but umat_main.f does no such thing around the native or external umat_user + // call: umat_user.f:86-104 requires results in the MATERIAL frame and tells + // the user to rotate with transformatrix(). Ignoring iorien would therefore + // return results in the wrong frame with no symptom, so refuse instead. + if (a.iorien != 0) { + calculix_zero_outputs(a.stress, a.stiff); + report_fatal(model_name, + "numsim CalculiX: *ORIENTATION (iorien != 0) is not supported " + "by this external behaviour — the native umat_user hook makes " + "frame rotation the material's own responsibility, and it is " + "not implemented here"); + return; + } + + // *INITIAL CONDITIONS,TYPE=STRESS (umat_user.f:52). Dropping it would make a + // preloaded model wrong from the first step, again with no symptom. + if (a.beta) { + for (std::size_t i = 0; i < calculix_ntens; ++i) { + if (a.beta[i] != 0.0) { + calculix_zero_outputs(a.stress, a.stiff); + report_fatal(model_name, + "numsim CalculiX: a nonzero initial stress (*INITIAL " + "CONDITIONS,TYPE=STRESS) is not supported by this external " + "behaviour"); + return; + } + } + } + + if (!a.emec || !a.emec0) { + calculix_zero_outputs(a.stress, a.stiff); + report_fatal(model_name, + "numsim CalculiX: the mechanical strain arrays are null"); + return; + } + + // Widened before negating so kode == INT_MIN cannot trap, and clamped because + // a malformed kode must not become a negative count. + const long long nconst_ll = -static_cast(a.kode) - 100; + const int nconst = nconst_ll > 0 ? static_cast(nconst_ll) : 0; + const int nstatv = a.nstatv > 0 ? a.nstatv : 0; + + // STATE VARIABLES. xstateini/xstate are the FULL arrays + // `real*8 xstate(nstate_, mi(1), *)` (umat_main.f:40), and umat_main.f passes + // their BASE to this hook (line 233) — unlike umat_abaqus.f:295, which slices + // `xstate(1,iint,iel)` before calling. So the offset of this point's block is + // the Fortran column-major stride. Getting this wrong does not fail loudly: it + // makes every integration point share element 1 / point 1's state. + std::size_t point_offset = 0; + if (nstatv > 0) { + if (a.iel < 1 || a.iint < 1 || a.mi1 < 1 || a.iint > a.mi1) { + calculix_zero_outputs(a.stress, a.stiff); + report_fatal(model_name, + "numsim CalculiX: iel/iint/mi(1) are inconsistent, so the " + "state-variable block for this integration point cannot be " + "located"); + return; + } + point_offset = static_cast(nstatv) * + (static_cast(a.iint - 1) + + static_cast(a.mi1) * + static_cast(a.iel - 1)); + } + + // Native emec is TENSORIAL; umat_dispatch consumes ENGINEERING shear + // (strain_from_buffer halves slots 3-5). Doubling the shear here makes the + // round-trip land back on the tensorial strain the models expect — the same + // conversion umat_abaqus.f:280-283 applies before an Abaqus UMAT. STRAN is the + // strain at the START (emec0); DSTRAN the increment, so STRAN+DSTRAN == emec. + double stran[calculix_ntens]; + double dstran[calculix_ntens]; + for (std::size_t i = 0; i < calculix_ntens; ++i) { + const double shear = (i < 3) ? 1.0 : 2.0; + stran[i] = shear * a.emec0[i]; + dstran[i] = shear * (a.emec[i] - a.emec0[i]); + } + + // The evaluator reads and writes ONE state array in place; CalculiX splits it + // into a read-only xstateini and a write-only xstate. Seed this point's block + // in the output from the committed input so the read sees last increment's + // state and the write lands where CalculiX expects the update. + double* statev_point = nullptr; + if (nstatv > 0 && a.statev_new) { + statev_point = a.statev_new + point_offset; + if (a.statev_old) + std::copy_n(a.statev_old + point_offset, nstatv, statev_point); + } + + // umat_dispatch fills a full 6x6; CalculiX wants only the packed upper + // triangle. Note this buffer is COLUMN-major (material_point_evaluator.h:82, + // narrow_matrix writes host[a + b*n]), which the packing below relies on. + double ddsdde36[calculix_ntens * calculix_ntens] = {0.0}; + + // TIME must be rebased onto the START of the increment. ccx hands over the + // step time at the END of the increment plus the total time at the start of + // the STEP, so mirror umat_abaqus.f:187-188 exactly: + // abqtime(1) = time - dtime + // abqtime(2) = ttime + time - dtime + // Passing {time, ttime} instead is correct only on the first increment of the + // first step, and silently wrong for every rate- or time-dependent model after + // that. + const double time2[2] = {a.time - a.dtime, + a.ttime + a.time - a.dtime}; + + dispatch_args d; + d.stress = a.stress; + d.statev = statev_point; + d.ddsdde = ddsdde36; + d.stran = stran; + d.dstran = dstran; + d.time = time2; + d.dtime = a.dtime; + d.pnewdt = a.pnewdt; + d.props = a.mprops; + d.nprops = nconst; + d.cmname = model_name; + d.cmname_len = std::char_traits::length(model_name); + d.ndi = 3; + d.nshr = 3; + d.ntens = static_cast(calculix_ntens); + d.nstatv = nstatv; + + umat_dispatch(d); + + // Pack into stiff(21): column-major upper triangle, so the 0-based (i, j) with + // i <= j lands at i + j*(j+1)/2. + // + // SYMMETRIZED, exactly as the reference umat_abaqus.f:335-355 does + // (stiff(2) = (ddsdde(1,2)+ddsdde(2,1))/2, and so on). Two reasons: ddsdde36 + // is column-major, so reading it row-major would transpose every off-diagonal + // for a major-ASYMMETRIC tangent (non-associative flow, damage); and stiff(21) + // has no room for an antisymmetric part anyway, so averaging is what CalculiX + // itself keeps. For a symmetric tangent this is the identity. + // + // icmd == 3 requests stress only and leaves stiff untouched, as umat_abaqus.f + // does; on an error path umat_dispatch has zeroed ddsdde36 already. + if (a.stiff && a.icmd != 3) { + for (std::size_t j = 0; j < calculix_ntens; ++j) + for (std::size_t i = 0; i <= j; ++i) + a.stiff[i + j * (j + 1) / 2] = + 0.5 * (ddsdde36[i + j * calculix_ntens] + + ddsdde36[j + i * calculix_ntens]); + } +} + +} // namespace numsim::materials::umat + +/// Emit a CalculiX-callable external-behaviour symbol @p FUNC bound to the +/// registered model @p MODELNAME. Place this in exactly ONE translation unit of +/// the shared library CalculiX dlopen's for `*MATERIAL, NAME=@LIB,FUNC`. +/// +/// The signature is the `calculixptr` prototype from CalculiX's +/// call_external_umat_user.c: every argument is passed by pointer (ccxint == +/// int, ccxreal == double), and the final `int` is the hidden Fortran length of +/// the leading character string. That trailing width is `int` because the C +/// typedef in call_external_umat_user.c fixes it as `int` — unlike the Abaqus +/// entry, where the Fortran compiler chooses and +/// NUMSIM_MATERIALS_FORTRAN_STRLEN exists to track it. Do not make it +/// configurable here; the C prototype is the contract. +/// +/// Arguments the small-strain path does not use are named in comments so their +/// absence reads as deliberate, not forgotten. +#define NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR(TRAITS, FUNC, MODELNAME) \ + extern "C" void FUNC( \ + const char* AMAT, const int* IEL, const int* IINT, \ + const int* KODE, const double* MPROPS, const double* EMEC, \ + const double* EMEC0, const double* BETA, const double* /*XOKL*/, \ + const double* /*VOJ*/, const double* /*XKL*/, const double* /*VJ*/, \ + const int* /*ITHERMAL*/, const double* /*T1L*/, const double* DTIME, \ + const double* TIME, const double* TTIME, const int* ICMD, \ + const int* IELAS, const int* MI, const int* NSTATV, \ + const double* STATEV0, double* STATEV1, double* STRESS, double* STIFF, \ + const int* IORIEN, const double* /*PGAUSS*/, const double* /*ORAB*/, \ + double* PNEWDT, const int* /*IPKON*/, const int /*SIZE*/) { \ + ::numsim::materials::umat::calculix_args a; \ + a.amat = AMAT; a.iel = *IEL; a.iint = *IINT; \ + a.kode = *KODE; a.mprops = MPROPS; \ + a.emec = EMEC; a.emec0 = EMEC0; a.beta = BETA; \ + a.dtime = *DTIME; a.time = *TIME; a.ttime = *TTIME; \ + a.icmd = *ICMD; a.ielas = *IELAS; a.mi1 = *MI; \ + a.nstatv = *NSTATV; \ + a.statev_old = STATEV0; a.statev_new = STATEV1; \ + a.stress = STRESS; a.stiff = STIFF; \ + a.iorien = *IORIEN; a.pnewdt = PNEWDT; \ + ::numsim::materials::umat::calculix_dispatch(a, MODELNAME); \ + } + +#endif // CALCULIX_INTERFACE_H diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c163bc0..0caebb8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -24,6 +24,7 @@ add_numsim_test(test_statev_map test_statev_map.cpp) add_numsim_test(test_material_point_evaluator test_material_point_evaluator.cpp) add_numsim_test(test_plane_stress_evaluator test_plane_stress_evaluator.cpp) add_numsim_test(test_umat_interface test_umat_interface.cpp) +add_numsim_test(test_calculix_interface test_calculix_interface.cpp) add_numsim_test(test_set_parameter test_set_parameter.cpp) add_numsim_test(test_tangent_generator test_tangent_generator.cpp) add_numsim_test(test_json_model test_json_model.cpp) diff --git a/tests/test_calculix_interface.cpp b/tests/test_calculix_interface.cpp new file mode 100644 index 0000000..da3daf8 --- /dev/null +++ b/tests/test_calculix_interface.cpp @@ -0,0 +1,735 @@ +#include + +#include +#include +#include +#include +#include + +#include + +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/materials/j2_plasticity.h" +#include "numsim-materials/materials/linear_elasticity.h" +#include "numsim-materials/materials/linear_isotropic_hardening.h" +#include "numsim-materials/solvers/local_newton.h" +#include "numsim-materials/umat/calculix_interface.h" +#include "numsim-materials/umat/external_state_source.h" +#include "numsim-materials/umat/umat_interface.h" + +// Both entry points live in this one translation unit: umat_ is the reference +// path the CalculiX adapter is checked against, and the clx_* symbols are the +// adapter under test. They are distinct extern "C" symbols, so they coexist. +NUMSIM_MATERIALS_DEFINE_UMAT(numsim::materials::material_policy_default) +NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( + numsim::materials::material_policy_default, clx_linear_elastic_, "LINELAS") +NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( + numsim::materials::material_policy_default, clx_j2_, "J2CLX") +NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( + numsim::materials::material_policy_default, clx_asym_, "ASYMTANGENT") +NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( + numsim::materials::material_policy_default, clx_missing_, "NOSUCHMODEL") +NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( + numsim::materials::material_policy_default, clx_time_, "TIMEPROBE") + +namespace { + +namespace nm = numsim::materials; +namespace u = numsim::materials::umat; + +using policy = nm::material_policy_default; +using T = policy::value_type; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; +using registry = u::umat_registry; + +// Deliberately asymmetric moduli so K and G cannot alias each other. +constexpr T K = 166.67; +constexpr T G = 76.92; +constexpr T sigma_0 = 50.0; +constexpr T H_mod = 1000.0; + +/// props[0] = K, props[1] = G. One builder, constants entirely from the deck. +void build_deck_elastic(ctx_type& ctx, std::span props) { + u::require_props(props, 2, "linelas"); + param_type p; + p.insert("name", "stepper"); + ctx.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("strain_producer_name", "stepper"); + p.insert("K", props[0]); + p.insert("G", props[1]); + ctx.create>(p); + ctx.finalize(); +} + +/// A genuinely stateful model, so the xstateini -> xstate plumbing has real +/// state to carry. Mirrors the sibling suite's J2 build. +void build_j2(ctx_type& ctx, std::span /*props*/) { + param_type p; + p.insert("name", "stepper"); + ctx.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("strain_producer_name", "stepper"); + p.insert("K", K); + p.insert("G", G); + ctx.create>(p); + p.clear(); + p.insert("name", "solver"); + ctx.create>(p); + p.clear(); + p.insert("name", "hardening"); + p.insert("source", "j2"); + p.insert("K", H_mod); + ctx.create>(p); + p.clear(); + p.insert("name", "j2"); + p.insert("K", K); + p.insert("hardening_source", "hardening"); + p.insert("strain_source", "stepper"); + p.insert("solver_source", "solver"); + p.insert("G", G); + p.insert("sigma_0", sigma_0); + ctx.create>(p); + ctx.finalize(); +} + +// The six canonical slots {11,22,33,12,13,23} as (i,j) index pairs. +constexpr int slot_i[6] = {0, 1, 2, 0, 0, 1}; +constexpr int slot_j[6] = {0, 1, 2, 1, 2, 2}; + +/// The value this probe puts at 6x6 entry (a, b). Deliberately NOT symmetric: +/// asym_value(a,b) != asym_value(b,a) for a != b, which is what makes a +/// transposed or one-sided read of the tangent observable. +constexpr T asym_value(int a, int b) { + return static_cast(100 * (a + 1) + (b + 1)); +} + +/// Publishes a known, minor-symmetric but major-ASYMMETRIC tangent. +/// +/// Minor symmetry is required (tangent_to_buffer asserts it); major symmetry is +/// not, and is exactly what a non-associative model lacks. Stress is left zero: +/// this probe exists only to pin the tangent packing. +template +class asym_tangent_material final + : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + using tensor2 = tmech::tensor; + using tensor4 = tmech::tensor; + + template + explicit asym_tangent_material(Args&&... args) + : base(std::forward(args)...), + m_stress(base::template add_output( + "stress", &asym_tangent_material::compute)), + m_tangent(base::template add_output("tangent")) { + fill_tangent(); + } + + static input_parameter_controller parameters() { return base::parameters(); } + + void compute() { + m_stress.fill(0.0); + fill_tangent(); + } + +private: + void fill_tangent() { + m_tangent.fill(0.0); + for (int a = 0; a < 6; ++a) { + for (int b = 0; b < 6; ++b) { + const value_type v = asym_value(a, b); + const int i = slot_i[a], j = slot_j[a]; + const int k = slot_i[b], l = slot_j[b]; + // Every minor-symmetric permutation gets the same value. + m_tangent(i, j, k, l) = v; + m_tangent(j, i, k, l) = v; + m_tangent(i, j, l, k) = v; + m_tangent(j, i, l, k) = v; + } + } + } + + tensor2& m_stress; + tensor4& m_tangent; +}; + +/// Reports the time it was bound with, so the TIME rebasing is observable from +/// outside. Without a probe like this, nothing reaching a material through the +/// shim depends on time at all and a wrong TIME(2) is undetectable at the ABI. +template +class time_probe_material final + : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + using tensor2 = tmech::tensor; + using tensor4 = tmech::tensor; + + template + explicit time_probe_material(Args&&... args) + : base(std::forward(args)...), + m_stress(base::template add_output( + "stress", &time_probe_material::compute)), + m_tangent(base::template add_output("tangent")), + m_time(base::template add_input_history( + nm::connection_source{"clock", "state"}, nm::EdgeKind::Global)) {} + + static input_parameter_controller parameters() { return base::parameters(); } + + void compute() { + m_stress.fill(0.0); + m_stress(0, 0) = m_time.old_value(); // start of incr + m_stress(1, 1) = m_time.new_value(); // end of incr + m_stress(2, 2) = m_time.new_value() - m_time.old_value(); // == dtime + } + +private: + tensor2& m_stress; + tensor4& m_tangent; + const nm::input_history& m_time; +}; + +void build_time_probe(ctx_type& ctx, std::span /*props*/) { + param_type p; + p.insert("name", "stepper"); + ctx.create>(p); + p.clear(); + p.insert("name", "clock"); + ctx.create>(p); + p.clear(); + p.insert("name", "probe"); + ctx.create>(p); + ctx.finalize(); +} + +void build_asym(ctx_type& ctx, std::span /*props*/) { + param_type p; + p.insert("name", "stepper"); + ctx.create>(p); + p.clear(); + p.insert("name", "probe"); + ctx.create>(p); + ctx.finalize(); +} + +struct Registration { + Registration() { + registry::config el; + el.strain_source = "stepper"; + el.stress_source = "elastic"; + registry::instance().register_model("LINELAS", build_deck_elastic, el); + + registry::config j2; + j2.strain_source = "stepper"; + j2.stress_source = "j2"; + registry::instance().register_model("J2CLX", build_j2, j2); + // Same graph, reached through the Abaqus entry for the cross-check. + registry::instance().register_model("J2REF", build_j2, j2); + + registry::config as; + as.strain_source = "stepper"; + as.stress_source = "probe"; + registry::instance().register_model("ASYMTANGENT", build_asym, as); + + registry::config tp; + tp.strain_source = "stepper"; + tp.stress_source = "probe"; + tp.time_source = "clock"; + registry::instance().register_model("TIMEPROBE", build_time_probe, tp); + } +}; +const Registration registration; + +/// CMNAME as Fortran passes it: character*80, blank padded, no NUL. +struct fortran_name { + char buf[80]; + explicit fortran_name(const std::string& s) { + for (std::size_t i = 0; i < 80; ++i) buf[i] = ' '; + for (std::size_t i = 0; i < s.size() && i < 80; ++i) buf[i] = s[i]; + } +}; + +/// A fatal fault must terminate the analysis. The handler is replaced so the +/// test can observe it instead of the runner being killed by XIT/abort. +struct FatalProbe { + static inline std::string last; + static inline int count = 0; + static void handler(const char* msg) { + last = msg; + ++count; + } + FatalProbe() { + last.clear(); + count = 0; + u::set_fatal_handler(&handler); + } + ~FatalProbe() { u::set_fatal_handler(nullptr); } +}; + +using clx_fn = void (*)(const char*, const int*, const int*, const int*, + const double*, const double*, const double*, + const double*, const double*, const double*, + const double*, const double*, const int*, const double*, + const double*, const double*, const double*, const int*, + const int*, const int*, const int*, const double*, + double*, double*, double*, const int*, const double*, + const double*, double*, const int*, const int); + +/// Everything a CalculiX call needs, with the defaults a well-formed solid-3D +/// call uses. Only what a test actually varies is set. +struct clx_call { + int iel = 1; + int iint = 1; + int mi1 = 1; + int nstatv = 0; + int icmd = 0; + int ielas = 0; + int iorien = 0; + T dtime = 1.0; + T time = 0.0; // step time at END of increment + T ttime = 0.0; // total time at START of step + const T* beta = nullptr; + const T* mprops = nullptr; + int nconst = 0; + T pnewdt = -1.0; + + void run(clx_fn fn, const std::string& amat, const T* emec, const T* emec0, + const T* statev_old, T* statev_new, T* stress, T* stiff) { + const fortran_name am(amat); + const int kode = -100 - nconst; + const int ipkon = 0; + const T zero6[6] = {0, 0, 0, 0, 0, 0}; + const T ident[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + const T voj = 1, vj = 1, t1l = 0; + const T pgauss[3] = {0, 0, 0}; + const T orab[7] = {0, 0, 0, 0, 0, 0, 0}; + const int ithermal = 0; + const T* beta_p = beta ? beta : zero6; + const T no_props = 0; + const T* props_p = mprops ? mprops : &no_props; + + fn(am.buf, &iel, &iint, &kode, props_p, emec, emec0, beta_p, ident, &voj, + ident, &vj, &ithermal, &t1l, &dtime, &time, &ttime, &icmd, &ielas, &mi1, + &nstatv, statev_old, statev_new, stress, stiff, &iorien, pgauss, orab, + &pnewdt, &ipkon, 80); + } +}; + +/// Drive the Abaqus umat_ reference for a solid-3D point. STRAN/DSTRAN are the +/// ENGINEERING-shear strain at the start and its increment. +void call_umat_reference(const std::string& name, T* stress, T* statev, + T* ddsdde, const T* stran, const T* dstran, + int nstatv, const T* props, int nprops, + T total_time = 0.0, T dtime_in = 1.0) { + const fortran_name cm(name); + T sse = 0, spd = 0, scd = 0, rpl = 0, drpldt = 0, pnewdt = -1; + T ddsddt[6] = {0}, drplde[6] = {0}; + const T time[2] = {total_time, total_time}; + const T dtime = dtime_in; + const T temp = 0, dtemp = 0, predef = 0, dpred = 0; + const T coords[3] = {0, 0, 0}; + const T identity[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + const T celent = 1.0; + int ndi = 3, nshr = 3, ntens = 6, nstatv_v = nstatv, nprops_v = nprops; + const int noel = 1, npt = 1, layer = 1, kspt = 1, jstep = 1, kinc = 1; + const T no_props = 0; + const T* props_p = props ? props : &no_props; + + umat_(stress, statev, ddsdde, &sse, &spd, &scd, &rpl, ddsddt, drplde, &drpldt, + stran, dstran, time, &dtime, &temp, &dtemp, &predef, &dpred, cm.buf, + &ndi, &nshr, &ntens, &nstatv_v, props_p, &nprops_v, coords, identity, + &pnewdt, &celent, identity, identity, &noel, &npt, &layer, &kspt, + &jstep, &kinc, 80); +} + +// --------------------------------------------------------------------------- +// Convention translation +// --------------------------------------------------------------------------- + +/// The load-bearing claim: the CalculiX adapter and the Abaqus entry, given the +/// SAME physical strain expressed in their respective conventions, produce the +/// same stress and the same tangent. A nonzero shear component makes the +/// tensorial->engineering doubling observable, and a NONZERO start strain makes +/// the stran/dstran split observable — with emec0 = 0 a swapped or dropped +/// emec0 is invisible. +TEST(CalculiXInterface, MatchesTheAbaqusEntryForTheSamePhysicalStrain) { + const T props[2] = {K, G}; + const int nconst = 2; + const int nstatv = static_cast( + registry::instance().nstatv("LINELAS", std::span(props, 2))); + + // Both start and end strain are nonzero in all six slots. + const T emec0[6] = {2.0e-4, 1.0e-4, -3.0e-4, 1.5e-4, 2.5e-4, -2.0e-4}; + const T emec[6] = {1.0e-3, -4.0e-4, 2.0e-4, 3.0e-4, -1.5e-4, 5.0e-4}; + + // --- Abaqus reference: engineering shear, STRAN = start, DSTRAN = increment. + T stran[6], dstran[6]; + for (int i = 0; i < 6; ++i) { + const T s = (i < 3) ? 1.0 : 2.0; + stran[i] = s * emec0[i]; + dstran[i] = s * (emec[i] - emec0[i]); + } + std::vector ref_statev(nstatv > 0 ? nstatv : 1, 0.0); + T ref_stress[6] = {0, 0, 0, 0, 0, 0}; + T ref_ddsdde[36] = {0}; + call_umat_reference("LINELAS", ref_stress, ref_statev.data(), ref_ddsdde, + stran, dstran, nstatv, props, nconst); + + // --- CalculiX adapter: native tensorial strain, split STATEV, packed stiff. + std::vector clx_old(nstatv > 0 ? nstatv : 1, 0.0); + std::vector clx_new(nstatv > 0 ? nstatv : 1, 0.0); + T clx_stress[6] = {0, 0, 0, 0, 0, 0}; + T clx_stiff[21] = {0}; + clx_call c; + c.nstatv = nstatv; + c.mprops = props; + c.nconst = nconst; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, clx_old.data(), + clx_new.data(), clx_stress, clx_stiff); + + for (int i = 0; i < 6; ++i) + EXPECT_NEAR(clx_stress[i], ref_stress[i], 1e-12) << "stress " << i; + + // ddsdde is COLUMN-major, so the reference entry (i, j) is ref[i + j*6]. + for (int j = 0; j < 6; ++j) + for (int i = 0; i <= j; ++i) + EXPECT_NEAR(clx_stiff[i + j * (j + 1) / 2], ref_ddsdde[i + j * 6], 1e-12) + << "stiff(" << i << "," << j << ")"; +} + +/// The shear conversion is not a no-op: feeding emec directly (tensorial) where +/// the adapter must double it would halve every shear-driven stress. +TEST(CalculiXInterface, TensorialShearIsConvertedToEngineering) { + const T props[2] = {K, G}; + const T e12 = 2.5e-3; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {0, 0, 0, e12, 0, 0}; + + T statev_old[1] = {0}, statev_new[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + clx_call c; + c.mprops = props; + c.nconst = 2; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, statev_old, statev_new, + stress, stiff); + + // sigma_12 = 2 * G * eps_12 (tensorial) = G * gamma_12. Dropping the + // conversion would give exactly half of this. + EXPECT_NEAR(stress[3], 2.0 * G * e12, 1e-12); + EXPECT_NEAR(stress[0], 0.0, 1e-12); + EXPECT_NEAR(stress[1], 0.0, 1e-12); + EXPECT_NEAR(stress[2], 0.0, 1e-12); + // The shear-shear tangent is G, with no engineering factor on the tangent. + EXPECT_NEAR(stiff[3 + 3 * (3 + 1) / 2], G, 1e-9); +} + +/// stiff(21) must carry the SYMMETRIZED tangent, as umat_abaqus.f:335-355 does. +/// +/// The buffer umat_dispatch fills is column-major, so a row-major read would +/// hand CalculiX C(j,i) instead of the average — invisible for every symmetric +/// material, which is why this probe is deliberately major-asymmetric. +TEST(CalculiXInterface, PacksTheSymmetrizedTangentIntoStiff21) { + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; + T statev_old[1] = {0}, statev_new[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + + clx_call c; + c.run(&clx_asym_, "ASYMTANGENT", emec, emec0, statev_old, statev_new, stress, + stiff); + + for (int j = 0; j < 6; ++j) { + for (int i = 0; i <= j; ++i) { + const T expected = 0.5 * (asym_value(i, j) + asym_value(j, i)); + EXPECT_NEAR(stiff[i + j * (j + 1) / 2], expected, 1e-9) + << "stiff(" << i << "," << j << ")"; + if (i != j) { + // Pin the bug this test exists for: a one-sided read would land on + // exactly one of these, and the average equals neither. + EXPECT_NE(expected, asym_value(i, j)); + EXPECT_NE(expected, asym_value(j, i)); + } + } + } +} + +// --------------------------------------------------------------------------- +// State variables: the full-array indexing contract +// --------------------------------------------------------------------------- + +/// CalculiX passes the WHOLE xstateini/xstate arrays — `(nstate_, mi(1), +/// #elements)`, umat_main.f:40 — and their base pointer, not a slice for this +/// point (contrast umat_abaqus.f:295). The adapter must therefore index by +/// (iint, iel) itself. +/// +/// Driven at element 2, integration point 3, this fails outright if the adapter +/// works at offset 0: the trajectory would be read back from, and written to, +/// the wrong block. The untouched blocks are checked too, because writing the +/// right values into the wrong place AND leaving neighbours dirty are different +/// bugs. +TEST(CalculiXInterface, IndexesStateByElementAndIntegrationPoint) { + const int nstatv = static_cast(registry::instance().nstatv("J2CLX")); + ASSERT_GT(nstatv, 0) << "this test is meaningless without real state"; + + constexpr int mi1 = 4; // integration points per element + constexpr int nelem = 3; // elements + constexpr int iel = 2; // 1-based + constexpr int iint = 3; // 1-based + const std::size_t total = static_cast(nstatv) * mi1 * nelem; + const std::size_t offset = + static_cast(nstatv) * ((iint - 1) + mi1 * (iel - 1)); + + // A strain ramp well past yield, so state genuinely accumulates. TENSORIAL, + // so the shear slot is half the engineering increment the sibling suite uses. + const T de[6] = {0.01, -0.0025, 0.0, 0.0025, 0.0, 0.0}; + constexpr int steps = 40; + + // --- CalculiX path, full arrays, driven at (iel=2, iint=3). --- + std::vector st_old(total, 0.0), st_new(total, 0.0); + T clx_stress[6] = {0, 0, 0, 0, 0, 0}; + T clx_stiff[21] = {0}; + T emec0[6] = {0, 0, 0, 0, 0, 0}; + T emec[6] = {0, 0, 0, 0, 0, 0}; + + for (int s = 0; s < steps; ++s) { + for (int i = 0; i < 6; ++i) { + emec0[i] = emec[i]; + emec[i] = emec0[i] + de[i]; + } + clx_call c; + c.iel = iel; + c.iint = iint; + c.mi1 = mi1; + c.nstatv = nstatv; + c.run(&clx_j2_, "J2CLX", emec, emec0, st_old.data(), st_new.data(), + clx_stress, clx_stiff); + // CalculiX commits xstate -> xstateini between increments. + st_old = st_new; + } + + // --- Abaqus reference, one point, its own array. --- + std::vector ref_statev(nstatv, 0.0); + T ref_stress[6] = {0, 0, 0, 0, 0, 0}; + T ref_ddsdde[36] = {0}; + T stran[6] = {0, 0, 0, 0, 0, 0}; + T dstran[6]; + for (int i = 0; i < 6; ++i) dstran[i] = (i < 3 ? 1.0 : 2.0) * de[i]; + for (int s = 0; s < steps; ++s) { + call_umat_reference("J2REF", ref_stress, ref_statev.data(), ref_ddsdde, + stran, dstran, nstatv, nullptr, 0); + for (int i = 0; i < 6; ++i) stran[i] += dstran[i]; + } + + for (int i = 0; i < 6; ++i) + EXPECT_NEAR(clx_stress[i], ref_stress[i], 1e-9) << "stress " << i; + + // The state landed in THIS point's block... + for (int i = 0; i < nstatv; ++i) + EXPECT_NEAR(st_new[offset + i], ref_statev[i], 1e-9) + << "statev " << i << " at (iel=2, iint=3)"; + + // ...and plasticity actually happened, so the block is not trivially zero. + T norm = 0; + for (int i = 0; i < nstatv; ++i) norm += std::abs(st_new[offset + i]); + EXPECT_GT(norm, 1e-8) << "no state accumulated; the test proves nothing"; + + // ...and nowhere else was touched. + for (std::size_t k = 0; k < total; ++k) { + if (k >= offset && k < offset + static_cast(nstatv)) continue; + EXPECT_EQ(st_new[k], 0.0) << "state written outside this point's block at " + << k; + } +} + +// --------------------------------------------------------------------------- +// TIME +// --------------------------------------------------------------------------- + +/// CalculiX hands over the step time at the END of the increment and the total +/// time at the START of the step; the Abaqus convention wants both rebased onto +/// the START of the increment (umat_abaqus.f:187-188). Passing {time, ttime} +/// through is right only on the first increment of the first step. +/// +/// The probe reports the bound times directly, so the exact rebasing is +/// asserted rather than inferred. The material sees the TOTAL time spanning the +/// increment: [ttime + time - dtime, ttime + time]. +TEST(CalculiXInterface, RebasesTimeOntoTheStartOfTheIncrement) { + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-4, 0, 0, 0, 0, 0}; + + auto times_seen = [&](T step_time_end, T total_time_step_start, T dt) { + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + clx_call c; + c.dtime = dt; + c.time = step_time_end; + c.ttime = total_time_step_start; + c.run(&clx_time_, "TIMEPROBE", emec, emec0, so, sn, stress, stiff); + // stress = {t_start, t_end, dt} in slots 0, 1, 2. + return std::array{stress[0], stress[1], stress[2]}; + }; + + // First increment of the first step: time = dtime, ttime = 0. This is the one + // case where passing {time, ttime} straight through would also be correct. + { + const auto t = times_seen(/*time=*/0.1, /*ttime=*/0.0, /*dt=*/0.1); + EXPECT_NEAR(t[0], 0.0, 1e-12) << "total time at the START of the increment"; + EXPECT_NEAR(t[1], 0.1, 1e-12) << "total time at the END of the increment"; + EXPECT_NEAR(t[2], 0.1, 1e-12) << "dtime"; + } + + // A later increment of a later step — where the naive mapping is wrong. + // ccx: step time ends at 0.9, the step began at total time 5.0, dt = 0.1. + // So the increment spans total time [5.0 + 0.9 - 0.1, 5.0 + 0.9] = [5.8, 5.9]. + { + const auto t = times_seen(/*time=*/0.9, /*ttime=*/5.0, /*dt=*/0.1); + EXPECT_NEAR(t[0], 5.8, 1e-12) + << "passing ttime straight through would give 5.0 here"; + EXPECT_NEAR(t[1], 5.9, 1e-12); + EXPECT_NEAR(t[2], 0.1, 1e-12); + } +} + +// --------------------------------------------------------------------------- +// Guards and error paths +// --------------------------------------------------------------------------- + +/// A local orientation is not applied around the native hook, so silently +/// ignoring it would return results in the wrong frame. It must be refused. +TEST(CalculiXInterface, RefusesALocalOrientation) { + FatalProbe probe; + const T props[2] = {K, G}; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {9, 9, 9, 9, 9, 9}; + T stiff[21]; + for (int i = 0; i < 21; ++i) stiff[i] = 9.0; + + clx_call c; + c.mprops = props; + c.nconst = 2; + c.iorien = 1; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, stiff); + + EXPECT_EQ(FatalProbe::count, 1); + EXPECT_NE(FatalProbe::last.find("ORIENTATION"), std::string::npos); + // Outputs are zeroed, so a handler that returns cannot leave the solver + // consuming the caller's uninitialised buffer. + for (int i = 0; i < 6; ++i) EXPECT_EQ(stress[i], 0.0); + for (int i = 0; i < 21; ++i) EXPECT_EQ(stiff[i], 0.0); +} + +/// A nonzero initial stress would make a preloaded model wrong from step 1. +TEST(CalculiXInterface, RefusesANonzeroInitialStress) { + FatalProbe probe; + const T props[2] = {K, G}; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; + const T beta[6] = {0, 0, 0, 0, 12.5, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + + clx_call c; + c.mprops = props; + c.nconst = 2; + c.beta = beta; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, stiff); + + EXPECT_EQ(FatalProbe::count, 1); + EXPECT_NE(FatalProbe::last.find("INITIAL"), std::string::npos); +} + +/// An all-zero beta is the normal case and must NOT trip the guard. +TEST(CalculiXInterface, AcceptsAZeroInitialStress) { + FatalProbe probe; + const T props[2] = {K, G}; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; + const T beta[6] = {0, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + + clx_call c; + c.mprops = props; + c.nconst = 2; + c.beta = beta; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, stiff); + + EXPECT_EQ(FatalProbe::count, 0); + EXPECT_GT(stress[0], 0.0); +} + +/// kode carries the constant count as -100 - nconst. Supplying fewer constants +/// than the model reads is a setup fault, and the decode is what surfaces it. +TEST(CalculiXInterface, DecodesTheConstantCountFromKode) { + FatalProbe probe; + const T props[2] = {K, G}; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + + clx_call c; + c.mprops = props; + c.nconst = 1; // kode = -101, but the model needs two constants + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, stiff); + + EXPECT_EQ(FatalProbe::count, 1); + EXPECT_NE(FatalProbe::last.find("constant"), std::string::npos); +} + +/// A symbol bound to a name nothing registered is unrecoverable. +TEST(CalculiXInterface, UnknownModelIsFatal) { + FatalProbe probe; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {5, 5, 5, 5, 5, 5}; + T stiff[21]; + for (int i = 0; i < 21; ++i) stiff[i] = 5.0; + + clx_call c; + c.run(&clx_missing_, "NOSUCHMODEL", emec, emec0, so, sn, stress, stiff); + + EXPECT_EQ(FatalProbe::count, 1); + for (int i = 0; i < 6; ++i) EXPECT_EQ(stress[i], 0.0); + // The zeroed 6x6 is packed through like any other result. + for (int i = 0; i < 21; ++i) EXPECT_EQ(stiff[i], 0.0); +} + +/// icmd == 3 requests stress only; stiff must be left as the caller supplied it. +TEST(CalculiXInterface, StressOnlyLeavesStiffUntouched) { + const T props[2] = {K, G}; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21]; + for (int i = 0; i < 21; ++i) stiff[i] = -999.0; + + clx_call c; + c.mprops = props; + c.nconst = 2; + c.icmd = 3; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, stiff); + + EXPECT_GT(stress[0], 0.0); + for (int i = 0; i < 21; ++i) + EXPECT_EQ(stiff[i], -999.0) << "stiff[" << i << "]"; +} + +} // namespace From e207412da8bc0f6ede54bec9eba208999b69bdbc Mon Sep 17 00:00:00 2001 From: petlenz Date: Mon, 14 Sep 2026 14:25:00 +0200 Subject: [PATCH 2/3] umat: cut the CalculiX interface comments down to the load-bearing facts --- .../umat/calculix_interface.h | 207 ++++++------------ tests/test_calculix_interface.cpp | 90 +++----- 2 files changed, 95 insertions(+), 202 deletions(-) diff --git a/include/numsim-materials/umat/calculix_interface.h b/include/numsim-materials/umat/calculix_interface.h index 743622e..f6c30b9 100644 --- a/include/numsim-materials/umat/calculix_interface.h +++ b/include/numsim-materials/umat/calculix_interface.h @@ -10,97 +10,58 @@ /// The CalculiX external-behaviour entry point (`call_external_umat_user`). /// -/// CalculiX has two external hooks, and they do NOT share a convention: +/// CalculiX has two external hooks. `call_external_umat` is reached from +/// umat_abaqus.f AFTER it has converted to the Abaqus convention, so +/// NUMSIM_MATERIALS_DEFINE_UMAT already serves it. `call_external_umat_user` is +/// the NATIVE hook, and this adapter translates it: /// -/// * `call_external_umat` is reached from umat_abaqus.f AFTER that wrapper has -/// already converted to the Abaqus UMAT convention (STRAN = strain at the -/// start with engineering shear, DSTRAN the increment, DDSDDE a full 6x6, -/// STATEV already sliced to the integration point). A material built with -/// NUMSIM_MATERIALS_DEFINE_UMAT already serves it — no adapter is needed. +/// kode #constants = -kode - 100 (umat_user.f:37-43) +/// emec/emec0 TENSORIAL {11,22,33,12,13,23}, end/start; doubled to the +/// engineering shear umat_dispatch expects (umat_abaqus.f:280) +/// stre PK2, in/out, unscaled +/// stiff(21) symmetrized upper triangle (umat_abaqus.f:335) +/// xstate* FULL (nstate_, mi(1), #elem) arrays, NOT a per-point slice, +/// so the callee indexes by iel/iint (umat_main.f:40,233) +/// TIME rebased onto the increment START (umat_abaqus.f:187) /// -/// * `call_external_umat_user` is the NATIVE hook (umat_user.f). It speaks -/// CalculiX's own convention, which this adapter translates to and from: +/// One emitted symbol serves one registered model: the deck picks a `@LIB,FUNC` +/// per material, so one FUNC implies one constant set. /// -/// - kode, not nprops: the 4th argument is `kode = -100 - #constants`, so -/// the constant count is `-kode - 100` (umat_user.f:37-43). -/// - emec / emec0 are the Lagrange mechanical strain at the END and START -/// of the increment, component order {11,22,33,12,13,23}, TENSORIAL -/// (no engineering doubling — umat_abaqus.f:280-283 doubles the shear -/// itself before handing strain to an Abaqus UMAT, which is exactly the -/// conversion this adapter performs). -/// - stre is the second Piola-Kirchhoff stress, in on entry / out on exit, -/// plain components with no scaling. -/// - stiff(21) is the upper triangle of the symmetric 6x6 tangent, packed -/// column-major, carrying no engineering factors. -/// - xstateini / xstate are the FULL three-dimensional state arrays -/// `(nstate_, mi(1), #elements)`, NOT a slice for this point. The callee -/// must index them itself — which is why iel, iint and mi are in the -/// signature. See the STATE VARIABLES note below. -/// - TIME must be rebased onto the START of the increment; ccx passes the -/// step time at its END. See the TIME note below. +/// SCOPE: geometrically linear. `emec` is Green-Lagrange and `stre` is PK2, +/// consumed as small-strain quantities. Exact for a linear `C:E` law (that IS +/// St-Venant-Kirchhoff); under NLGEOM an inelastic model would get a wrong +/// CONVERGED stress, not merely a slow tangent. /// -/// One emitted symbol maps to one registered model (the CalculiX deck selects a -/// specific `@LIB,FUNC` per material), so the model name is baked into the macro -/// rather than taken from `amat`. That also means one FUNC implies one constant -/// set: the registry's props-invariance check is what enforces it, and its -/// message is phrased for Abaqus ("use distinct *MATERIAL names"), which in ccx -/// means a distinct FUNC. -/// -/// SCOPE — GEOMETRICALLY LINEAR (SMALL STRAIN) ONLY. `emec` is the Green-Lagrange -/// strain and `stre` is PK2, while the models below consume and produce these as -/// small-strain quantities. For a linear `C:E` law the two coincide (that pairing -/// IS St-Venant-Kirchhoff), so the elastic case is exactly right. It is NOT -/// merely a slower tangent for anything else: under NLGEOM ccx passes large -/// Green-Lagrange strains, and feeding those to a small- or logarithmic-strain -/// return map produces a wrong CONVERGED stress, silently. Do not use this hook -/// for an inelastic model under finite deformation without first choosing a -/// conjugate stress-strain pair and converting here. -/// -/// INHERITED ERROR SEMANTICS. Fault handling belongs to umat_dispatch, not to -/// this file: a setup fault zeroes the outputs and terminates via the fatal -/// handler, and any other exception zeroes the outputs and requests a cutback -/// with PNEWDT = 0.25 — which is a valid ccx pnewdt (0 < pnewdt < 1), so the -/// Abaqus-side convention transfers unchanged. The zeroed 6x6 is packed into -/// stiff(21) like any other result. -/// -/// NOT handled, and rejected rather than ignored: a local material orientation -/// (`iorien != 0`) and a nonzero initial stress (`beta`). Both would otherwise be -/// silently wrong — see the guards in calculix_dispatch. Deliberately IGNORED: -/// the deformation gradients and Jacobians (small-strain scope, above), -/// temperature (constants are read once when the graph is built), and `ielas`, -/// whose elastic-iteration request an inelastic model would need to honour. +/// Errors belong to umat_dispatch: a setup fault zeroes outputs and terminates, +/// anything else zeroes them and asks for a cutback with PNEWDT = 0.25 (a valid +/// ccx pnewdt). REFUSED rather than ignored: `iorien != 0` and a nonzero `beta`. +/// IGNORED: deformation gradients, temperature, and `ielas`. namespace numsim::materials::umat { -/// The native hook is always full 3D: six components, and the symmetric 6x6 -/// tangent packed as 21 upper-triangular entries. +/// The native hook is always full 3D. inline constexpr std::size_t calculix_ntens = 6; inline constexpr std::size_t calculix_nstiff = calculix_ntens * (calculix_ntens + 1) / 2; -/// Every argument the CalculiX shim forwards, named. -/// -/// The Abaqus side funnels through dispatch_args for the same reason: a long -/// positional list of same-typed pointers lets a transposition compile cleanly, -/// and `time`/`ttime` or `emec`/`emec0` are exactly the pairs that would be -/// swapped. Scalars are held by value; CalculiX passes them by reference as -/// Fortran always does, and the macro dereferences them, so a null there would -/// be a CalculiX defect rather than a case to handle. +/// Every argument the shim forwards, named — a positional list of same-typed +/// pointers lets `time`/`ttime` or `emec`/`emec0` transpose silently. Scalars +/// are by value; CalculiX always passes them by reference, so the macro derefs. struct calculix_args { const char* amat{nullptr}; - int iel{0}; ///< element number, 1-based - int iint{0}; ///< integration point number, 1-based - int kode{0}; ///< -100 - #constants + int iel{0}; ///< 1-based + int iint{0}; ///< 1-based + int kode{0}; ///< -100 - #constants const double* mprops{nullptr}; - const double* emec{nullptr}; ///< tensorial strain at the END - const double* emec0{nullptr}; ///< tensorial strain at the START + const double* emec{nullptr}; ///< tensorial strain, END of increment + const double* emec0{nullptr}; ///< tensorial strain, START of increment const double* beta{nullptr}; ///< initial (residual) stress double dtime{0}; double time{0}; ///< STEP time at the END of the increment double ttime{0}; ///< TOTAL time at the START of the step - int icmd{0}; ///< 3 => stress only, leave stiff alone + int icmd{0}; ///< 3 => stress only int ielas{0}; - int mi1{0}; ///< mi(1): max integration points per element - int nstatv{0}; ///< nstate_: state variables PER integration point + int mi1{0}; ///< mi(1): max integration points per element + int nstatv{0}; ///< nstate_: state variables PER integration point const double* statev_old{nullptr}; ///< xstateini, FULL array double* statev_new{nullptr}; ///< xstate, FULL array double* stress{nullptr}; @@ -109,7 +70,6 @@ struct calculix_args { double* pnewdt{nullptr}; }; -/// Zero what CalculiX reads back, in CalculiX's own shapes. inline void calculix_zero_outputs(double* stress, double* stiff) noexcept { if (stress) for (std::size_t i = 0; i < calculix_ntens; ++i) stress[i] = 0.0; @@ -119,22 +79,18 @@ inline void calculix_zero_outputs(double* stress, double* stiff) noexcept { /// Shared implementation behind the CalculiX external symbol. /// -/// Everything here is fixed-size stack work with no allocation, so it cannot -/// throw before reaching umat_dispatch, which owns the try/catch that keeps a -/// C++ exception from unwinding into CalculiX's Fortran. The guards below -/// report through the same fatal handler rather than throwing, for that reason. +/// Fixed-size stack work only, so nothing here can throw before umat_dispatch, +/// which owns the try/catch keeping exceptions out of Fortran. The guards report +/// through the fatal handler rather than throwing, for the same reason. template void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept { using T = typename Traits::value_type; static_assert(std::is_same_v, "the CalculiX external interface is double-precision (ccxreal)"); - // A local orientation is NOT applied around this hook. umat_abaqus.f rotates - // strain in and stress/stiffness out around its call (lines 202-276, 305-330), - // but umat_main.f does no such thing around the native or external umat_user - // call: umat_user.f:86-104 requires results in the MATERIAL frame and tells - // the user to rotate with transformatrix(). Ignoring iorien would therefore - // return results in the wrong frame with no symptom, so refuse instead. + // umat_main.f does not rotate around this hook, unlike umat_abaqus.f; the + // material owes results in the material frame (umat_user.f:86-104). Ignoring + // iorien would return the wrong frame with no symptom. if (a.iorien != 0) { calculix_zero_outputs(a.stress, a.stiff); report_fatal(model_name, @@ -145,8 +101,8 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept return; } - // *INITIAL CONDITIONS,TYPE=STRESS (umat_user.f:52). Dropping it would make a - // preloaded model wrong from the first step, again with no symptom. + // *INITIAL CONDITIONS,TYPE=STRESS: dropping it makes a preloaded model wrong + // from the first step. if (a.beta) { for (std::size_t i = 0; i < calculix_ntens; ++i) { if (a.beta[i] != 0.0) { @@ -167,18 +123,14 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept return; } - // Widened before negating so kode == INT_MIN cannot trap, and clamped because - // a malformed kode must not become a negative count. + // Widened before negating so kode == INT_MIN cannot trap; clamped so a + // malformed kode cannot become a negative count. const long long nconst_ll = -static_cast(a.kode) - 100; const int nconst = nconst_ll > 0 ? static_cast(nconst_ll) : 0; const int nstatv = a.nstatv > 0 ? a.nstatv : 0; - // STATE VARIABLES. xstateini/xstate are the FULL arrays - // `real*8 xstate(nstate_, mi(1), *)` (umat_main.f:40), and umat_main.f passes - // their BASE to this hook (line 233) — unlike umat_abaqus.f:295, which slices - // `xstate(1,iint,iel)` before calling. So the offset of this point's block is - // the Fortran column-major stride. Getting this wrong does not fail loudly: it - // makes every integration point share element 1 / point 1's state. + // Fortran column-major stride into the full state array. Getting this wrong is + // silent: every point would share element 1 / point 1's state. std::size_t point_offset = 0; if (nstatv > 0) { if (a.iel < 1 || a.iint < 1 || a.mi1 < 1 || a.iint > a.mi1) { @@ -195,11 +147,8 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept static_cast(a.iel - 1)); } - // Native emec is TENSORIAL; umat_dispatch consumes ENGINEERING shear - // (strain_from_buffer halves slots 3-5). Doubling the shear here makes the - // round-trip land back on the tensorial strain the models expect — the same - // conversion umat_abaqus.f:280-283 applies before an Abaqus UMAT. STRAN is the - // strain at the START (emec0); DSTRAN the increment, so STRAN+DSTRAN == emec. + // Tensorial -> engineering, so strain_from_buffer's halving lands back on the + // tensorial strain. STRAN is the start, so STRAN + DSTRAN == emec. double stran[calculix_ntens]; double dstran[calculix_ntens]; for (std::size_t i = 0; i < calculix_ntens; ++i) { @@ -208,10 +157,8 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept dstran[i] = shear * (a.emec[i] - a.emec0[i]); } - // The evaluator reads and writes ONE state array in place; CalculiX splits it - // into a read-only xstateini and a write-only xstate. Seed this point's block - // in the output from the committed input so the read sees last increment's - // state and the write lands where CalculiX expects the update. + // The evaluator updates one array in place; CalculiX splits read and write, so + // seed this point's block from the committed state. double* statev_point = nullptr; if (nstatv > 0 && a.statev_new) { statev_point = a.statev_new + point_offset; @@ -219,21 +166,12 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept std::copy_n(a.statev_old + point_offset, nstatv, statev_point); } - // umat_dispatch fills a full 6x6; CalculiX wants only the packed upper - // triangle. Note this buffer is COLUMN-major (material_point_evaluator.h:82, - // narrow_matrix writes host[a + b*n]), which the packing below relies on. + // COLUMN-major (material_point_evaluator.h:82), which the packing below needs. double ddsdde36[calculix_ntens * calculix_ntens] = {0.0}; - // TIME must be rebased onto the START of the increment. ccx hands over the - // step time at the END of the increment plus the total time at the start of - // the STEP, so mirror umat_abaqus.f:187-188 exactly: - // abqtime(1) = time - dtime - // abqtime(2) = ttime + time - dtime - // Passing {time, ttime} instead is correct only on the first increment of the - // first step, and silently wrong for every rate- or time-dependent model after - // that. - const double time2[2] = {a.time - a.dtime, - a.ttime + a.time - a.dtime}; + // ccx gives step time at the increment END and total time at the STEP start. + // Passing {time, ttime} through is right only on the very first increment. + const double time2[2] = {a.time - a.dtime, a.ttime + a.time - a.dtime}; dispatch_args d; d.stress = a.stress; @@ -255,18 +193,10 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept umat_dispatch(d); - // Pack into stiff(21): column-major upper triangle, so the 0-based (i, j) with - // i <= j lands at i + j*(j+1)/2. - // - // SYMMETRIZED, exactly as the reference umat_abaqus.f:335-355 does - // (stiff(2) = (ddsdde(1,2)+ddsdde(2,1))/2, and so on). Two reasons: ddsdde36 - // is column-major, so reading it row-major would transpose every off-diagonal - // for a major-ASYMMETRIC tangent (non-associative flow, damage); and stiff(21) - // has no room for an antisymmetric part anyway, so averaging is what CalculiX - // itself keeps. For a symmetric tangent this is the identity. - // - // icmd == 3 requests stress only and leaves stiff untouched, as umat_abaqus.f - // does; on an error path umat_dispatch has zeroed ddsdde36 already. + // Column-major upper triangle. SYMMETRIZED as umat_abaqus.f:335-355 does: a + // row-major read would transpose every off-diagonal of an asymmetric tangent, + // and stiff(21) has no room for an antisymmetric part anyway. icmd == 3 leaves + // stiff alone; on an error path ddsdde36 is already zeroed. if (a.stiff && a.icmd != 3) { for (std::size_t j = 0; j < calculix_ntens; ++j) for (std::size_t i = 0; i <= j; ++i) @@ -278,21 +208,14 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept } // namespace numsim::materials::umat -/// Emit a CalculiX-callable external-behaviour symbol @p FUNC bound to the -/// registered model @p MODELNAME. Place this in exactly ONE translation unit of -/// the shared library CalculiX dlopen's for `*MATERIAL, NAME=@LIB,FUNC`. -/// -/// The signature is the `calculixptr` prototype from CalculiX's -/// call_external_umat_user.c: every argument is passed by pointer (ccxint == -/// int, ccxreal == double), and the final `int` is the hidden Fortran length of -/// the leading character string. That trailing width is `int` because the C -/// typedef in call_external_umat_user.c fixes it as `int` — unlike the Abaqus -/// entry, where the Fortran compiler chooses and -/// NUMSIM_MATERIALS_FORTRAN_STRLEN exists to track it. Do not make it -/// configurable here; the C prototype is the contract. +/// Emit a CalculiX-callable symbol @p FUNC bound to model @p MODELNAME. Place in +/// exactly ONE translation unit of the library ccx loads for `@LIB,FUNC`. /// -/// Arguments the small-strain path does not use are named in comments so their -/// absence reads as deliberate, not forgotten. +/// Signature is the `calculixptr` prototype (call_external_umat_user.c): all +/// arguments by pointer (ccxint == int, ccxreal == double), with the hidden +/// Fortran string length last. That length is `int` because the C typedef fixes +/// it — unlike the Abaqus entry, where the Fortran compiler chooses. Do not make +/// it configurable. #define NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR(TRAITS, FUNC, MODELNAME) \ extern "C" void FUNC( \ const char* AMAT, const int* IEL, const int* IINT, \ diff --git a/tests/test_calculix_interface.cpp b/tests/test_calculix_interface.cpp index da3daf8..8f9833b 100644 --- a/tests/test_calculix_interface.cpp +++ b/tests/test_calculix_interface.cpp @@ -17,9 +17,7 @@ #include "numsim-materials/umat/external_state_source.h" #include "numsim-materials/umat/umat_interface.h" -// Both entry points live in this one translation unit: umat_ is the reference -// path the CalculiX adapter is checked against, and the clx_* symbols are the -// adapter under test. They are distinct extern "C" symbols, so they coexist. +// umat_ is the reference path; the clx_* symbols are the adapter under test. NUMSIM_MATERIALS_DEFINE_UMAT(numsim::materials::material_policy_default) NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( numsim::materials::material_policy_default, clx_linear_elastic_, "LINELAS") @@ -64,8 +62,7 @@ void build_deck_elastic(ctx_type& ctx, std::span props) { ctx.finalize(); } -/// A genuinely stateful model, so the xstateini -> xstate plumbing has real -/// state to carry. Mirrors the sibling suite's J2 build. +/// A stateful model, so the xstateini -> xstate plumbing has real state. void build_j2(ctx_type& ctx, std::span /*props*/) { param_type p; p.insert("name", "stepper"); @@ -100,18 +97,13 @@ void build_j2(ctx_type& ctx, std::span /*props*/) { constexpr int slot_i[6] = {0, 1, 2, 0, 0, 1}; constexpr int slot_j[6] = {0, 1, 2, 1, 2, 2}; -/// The value this probe puts at 6x6 entry (a, b). Deliberately NOT symmetric: -/// asym_value(a,b) != asym_value(b,a) for a != b, which is what makes a -/// transposed or one-sided read of the tangent observable. +/// Deliberately NOT symmetric, so a transposed or one-sided read is visible. constexpr T asym_value(int a, int b) { return static_cast(100 * (a + 1) + (b + 1)); } -/// Publishes a known, minor-symmetric but major-ASYMMETRIC tangent. -/// -/// Minor symmetry is required (tangent_to_buffer asserts it); major symmetry is -/// not, and is exactly what a non-associative model lacks. Stress is left zero: -/// this probe exists only to pin the tangent packing. +/// Minor-symmetric (which tangent_to_buffer asserts) but major-ASYMMETRIC, as a +/// non-associative model is. Stress stays zero; this pins the tangent packing. template class asym_tangent_material final : public nm::material_base, Traits> { @@ -159,9 +151,8 @@ class asym_tangent_material final tensor4& m_tangent; }; -/// Reports the time it was bound with, so the TIME rebasing is observable from -/// outside. Without a probe like this, nothing reaching a material through the -/// shim depends on time at all and a wrong TIME(2) is undetectable at the ABI. +/// Reports the time it was bound with; nothing else through the shim depends on +/// time, so without this a wrong TIME(2) is undetectable at the ABI. template class time_probe_material final : public nm::material_base, Traits> { @@ -282,8 +273,7 @@ using clx_fn = void (*)(const char*, const int*, const int*, const int*, double*, double*, double*, const int*, const double*, const double*, double*, const int*, const int); -/// Everything a CalculiX call needs, with the defaults a well-formed solid-3D -/// call uses. Only what a test actually varies is set. +/// Defaults for a well-formed solid-3D call; tests set only what they vary. struct clx_call { int iel = 1; int iint = 1; @@ -353,12 +343,9 @@ void call_umat_reference(const std::string& name, T* stress, T* statev, // Convention translation // --------------------------------------------------------------------------- -/// The load-bearing claim: the CalculiX adapter and the Abaqus entry, given the -/// SAME physical strain expressed in their respective conventions, produce the -/// same stress and the same tangent. A nonzero shear component makes the -/// tensorial->engineering doubling observable, and a NONZERO start strain makes -/// the stran/dstran split observable — with emec0 = 0 a swapped or dropped -/// emec0 is invisible. +/// The same physical strain, in each hook's own convention, must give the same +/// stress and tangent. Nonzero shear exposes the engineering doubling; nonzero +/// emec0 exposes the stran/dstran split (at emec0 = 0 a swap is invisible). TEST(CalculiXInterface, MatchesTheAbaqusEntryForTheSamePhysicalStrain) { const T props[2] = {K, G}; const int nconst = 2; @@ -404,8 +391,7 @@ TEST(CalculiXInterface, MatchesTheAbaqusEntryForTheSamePhysicalStrain) { << "stiff(" << i << "," << j << ")"; } -/// The shear conversion is not a no-op: feeding emec directly (tensorial) where -/// the adapter must double it would halve every shear-driven stress. +/// Feeding emec straight through would halve every shear-driven stress. TEST(CalculiXInterface, TensorialShearIsConvertedToEngineering) { const T props[2] = {K, G}; const T e12 = 2.5e-3; @@ -431,11 +417,9 @@ TEST(CalculiXInterface, TensorialShearIsConvertedToEngineering) { EXPECT_NEAR(stiff[3 + 3 * (3 + 1) / 2], G, 1e-9); } -/// stiff(21) must carry the SYMMETRIZED tangent, as umat_abaqus.f:335-355 does. -/// -/// The buffer umat_dispatch fills is column-major, so a row-major read would -/// hand CalculiX C(j,i) instead of the average — invisible for every symmetric -/// material, which is why this probe is deliberately major-asymmetric. +/// stiff(21) carries the SYMMETRIZED tangent (umat_abaqus.f:335-355). The buffer +/// is column-major, so a row-major read would hand ccx C(j,i) — invisible for +/// any symmetric material, hence the asymmetric probe. TEST(CalculiXInterface, PacksTheSymmetrizedTangentIntoStiff21) { const T emec0[6] = {0, 0, 0, 0, 0, 0}; const T emec[6] = {1.0e-3, 0, 0, 0, 0, 0}; @@ -453,8 +437,7 @@ TEST(CalculiXInterface, PacksTheSymmetrizedTangentIntoStiff21) { EXPECT_NEAR(stiff[i + j * (j + 1) / 2], expected, 1e-9) << "stiff(" << i << "," << j << ")"; if (i != j) { - // Pin the bug this test exists for: a one-sided read would land on - // exactly one of these, and the average equals neither. + // A one-sided read lands on one of these; the average equals neither. EXPECT_NE(expected, asym_value(i, j)); EXPECT_NE(expected, asym_value(j, i)); } @@ -466,16 +449,11 @@ TEST(CalculiXInterface, PacksTheSymmetrizedTangentIntoStiff21) { // State variables: the full-array indexing contract // --------------------------------------------------------------------------- -/// CalculiX passes the WHOLE xstateini/xstate arrays — `(nstate_, mi(1), -/// #elements)`, umat_main.f:40 — and their base pointer, not a slice for this -/// point (contrast umat_abaqus.f:295). The adapter must therefore index by -/// (iint, iel) itself. -/// -/// Driven at element 2, integration point 3, this fails outright if the adapter -/// works at offset 0: the trajectory would be read back from, and written to, -/// the wrong block. The untouched blocks are checked too, because writing the -/// right values into the wrong place AND leaving neighbours dirty are different -/// bugs. +/// ccx passes the WHOLE state arrays' base (umat_main.f:40,233), not a per-point +/// slice as umat_abaqus.f:295 does, so the adapter must index by (iint, iel). +/// Driven at element 2, point 3, this fails outright at offset 0. Neighbouring +/// blocks are checked too: right values in the wrong place and dirty neighbours +/// are different bugs. TEST(CalculiXInterface, IndexesStateByElementAndIntegrationPoint) { const int nstatv = static_cast(registry::instance().nstatv("J2CLX")); ASSERT_GT(nstatv, 0) << "this test is meaningless without real state"; @@ -488,8 +466,7 @@ TEST(CalculiXInterface, IndexesStateByElementAndIntegrationPoint) { const std::size_t offset = static_cast(nstatv) * ((iint - 1) + mi1 * (iel - 1)); - // A strain ramp well past yield, so state genuinely accumulates. TENSORIAL, - // so the shear slot is half the engineering increment the sibling suite uses. + // Past yield, so state accumulates. TENSORIAL: shear is half the engineering. const T de[6] = {0.01, -0.0025, 0.0, 0.0025, 0.0, 0.0}; constexpr int steps = 40; @@ -512,7 +489,7 @@ TEST(CalculiXInterface, IndexesStateByElementAndIntegrationPoint) { c.nstatv = nstatv; c.run(&clx_j2_, "J2CLX", emec, emec0, st_old.data(), st_new.data(), clx_stress, clx_stiff); - // CalculiX commits xstate -> xstateini between increments. + // ccx commits xstate -> xstateini between increments. st_old = st_new; } @@ -559,9 +536,8 @@ TEST(CalculiXInterface, IndexesStateByElementAndIntegrationPoint) { /// the START of the increment (umat_abaqus.f:187-188). Passing {time, ttime} /// through is right only on the first increment of the first step. /// -/// The probe reports the bound times directly, so the exact rebasing is -/// asserted rather than inferred. The material sees the TOTAL time spanning the -/// increment: [ttime + time - dtime, ttime + time]. +/// The probe reports the bound times, so the rebasing is asserted exactly: the +/// material sees [ttime + time - dtime, ttime + time]. TEST(CalculiXInterface, RebasesTimeOntoTheStartOfTheIncrement) { const T emec0[6] = {0, 0, 0, 0, 0, 0}; const T emec[6] = {1.0e-4, 0, 0, 0, 0, 0}; @@ -579,8 +555,7 @@ TEST(CalculiXInterface, RebasesTimeOntoTheStartOfTheIncrement) { return std::array{stress[0], stress[1], stress[2]}; }; - // First increment of the first step: time = dtime, ttime = 0. This is the one - // case where passing {time, ttime} straight through would also be correct. + // First increment of the first step — the one case the naive mapping gets right. { const auto t = times_seen(/*time=*/0.1, /*ttime=*/0.0, /*dt=*/0.1); EXPECT_NEAR(t[0], 0.0, 1e-12) << "total time at the START of the increment"; @@ -588,9 +563,7 @@ TEST(CalculiXInterface, RebasesTimeOntoTheStartOfTheIncrement) { EXPECT_NEAR(t[2], 0.1, 1e-12) << "dtime"; } - // A later increment of a later step — where the naive mapping is wrong. - // ccx: step time ends at 0.9, the step began at total time 5.0, dt = 0.1. - // So the increment spans total time [5.0 + 0.9 - 0.1, 5.0 + 0.9] = [5.8, 5.9]. + // A later increment of a later step: spans [5.0+0.9-0.1, 5.0+0.9] = [5.8, 5.9]. { const auto t = times_seen(/*time=*/0.9, /*ttime=*/5.0, /*dt=*/0.1); EXPECT_NEAR(t[0], 5.8, 1e-12) @@ -604,8 +577,7 @@ TEST(CalculiXInterface, RebasesTimeOntoTheStartOfTheIncrement) { // Guards and error paths // --------------------------------------------------------------------------- -/// A local orientation is not applied around the native hook, so silently -/// ignoring it would return results in the wrong frame. It must be refused. +/// Nothing rotates around the native hook, so ignoring iorien is wrong-frame. TEST(CalculiXInterface, RefusesALocalOrientation) { FatalProbe probe; const T props[2] = {K, G}; @@ -624,8 +596,7 @@ TEST(CalculiXInterface, RefusesALocalOrientation) { EXPECT_EQ(FatalProbe::count, 1); EXPECT_NE(FatalProbe::last.find("ORIENTATION"), std::string::npos); - // Outputs are zeroed, so a handler that returns cannot leave the solver - // consuming the caller's uninitialised buffer. + // Zeroed, so a returning handler cannot leave the solver on a stale buffer. for (int i = 0; i < 6; ++i) EXPECT_EQ(stress[i], 0.0); for (int i = 0; i < 21; ++i) EXPECT_EQ(stiff[i], 0.0); } @@ -672,8 +643,7 @@ TEST(CalculiXInterface, AcceptsAZeroInitialStress) { EXPECT_GT(stress[0], 0.0); } -/// kode carries the constant count as -100 - nconst. Supplying fewer constants -/// than the model reads is a setup fault, and the decode is what surfaces it. +/// kode carries the count as -100 - nconst; too few constants is a setup fault. TEST(CalculiXInterface, DecodesTheConstantCountFromKode) { FatalProbe probe; const T props[2] = {K, G}; From 1f67c66057d63522fe879cc635a6299a73f2bd90 Mon Sep 17 00:00:00 2001 From: petlenz Date: Mon, 14 Sep 2026 18:17:24 +0200 Subject: [PATCH 3/3] umat: close the CalculiX review's open items Four gaps found by re-reviewing this branch against the real ccx 2.22 source. The conversions themselves held up -- all eight rows of the table, the 31-arg ABI, the stiff(21) packing and the TIME rebasing were re-derived from the Fortran and are correct -- so these are the edges. ielas is now refused rather than ignored. ielas = 1 asks for a response with no irreversible deformation; arpack.c, arpackbu.c and arpackcs.c all set it, which is every *BUCKLE and *FREQUENCY step. Nothing here can suppress plastic flow on request, so answering anyway folds irreversible effects into the tangent and skews the extracted eigenvalues with no symptom. It joins iorien and beta in the refuse-rather-than-guess group. The emec/emec0 null guard had no test. REVIEW.md's resolution table listed M6 as fixed, and deleting the guard left all eleven tests passing -- a claim, not a behaviour. With the new test the binary segfaults instead, which is the guard's whole point. A fatal tripped outside a FatalProbe aborted the entire test binary. The default handler calls std::abort(), which is right in a solver -- returning would hand the host a silently wrong material response -- and wrong in a test run, where it masks every other result. Mutating the kode decode by one was enough to lose the whole suite to a SIGABRT. A binary-wide handler now fails the current test instead; verified by mutation to give exit 1 with the other 15 tests still running. FatalProbe restores that handler rather than the library default, so the first probe in a run no longer re-arms abort() for everything after it. The header said nothing about how the .so is found. ccx prepends "lib", appends ".so" and dlsyms FUNC verbatim, which is why FUNC carries no trailing underscore here while the Abaqus entry point is umat_ -- get it backwards and it fails at run time with "unable to load function" and nothing at build time. Also recorded: the symbol must stay exported (-fvisibility=hidden hides it with no other symptom), and ccx builds the library name into a fixed char b[80] with an unchecked memcpy, so the LIB name must be short. Two tests written during the review are ported in, closing the last two items REVIEW.md listed as still open: a throwing material proves PNEWDT = 0.25 lands on the CALLER's buffer for a cutback, and eight threads driving the entry point at distinct (iel, iint) must agree with the same sequences run serially. Every one of the four is mutation-verified: removing the ielas guard, the null guard, the pnewdt forwarding, or the isolation handler each fails exactly its own test. 16/16 pass. --- .../umat/calculix_interface.h | 49 ++++- tests/test_calculix_interface.cpp | 206 +++++++++++++++++- 2 files changed, 252 insertions(+), 3 deletions(-) diff --git a/include/numsim-materials/umat/calculix_interface.h b/include/numsim-materials/umat/calculix_interface.h index f6c30b9..48ae1d3 100644 --- a/include/numsim-materials/umat/calculix_interface.h +++ b/include/numsim-materials/umat/calculix_interface.h @@ -34,8 +34,37 @@ /// /// Errors belong to umat_dispatch: a setup fault zeroes outputs and terminates, /// anything else zeroes them and asks for a cutback with PNEWDT = 0.25 (a valid -/// ccx pnewdt). REFUSED rather than ignored: `iorien != 0` and a nonzero `beta`. -/// IGNORED: deformation gradients, temperature, and `ielas`. +/// ccx pnewdt). REFUSED rather than ignored: `iorien != 0`, a nonzero `beta`, +/// and `ielas != 0`. IGNORED: deformation gradients and temperature. +/// +/// ## Building and naming the .so +/// +/// ccx resolves the plugin itself, and both halves of the name are its choice, +/// not ours. `external.c` prepends "lib" and appends ".so" to the LIB part of +/// `*MATERIAL, NAME=@LIB,FUNC`, then `dlsym`s FUNC **verbatim** -- no Fortran +/// mangling, no trailing underscore. So for +/// +/// *MATERIAL, NAME=@numsimmat,my_model +/// +/// the shared object must be `libnumsimmat.so` on ccx's library path, and this +/// macro must be instantiated with FUNC spelled exactly `my_model`: +/// +/// NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR(my_traits, my_model, "MYMODEL") +/// +/// That is why FUNC carries no underscore here while the Abaqus entry point in +/// umat_interface.h is `umat_` -- that one is called from Fortran directly, this +/// one through dlsym. Getting it backwards fails at run time with ccx's +/// "unable to load function" and produces nothing at build time. +/// +/// Two further constraints, both ccx's: +/// +/// * The symbol must be EXPORTED. `extern "C"` alone is enough only while the +/// build leaves default visibility; a target compiled `-fvisibility=hidden` +/// hides it and ccx reports "unable to load function" with no other symptom. +/// Add `-fvisibility=default` for this translation unit if that changes. +/// * Keep LIB short. `external.c` builds the file name into a fixed `char +/// b[80]` with an unchecked `memcpy`, so a long name overwrites its stack. +/// namespace numsim::materials::umat { /// The native hook is always full 3D. @@ -91,6 +120,22 @@ void calculix_dispatch(const calculix_args& a, const char* model_name) noexcept // umat_main.f does not rotate around this hook, unlike umat_abaqus.f; the // material owes results in the material frame (umat_user.f:86-104). Ignoring // iorien would return the wrong frame with no symptom. + // ielas = 1 is an ELASTIC iteration: ccx is asking for a response with no + // irreversible deformation. arpack.c, arpackbu.c and arpackcs.c set it, i.e. + // every eigenvalue and buckling analysis. A material that ignores it returns + // a tangent with plastic flow folded in, and the extracted eigenvalues are + // quietly wrong. Nothing in this library can suppress irreversible effects on + // request, so refuse rather than answer the wrong question. + if (a.ielas != 0) { + calculix_zero_outputs(a.stress, a.stiff); + report_fatal(model_name, + "numsim CalculiX: an elastic iteration (ielas != 0) was " + "requested -- *BUCKLE and *FREQUENCY need a response with no " + "irreversible deformation, which this material cannot " + "produce. Use a linear elastic material for those steps."); + return; + } + if (a.iorien != 0) { calculix_zero_outputs(a.stress, a.stiff); report_fatal(model_name, diff --git a/tests/test_calculix_interface.cpp b/tests/test_calculix_interface.cpp index 8f9833b..a08390d 100644 --- a/tests/test_calculix_interface.cpp +++ b/tests/test_calculix_interface.cpp @@ -1,9 +1,12 @@ #include #include +#include #include #include +#include #include +#include #include #include @@ -29,6 +32,8 @@ NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( numsim::materials::material_policy_default, clx_missing_, "NOSUCHMODEL") NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( numsim::materials::material_policy_default, clx_time_, "TIMEPROBE") +NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( + numsim::materials::material_policy_default, clx_throws_, "THROWSCLX") namespace { @@ -187,6 +192,42 @@ class time_probe_material final const nm::input_history& m_time; }; +/// Always throws a plain exception, to drive the cutback path. +template +class throwing_material final + : public nm::material_base, Traits> { + public: + using base = nm::material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + using tensor2 = tmech::tensor; + using tensor4 = tmech::tensor; + + template + explicit throwing_material(Args&&... args) + : base(std::forward(args)...), + m_stress(base::template add_output( + "stress", &throwing_material::compute)), + m_tangent(base::template add_output("tangent")) {} + + static input_parameter_controller parameters() { return base::parameters(); } + void compute() { throw std::runtime_error("deliberate non-convergence"); } + + private: + tensor2& m_stress; + tensor4& m_tangent; +}; + +void build_throws(ctx_type& ctx, std::span /*props*/) { + param_type p; + p.insert("name", "stepper"); + ctx.create>(p); + p.clear(); + p.insert("name", "probe"); + ctx.create>(p); + ctx.finalize(); +} + void build_time_probe(ctx_type& ctx, std::span /*props*/) { param_type p; p.insert("name", "stepper"); @@ -234,6 +275,11 @@ struct Registration { tp.stress_source = "probe"; tp.time_source = "clock"; registry::instance().register_model("TIMEPROBE", build_time_probe, tp); + + registry::config th; + th.strain_source = "stepper"; + th.stress_source = "probe"; + registry::instance().register_model("THROWSCLX", build_throws, th); } }; const Registration registration; @@ -249,6 +295,20 @@ struct fortran_name { /// A fatal fault must terminate the analysis. The handler is replaced so the /// test can observe it instead of the runner being killed by XIT/abort. +/// The default fatal handler calls std::abort(), which is right in a solver -- +/// returning would hand the host a silently wrong material response -- and +/// wrong in a test binary, where it takes every remaining test with it. This +/// installs a binary-wide handler that fails the current test instead. A test +/// that EXPECTS a fatal wraps itself in FatalProbe. +struct UnexpectedFatalsFailTheTest : ::testing::Environment { + static void handler(const char* msg) { + ADD_FAILURE() << "unexpected fatal from the adapter: " << msg; + } + void SetUp() override { u::set_fatal_handler(&handler); } +}; +const auto* const fatal_env = + ::testing::AddGlobalTestEnvironment(new UnexpectedFatalsFailTheTest); + struct FatalProbe { static inline std::string last; static inline int count = 0; @@ -261,7 +321,9 @@ struct FatalProbe { count = 0; u::set_fatal_handler(&handler); } - ~FatalProbe() { u::set_fatal_handler(nullptr); } + /// Restore the binary-wide handler, NOT the library default -- otherwise the + /// first probe in a run re-arms std::abort() for every test after it. + ~FatalProbe() { u::set_fatal_handler(&UnexpectedFatalsFailTheTest::handler); } }; using clx_fn = void (*)(const char*, const int*, const int*, const int*, @@ -663,6 +725,148 @@ TEST(CalculiXInterface, DecodesTheConstantCountFromKode) { } /// A symbol bound to a name nothing registered is unrecoverable. +/// M6's guard, which the review table lists as resolved and nothing exercised. +/// Removing the null check leaves every other test passing; with this one, the +/// binary segfaults instead, which is the point -- the guard exists so a wiring +/// mistake degrades rather than crashing. +TEST(CalculiXInterface, RefusesANullStrainPointer) { + const T props[2] = {166.67, 76.92}; + const T emec[6] = {1e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + + for (int which = 0; which < 2; ++which) { + T stress[6] = {9, 9, 9, 9, 9, 9}; + T stiff[21]; + for (auto& v : stiff) v = 5.0; + + FatalProbe probe; + clx_call c; + c.mprops = props; + c.nconst = 2; + c.run(&clx_linear_elastic_, "LINELAS", which == 0 ? nullptr : emec, + which == 0 ? emec : nullptr, so, sn, stress, stiff); + + EXPECT_EQ(FatalProbe::count, 1) + << (which == 0 ? "null emec" : "null emec0") << " was not refused"; + for (int i = 0; i < 6; ++i) EXPECT_EQ(stress[i], 0.0) << "i=" << i; + for (int i = 0; i < 21; ++i) EXPECT_EQ(stiff[i], 0.0) << "i=" << i; + } +} + +/// ielas = 1 asks for a response with no irreversible deformation -- what +/// *BUCKLE and *FREQUENCY need; arpack.c, arpackbu.c and arpackcs.c all set it. +/// Nothing here can suppress plastic flow on request, so answering anyway would +/// fold irreversible effects into the tangent and quietly skew the eigenvalues. +TEST(CalculiXInterface, RefusesAnElasticIteration) { + const T props[2] = {166.67, 76.92}; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {9, 9, 9, 9, 9, 9}; + T stiff[21]; + for (auto& v : stiff) v = 5.0; + + FatalProbe probe; + clx_call c; + c.mprops = props; + c.nconst = 2; + c.ielas = 1; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, stiff); + + EXPECT_EQ(FatalProbe::count, 1) << "an elastic iteration must be refused"; + EXPECT_NE(FatalProbe::last.find("ielas"), std::string::npos) + << FatalProbe::last; + for (int i = 0; i < 6; ++i) EXPECT_EQ(stress[i], 0.0) << "i=" << i; + for (int i = 0; i < 21; ++i) EXPECT_EQ(stiff[i], 0.0) << "i=" << i; +} + +/// And ielas = 0 -- every static increment -- must still go through. +TEST(CalculiXInterface, AcceptsTheOrdinaryNonElasticIteration) { + const T props[2] = {166.67, 76.92}; + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + + clx_call c; + c.mprops = props; + c.nconst = 2; + c.ielas = 0; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, stiff); + EXPECT_GT(std::abs(stress[0]), 1e-12) << "the ordinary path was refused too"; +} + +/// A cutback must reach ccx through PNEWDT on the CALLER's buffer. umat_user.f: +/// pnewdt "should exceed zero but be less than 1. Default is -1 indicating that +/// the user routine has converged", and checkconvergence.c multiplies the step +/// by it. Drop the write and a diverged point is reported as converged. +TEST(CalculiXInterface, ACutbackPropagatesThroughTheAdapter) { + const T emec0[6] = {0, 0, 0, 0, 0, 0}; + const T emec[6] = {1e-3, 0, 0, 0, 0, 0}; + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {1, 1, 1, 1, 1, 1}; + T stiff[21]; + for (auto& v : stiff) v = 1.0; + + clx_call c; + c.pnewdt = -1.0; // ccx's "converged" default + c.run(&clx_throws_, "THROWSCLX", emec, emec0, so, sn, stress, stiff); + + EXPECT_DOUBLE_EQ(c.pnewdt, 0.25) + << "a throwing material must ask ccx for a smaller increment"; + for (int i = 0; i < 6; ++i) EXPECT_EQ(stress[i], 0.0) << "i=" << i; + for (int i = 0; i < 21; ++i) EXPECT_EQ(stiff[i], 0.0) << "i=" << i; +} + +/// ccx drives one .so from a threaded element loop, so the entry point is +/// re-entered concurrently at distinct (iel, iint). umat_interface caches one +/// context per thread; this pins that the adapter adds no sharing of its own. +TEST(CalculiXInterface, ConcurrentCallsAtDistinctPointsAgreeWithSerialOnes) { + constexpr int nthreads = 8; + constexpr int steps = 50; + const T props[2] = {166.67, 76.92}; + + auto drive = [&](int thread_index) { + T emec0[6] = {0, 0, 0, 0, 0, 0}, emec[6] = {0, 0, 0, 0, 0, 0}; + T last = 0; + for (int step = 0; step < steps; ++step) { + for (int i = 0; i < 6; ++i) { + emec0[i] = emec[i]; + emec[i] = emec0[i] + 1e-5 * (thread_index + 1); + } + T so[1] = {0}, sn[1] = {0}; + T stress[6] = {0, 0, 0, 0, 0, 0}; + T stiff[21] = {0}; + clx_call c; + c.iel = thread_index + 1; + c.iint = 1; + c.mi1 = 4; + c.mprops = props; + c.nconst = 2; + c.run(&clx_linear_elastic_, "LINELAS", emec, emec0, so, sn, stress, + stiff); + last = stress[0]; + } + return last; + }; + + std::vector serial(nthreads); + for (int t = 0; t < nthreads; ++t) serial[t] = drive(t); + + std::vector concurrent(nthreads); + std::vector pool; + for (int t = 0; t < nthreads; ++t) + pool.emplace_back([&, t] { concurrent[t] = drive(t); }); + for (auto& th : pool) th.join(); + + for (int t = 0; t < nthreads; ++t) { + ASSERT_TRUE(std::isfinite(serial[t])) << "t=" << t; + EXPECT_DOUBLE_EQ(concurrent[t], serial[t]) + << "thread " << t << " disagreed with the same sequence run serially"; + } +} + TEST(CalculiXInterface, UnknownModelIsFatal) { FatalProbe probe; const T emec0[6] = {0, 0, 0, 0, 0, 0};