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..48ae1d3 --- /dev/null +++ b/include/numsim-materials/umat/calculix_interface.h @@ -0,0 +1,289 @@ +#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. `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: +/// +/// 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) +/// +/// One emitted symbol serves one registered model: the deck picks a `@LIB,FUNC` +/// per material, so one FUNC implies one constant set. +/// +/// 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. +/// +/// 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`, 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. +inline constexpr std::size_t calculix_ntens = 6; +inline constexpr std::size_t calculix_nstiff = + calculix_ntens * (calculix_ntens + 1) / 2; + +/// 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}; ///< 1-based + int iint{0}; ///< 1-based + int kode{0}; ///< -100 - #constants + const double* mprops{nullptr}; + 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 + 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}; +}; + +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. +/// +/// 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)"); + + // 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, + "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: 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) { + 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; 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; + + // 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) { + 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)); + } + + // 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) { + 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 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; + if (a.statev_old) + std::copy_n(a.statev_old + point_offset, nstatv, statev_point); + } + + // COLUMN-major (material_point_evaluator.h:82), which the packing below needs. + double ddsdde36[calculix_ntens * calculix_ntens] = {0.0}; + + // 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; + 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); + + // 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) + 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 symbol @p FUNC bound to model @p MODELNAME. Place in +/// exactly ONE translation unit of the library ccx loads for `@LIB,FUNC`. +/// +/// 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, \ + 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..a08390d --- /dev/null +++ b/tests/test_calculix_interface.cpp @@ -0,0 +1,909 @@ +#include + +#include +#include +#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" + +// 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") +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") +NUMSIM_MATERIALS_DEFINE_CALCULIX_BEHAVIOUR( + numsim::materials::material_policy_default, clx_throws_, "THROWSCLX") + +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 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"); + 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}; + +/// 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)); +} + +/// 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> { +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; 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> { +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; +}; + +/// 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"); + 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); + + registry::config th; + th.strain_source = "stepper"; + th.stress_source = "probe"; + registry::instance().register_model("THROWSCLX", build_throws, th); + } +}; +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. +/// 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; + static void handler(const char* msg) { + last = msg; + ++count; + } + FatalProbe() { + last.clear(); + count = 0; + u::set_fatal_handler(&handler); + } + /// 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*, + 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); + +/// Defaults for a well-formed solid-3D call; tests set only what they vary. +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 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; + 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 << ")"; +} + +/// 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; + 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) 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}; + 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) { + // 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)); + } + } + } +} + +// --------------------------------------------------------------------------- +// State variables: the full-array indexing contract +// --------------------------------------------------------------------------- + +/// 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"; + + 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)); + + // 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; + + // --- 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); + // ccx 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, 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}; + + 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 — 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"; + 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: 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) + << "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 +// --------------------------------------------------------------------------- + +/// Nothing rotates around the native hook, so ignoring iorien is wrong-frame. +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); + // 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); +} + +/// 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 count as -100 - nconst; too few constants is a setup fault. +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. +/// 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}; + 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