CalculiX target: dlopen'd .so material (no ccx recompile per material)

CMakeLists.txt

@@ -181,6 +181,7 @@ add_library(numsim_codegen STATIC
     src/targets/standalone_cxx.cpp
     src/targets/moose_material.cpp
     src/targets/numsim_material.cpp
+    src/targets/calculix_external.cpp
     src/targets/target_factory.cpp)
 add_library(numsim::codegen ALIAS numsim_codegen)
 target_link_libraries(numsim_codegen

README.md

@@ -28,6 +28,7 @@ boundary using tmech's adaptors (`full`, `voigt`, `abq_std`).
 | `StandaloneCxxTarget` | ✓ | Single inline header with the generic compute function |
 | `MooseMaterialTarget` | ✓ | `.h` + `.C` pair: Material class with `validParams`, constructor, `computeQpProperties`, optional `Jacobian_mult` consistent tangent |
 | `NumSimMaterialTarget` | ✓ | numsim-materials material header + JSON config. Two sub-contracts: rate path (exactly one scalar state variable + one evolution equation — the rk_integrator contract); return-map path (one scalar Newton unknown, plus additional scalar or tensor history state via update equations) |
+| `CalculiXExternalTarget` (`calculix`) | ✓ elastic | `<Model>_ext.cpp` → `lib<MODEL>.so`: CalculiX external behaviour loaded at RUNTIME via `dlopen` (deck `NAME=@<MODEL>_NCG_UMAT`) — build ccx **once**, then **no recompile per material**. tmech `abq_std` Voigt boundary (native `emec`/`stre`/`stiff(21)`); per-thread evaluator. CI: ABI gate vs. an independent oracle (`tests/generated/calculix_check_driver.cpp`). Real-ccx golden-file validation lives in the `tests/calculix/` harness (issue #128) |
 | `AbaqusUMATTarget`    | planned | Fortran-callable `extern "C"` UMAT with Voigt boundary |
 | `AnsysUSERMATTarget`  | planned | Fortran-callable USERMAT |
 | `LSDynaUMATTarget`    | planned | LS-DYNA convention |

include/numsim_codegen/targets/calculix_external.h

@@ -0,0 +1,39 @@
+#ifndef NUMSIM_CODEGEN_TARGETS_CALCULIX_EXTERNAL_H
+#define NUMSIM_CODEGEN_TARGETS_CALCULIX_EXTERNAL_H
+
+#include <numsim_codegen/code_emit/linear_algebra_emitter.h>
+#include <numsim_codegen/targets/target.h>
+
+namespace numsim::codegen {
+
+// CalculiX external-behaviour target: emits a `.cpp` compiling to `lib<MODEL>.so`,
+// which ccx dlopens at RUNTIME — build ccx once with
+// `-DCALCULIX_EXTERNAL_BEHAVIOURS_SUPPORT -ldl`, then select per material with
+// `*MATERIAL, NAME=@<MODEL>_NCG_UMAT` (ccx uppercases, splits `@<LIB>_<FUNC>`,
+// dlopens `lib<LIB>.so`, dlsyms `<FUNC>`). No recompile per material.
+//
+// The hook passes NATIVE quantities under Abaqus-flavoured names — STRAN1=emec
+// (tensorial), STRESS=stre, DDSDDE=stiff(21) — so the boundary is a tmech
+// `abq_std` adaptor plus the stiff(21) packing in calculix_boundary.h. The
+// evaluator is thread_local: ccx runs the element loop multi-threaded.
+//
+// SCOPE: stateless recipes only — one rank-2 strain in, one rank-2 stress out,
+// one rank-4 tangent, scalar parameters (the *USER MATERIAL constants). The
+// STATEV round-trip is the numsim-materials-backed follow-up (see #160).
+class CalculiXExternalTarget : public Target {
+public:
+  explicit CalculiXExternalTarget(
+      LinearAlgebraEmitter const &la = default_linear_algebra_emitter())
+      : m_la(la) {}
+  CalculiXExternalTarget(LinearAlgebraEmitter const &&) = delete;
+  [[nodiscard]] auto emit(ConstitutiveModel const &model) const
+      -> std::vector<EmittedFile> override;
+  [[nodiscard]] auto target_name() const -> std::string override;
+
+private:
+  LinearAlgebraEmitter const &m_la;
+};
+
+} // namespace numsim::codegen
+
+#endif // NUMSIM_CODEGEN_TARGETS_CALCULIX_EXTERNAL_H

src/targets/calculix_boundary.h

@@ -0,0 +1,147 @@
+#ifndef NUMSIM_CODEGEN_SRC_TARGETS_CALCULIX_BOUNDARY_H
+#define NUMSIM_CODEGEN_SRC_TARGETS_CALCULIX_BOUNDARY_H
+
+// Shared by both CalculiX targets (linked-in umat_user_ and the dlopen'd
+// plugin): the stateless-elastic scope rules and the `stiff(21)` packing are one
+// contract, so they live here rather than being copy-pasted. The lifecycle stays
+// per-target.
+
+#include <numsim_codegen/recipe.h>
+
+#include <cstddef>
+#include <optional>
+#include <ostream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace numsim::codegen::detail {
+
+struct CalculiXTensorArg {
+  std::string name;
+  std::size_t dim = 0;
+  std::size_t rank = 0;
+};
+
+// The collected, validated boundary variables of a stateless-elastic recipe.
+struct CalculiXScope {
+  std::vector<std::string> params; // → *USER MATERIAL constants, in order
+  CalculiXTensorArg strain;        // the one rank-2 (symmetric) tensor input
+  CalculiXTensorArg stress;        // the one rank-2 (symmetric) tensor output
+  CalculiXTensorArg tangent;       // the one rank-4 consistent tangent
+};
+
+// Scan the canonical argument list (post-emit order, issue #77) and enforce the
+// stateless-elastic scope. `label` prefixes diagnostics. Throws on violation.
+[[nodiscard]] inline auto
+scan_calculix_scope(ConstitutiveModel const &model, char const *label)
+    -> CalculiXScope {
+  CalculiXScope scope;
+  std::optional<CalculiXTensorArg> strain, stress, tangent;
+
+  auto const reject = [&](std::string const &what) {
+    throw std::runtime_error(
+        std::string(label) + ": recipe '" + model.name() + "' " + what +
+        ". This target's first cut supports stateless materials only (one "
+        "symmetric strain input, one stress output, one consistent tangent, "
+        "plus scalar parameters → the *USER MATERIAL constants); state "
+        "variables, scalar inputs and rate/implicit forms are a follow-up.");
+  };
+
+  // abq_std is a symmetric 6-component adaptor, so a non-symmetric leaf (e.g.
+  // roles::DeformationGradient) would be silently truncated.
+  auto input_is_symmetric = [&](std::string const &name) {
+    for (auto const &s : model.inputs())
+      if (s.name == name) return s.role.is_symmetric;
+    return false;
+  };
+  auto output_is_symmetric = [&](std::string const &name) {
+    for (auto const &o : model.outputs())
+      if (o.name == name) return o.role.is_symmetric;
+    return false;
+  };
+
+  for (auto const &a : canonical_arguments(RecipeView{model})) {
+    switch (a.role) {
+    case ArgSpec::Role::ScalarParam:
+      scope.params.push_back(a.name);
+      break;
+    case ArgSpec::Role::TensorInput:
+      if (strain) reject("has more than one tensor input");
+      if (a.dim != 3 || a.rank != 2)
+        reject("has a strain input '" + a.name +
+               "' that is not a 3D rank-2 tensor");
+      if (!input_is_symmetric(a.name))
+        reject("has a non-symmetric tensor input '" + a.name +
+               "' (e.g. a deformation gradient); CalculiX's abq_std Voigt "
+               "boundary is symmetric — use a symmetric strain measure "
+               "(roles::Strain)");
+      strain = CalculiXTensorArg{a.name, a.dim, a.rank};
+      break;
+    case ArgSpec::Role::TensorOutput:
+      if (stress) reject("has more than one tensor output");
+      if (a.dim != 3 || a.rank != 2)
+        reject("has a stress output '" + a.name +
+               "' that is not a 3D rank-2 tensor");
+      if (!output_is_symmetric(a.name))
+        reject("has a non-symmetric tensor output '" + a.name +
+               "' (CalculiX stress storage is symmetric)");
+      stress = CalculiXTensorArg{a.name, a.dim, a.rank};
+      break;
+    case ArgSpec::Role::TensorTangentOutput:
+      if (tangent) reject("has more than one consistent tangent");
+      if (a.dim != 3) reject("has a consistent tangent that is not 3D");
+      tangent = CalculiXTensorArg{a.name, a.dim, a.rank};
+      break;
+    case ArgSpec::Role::ScalarInput:
+      reject("has a scalar input '" + a.name + "'");
+      break;
+    case ArgSpec::Role::TimeStep:
+      reject("uses the time step (rate/implicit form)");
+      break;
+    case ArgSpec::Role::ScalarOutput:
+      reject("has a scalar output '" + a.name + "'");
+      break;
+    case ArgSpec::Role::StateOld:
+    case ArgSpec::Role::StateCurrentRead:
+    case ArgSpec::Role::NewtonStateOut:
+      reject("has an internal state variable '" + a.name + "'");
+      break;
+    }
+  }
+
+  if (!strain) reject("has no tensor (strain) input");
+  if (!stress) reject("has no tensor (stress) output");
+  if (!tangent)
+    reject("has no consistent tangent (add_algorithmic_tangent is required so "
+           "CalculiX gets the material stiffness)");
+
+  scope.strain = *strain;
+  scope.stress = *stress;
+  scope.tangent = *tangent;
+  return scope;
+}
+
+// Pack a 6x6 (`d6_name`) into stiff(21): column-major upper triangle,
+// k = i + j*(j+1)/2 (0-based i<=j), symmetrized as umat_abaqus.f:335-355 does.
+// An asymmetric tangent (non-associative plasticity) loses its antisymmetric
+// part — stiff(21) has no room for it. ccx passes icmd==3 for stress only.
+inline void emit_stiff21_packing(std::ostream &os, std::string const &d6_name,
+                                 std::string const &stiff_ptr,
+                                 std::string const &icmd_expr,
+                                 std::string const &indent) {
+  os << indent << "// Pack the 6x6 into CalculiX's symmetric stiff(21):\n";
+  os << indent << "// column-major upper, k = i + j*(j+1)/2 (0-based i<=j),\n";
+  os << indent << "// major-symmetrized. icmd==3 → CalculiX wants stress only.\n";
+  os << indent << "if (" << icmd_expr << " != 3) {\n";
+  os << indent << "  for (int j = 0; j < 6; ++j)\n";
+  os << indent << "    for (int i = 0; i <= j; ++i)\n";
+  os << indent << "      " << stiff_ptr << "[i + j * (j + 1) / 2] =\n";
+  os << indent << "          0.5 * (" << d6_name << "[i * 6 + j] + " << d6_name
+     << "[j * 6 + i]);\n";
+  os << indent << "}\n";
+}
+
+} // namespace numsim::codegen::detail
+
+#endif // NUMSIM_CODEGEN_SRC_TARGETS_CALCULIX_BOUNDARY_H

src/targets/calculix_external.cpp

@@ -0,0 +1,161 @@
+#include <numsim_codegen/targets/calculix_external.h>
+
+#include "calculix_boundary.h"
+
+#include <numsim_codegen/recipe.h>
+
+#include <cctype>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace numsim::codegen {
+
+namespace {
+
+// The `calculixptr` typedef from call_external_umat_user.c. Quantities are
+// NATIVE: STRAN1=emec (tensorial, {11,22,33,12,13,23}), STRESS=stre,
+// DDSDDE=stiff(21), MPROPS=the *USER MATERIAL constants (NPROPS is ccx's `kode`,
+// not a count). Trailing `int size` is the amat length, BY VALUE. Only the
+// arguments we use are named.
+constexpr char const *kExternalSignature = R"(extern "C" void NCG_UMAT(
+    char const *   /*amat*/,   int const *    /*iel*/,    int const *    /*iint*/,
+    int const *    /*NPROPS*/, double const * MPROPS,     double const * STRAN1,
+    double const * /*STRAN0*/, double const * /*beta*/,   double const * /*F0*/,
+    double const * /*voj*/,    double const * /*F1*/,     double const * /*vj*/,
+    int const *    /*ithermal*/, double const * /*TEMP1*/, double const * /*DTIME*/,
+    double const * /*time*/,   double const * /*ttime*/,  int const *    icmd,
+    int const *    /*ielas*/,  int const *    /*mi*/,     int const *    /*NSTATV*/,
+    [[maybe_unused]] double const * STATEV0,
+    [[maybe_unused]] double *       STATEV1,
+    double *       STRESS,
+    double *       DDSDDE,      int const *    /*iorien*/, double const * /*pgauss*/,
+    double const * /*orab*/,   double *       /*PNEWDT*/, int const *    /*ipkon*/,
+    int            /*size*/))";
+
+// The library name ccx parses out of `@<LIB>_NCG_UMAT`: uppercase, alphanumeric
+// only. Underscores are dropped so the first `_` splits LIB from FUNC.
+// LOSSY — Mat_1, MAT1 and mat1 all map to MAT1, and the .so filename is the only
+// namespace (every plugin exports NCG_UMAT), so callers must avoid collisions.
+auto library_name(std::string const &model) -> std::string {
+  std::string out;
+  for (char c : model)
+    if (std::isalnum(static_cast<unsigned char>(c)))
+      out.push_back(
+          static_cast<char>(std::toupper(static_cast<unsigned char>(c))));
+  return out;
+}
+
+} // namespace
+
+auto CalculiXExternalTarget::emit(ConstitutiveModel const &model) const
+    -> std::vector<EmittedFile> {
+  std::string const body = model.emit_compute_function(m_la);
+  bool const needs_la = body.find(m_la.usage_marker()) != std::string::npos;
+
+  auto const scope =
+      detail::scan_calculix_scope(model, "CalculiXExternalTarget");
+
+  auto const lib = library_name(model.name());
+  if (lib.empty()) {
+    throw std::runtime_error(
+        "CalculiXExternalTarget: recipe '" + model.name() +
+        "' has no alphanumeric characters in its name, so it maps to an empty "
+        "library name (ccx would reject the deck). Rename the model.");
+  }
+
+  std::ostringstream os;
+  os << "// Auto-generated by numsim-codegen. Do not edit.\n";
+  os << "//\n";
+  os << "// CalculiX external behaviour, dlopen'd at runtime. Build as:\n";
+  os << "//   g++ -std=c++23 -O2 -fPIC -shared this.cpp -I<tmech/include> -o lib"
+     << lib << ".so\n";
+  os << "// Deck:  *MATERIAL, NAME=@" << lib << "_NCG_UMAT\n";
+  os << "//        *USER MATERIAL, CONSTANTS=" << scope.params.size() << "\n";
+  os << "// CONSTANTS must equal the parameter count (" << scope.params.size()
+     << "): MPROPS is read positionally, unchecked.\n";
+  os << "//\n";
+  os << "// Each call reads constants + strain, evaluates, packs stress +\n";
+  os << "// tangent. Constants are read FRESH per call: ccx interpolates them by\n";
+  os << "// temperature, so they can vary.\n\n";
+  os << "#include <tmech/tmech.h>\n";
+  os << "#include <cmath>\n";
+  if (needs_la) {
+    for (auto const &inc : m_la.includes()) os << "#include " << inc << "\n";
+  }
+  os << "\n";
+  os << body;
+  os << "\n";
+
+  // A per-thread stateless evaluator: no cached constants.
+  os << "namespace {\n\n";
+  os << "// Stateless, cached per thread (ccx threads the element loop). Holds no\n";
+  os << "// constants: they are read from MPROPS each call, since ccx may vary them\n";
+  os << "// and one library may serve several *MATERIAL blocks.\n";
+  os << "struct ncg_material {\n";
+  os << "  void evaluate([[maybe_unused]] double const *mprops,\n";
+  os << "                double const *strain_in,\n";
+  os << "                double *stress_out, double *tangent_stiff, int icmd) const {\n";
+  for (std::size_t k = 0; k < scope.params.size(); ++k) {
+    os << "    double const " << scope.params[k] << " = mprops[" << k << "];\n";
+  }
+  os << "    tmech::adaptor<double const, 3, 2, tmech::abq_std<3, false>> "
+     << scope.strain.name << "_ad(strain_in);\n";
+  os << "    tmech::adaptor<double, 3, 2, tmech::abq_std<3, false>> "
+     << scope.stress.name << "_ad(stress_out);\n";
+  os << "    double " << scope.tangent.name << "_D6[36];\n";
+  os << "    tmech::adaptor<double, 3, 4, tmech::abq_std<3, false>> "
+     << scope.tangent.name << "_ad(" << scope.tangent.name << "_D6);\n";
+  os << "\n    " << model.name() << "_compute(\n";
+  bool first = true;
+  for (auto const &a : canonical_arguments(RecipeView{model})) {
+    if (!first) os << ",\n";
+    first = false;
+    switch (a.role) {
+    case ArgSpec::Role::ScalarParam:
+      os << "        " << a.name;
+      break;
+    case ArgSpec::Role::TensorInput:
+    case ArgSpec::Role::TensorOutput:
+    case ArgSpec::Role::TensorTangentOutput:
+      os << "        " << a.name << "_ad";
+      break;
+    case ArgSpec::Role::ScalarInput:
+    case ArgSpec::Role::StateOld:
+    case ArgSpec::Role::StateCurrentRead:
+    case ArgSpec::Role::TimeStep:
+    case ArgSpec::Role::ScalarOutput:
+    case ArgSpec::Role::NewtonStateOut:
+      throw std::runtime_error(
+          "CalculiXExternalTarget: internal error — unrejected role reached the "
+          "call site for '" + a.name + "'");
+    }
+  }
+  os << ");\n";
+  detail::emit_stiff21_packing(os, scope.tangent.name + "_D6", "tangent_stiff",
+                               "icmd", "    ");
+  os << "  }\n";
+  os << "};\n\n";
+
+  os << "ncg_material const &thread_state() {\n";
+  os << "  thread_local ncg_material const evaluator;\n";
+  os << "  return evaluator;\n";
+  os << "}\n\n";
+  os << "} // namespace\n\n";
+
+  // The exported entry point. STATEV0/STATEV1 are [[maybe_unused]] but stay
+  // NAMED for the future history round-trip.
+  os << kExternalSignature << " {\n";
+  os << "  thread_state().evaluate(MPROPS, STRAN1, STRESS, DDSDDE, *icmd);\n";
+  os << "}\n";
+
+  return {EmittedFile{model.name() + "_ext.cpp", os.str(), "",
+                      EmittedFile::Kind::Source}};
+}
+
+auto CalculiXExternalTarget::target_name() const -> std::string {
+  return "CalculiXExternal";
+}
+
+} // namespace numsim::codegen