diff --git a/CMakeLists.txt b/CMakeLists.txt index 3db9eb9..d157f7c 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/README.md b/README.md index da4e73e..f1e5eb2 100644 --- a/README.md +++ b/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 | `_ext.cpp` → `lib.so`: CalculiX external behaviour loaded at RUNTIME via `dlopen` (deck `NAME=@_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 | diff --git a/include/numsim_codegen/targets/calculix_external.h b/include/numsim_codegen/targets/calculix_external.h new file mode 100644 index 0000000..56cf928 --- /dev/null +++ b/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 +#include + +namespace numsim::codegen { + +// CalculiX external-behaviour target: emits a `.cpp` compiling to `lib.so`, +// which ccx dlopens at RUNTIME — build ccx once with +// `-DCALCULIX_EXTERNAL_BEHAVIOURS_SUPPORT -ldl`, then select per material with +// `*MATERIAL, NAME=@_NCG_UMAT` (ccx uppercases, splits `@_`, +// dlopens `lib.so`, dlsyms ``). 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 override; + [[nodiscard]] auto target_name() const -> std::string override; + +private: + LinearAlgebraEmitter const &m_la; +}; + +} // namespace numsim::codegen + +#endif // NUMSIM_CODEGEN_TARGETS_CALCULIX_EXTERNAL_H diff --git a/src/targets/calculix_boundary.h b/src/targets/calculix_boundary.h new file mode 100644 index 0000000..c3d649c --- /dev/null +++ b/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 + +#include +#include +#include +#include +#include +#include + +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 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 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 diff --git a/src/targets/calculix_external.cpp b/src/targets/calculix_external.cpp new file mode 100644 index 0000000..829cc04 --- /dev/null +++ b/src/targets/calculix_external.cpp @@ -0,0 +1,161 @@ +#include + +#include "calculix_boundary.h" + +#include + +#include +#include +#include +#include +#include + +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 `@_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(c))) + out.push_back( + static_cast(std::toupper(static_cast(c)))); + return out; +} + +} // namespace + +auto CalculiXExternalTarget::emit(ConstitutiveModel const &model) const + -> std::vector { + 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 -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 \n"; + os << "#include \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> " + << scope.strain.name << "_ad(strain_in);\n"; + os << " tmech::adaptor> " + << scope.stress.name << "_ad(stress_out);\n"; + os << " double " << scope.tangent.name << "_D6[36];\n"; + os << " tmech::adaptor> " + << 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 diff --git a/src/targets/target_factory.cpp b/src/targets/target_factory.cpp index 1c288ff..83772f5 100644 --- a/src/targets/target_factory.cpp +++ b/src/targets/target_factory.cpp @@ -1,35 +1,64 @@ #include +#include #include #include #include +#include #include #include namespace numsim::codegen { +namespace { + +// Single source of truth for the target registry — both `target_names()` and +// `make_target()` iterate this, so a target can't appear in one list but not +// the other (review: arch #5). Default first (matches `default_target_name`). +// The LA-backed targets default-construct their `LinearAlgebraEmitter const&` +// from the static `default_linear_algebra_emitter()` accessor (safe lifetime, +// not the `=delete`'d rvalue ctor). +struct Entry { + std::string_view name; + std::unique_ptr (*make)(); +}; + +auto registry() -> std::vector const & { + static std::vector const entries{ + {"numsim_material", + [] { return std::unique_ptr(std::make_unique()); }}, + {"standalone", + [] { return std::unique_ptr(std::make_unique()); }}, + {"moose", + [] { return std::unique_ptr(std::make_unique()); }}, + {"calculix", + [] { return std::unique_ptr(std::make_unique()); }}, + }; + return entries; +} + +} // namespace + auto target_names() -> std::vector const & { - // Default first (kept in sync with make_target + default_target_name). - static std::vector const names{"numsim_material", - "standalone", "moose"}; + static std::vector const names = [] { + std::vector n; + for (auto const &e : registry()) n.push_back(e.name); + return n; + }(); return names; } auto make_target(std::string_view name) -> std::unique_ptr { - // The LA-backed targets default-construct their `LinearAlgebraEmitter const&` - // from the static `default_linear_algebra_emitter()` accessor — safe lifetime, - // and not the `=delete`'d rvalue ctor. - if (name == "numsim_material") return std::make_unique(); - if (name == "standalone") return std::make_unique(); - if (name == "moose") return std::make_unique(); + for (auto const &e : registry()) + if (e.name == name) return e.make(); std::string msg = "make_target: unknown target '"; msg.append(name); msg += "'. Known:"; - for (auto const known : target_names()) { + for (auto const &e : registry()) { msg += ' '; - msg.append(known); + msg.append(e.name); } throw std::runtime_error(msg); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f0aa040..322cd63 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,6 +35,7 @@ add_executable(numsim_codegen_tests LinearAlgebraEmitterTest.cpp NumSimMaterialTargetTest.cpp TargetFactoryTest.cpp + CalculiXTargetTest.cpp ) target_link_libraries(numsim_codegen_tests @@ -389,6 +390,46 @@ else() "NUMSIM_CODEGEN_FETCH_MATERIALS=OFF") endif() +# ─── CalculiX external-behaviour end-to-end gate ───────────────────────────── +# +# Emit full isotropic linear elasticity as (a) a StandaloneCxx header (for the +# FD tangent + stress oracle) and (b) a CalculiXExternal `NCG_UMAT` .cpp, then +# compile a driver that calls the emitted `NCG_UMAT` exactly as `ccx` would +# (via its external/dlopen ABI) and checks stre(6)/stiff(21) against an +# INDEPENDENT isotropic oracle. Proves the abq_std Voigt boundary + column-major +# stiff packing with zero external deps — the correctness lock behind the real +# `ccx` run (examples/calculix/). +add_executable(generate_calculix_check + generated/generate_calculix_check.cpp) +target_link_libraries(generate_calculix_check + PRIVATE numsim::codegen numsim_codegen_warnings) + +set(GENERATED_ELASTIC_HEADER + ${CMAKE_CURRENT_BINARY_DIR}/generated/LinearElastic.h) +set(GENERATED_ELASTIC_EXT + ${CMAKE_CURRENT_BINARY_DIR}/generated/LinearElastic_ext.cpp) +add_custom_command( + OUTPUT ${GENERATED_ELASTIC_HEADER} ${GENERATED_ELASTIC_EXT} + COMMAND ${CMAKE_COMMAND} -E make_directory + ${CMAKE_CURRENT_BINARY_DIR}/generated + COMMAND $ + ${GENERATED_ELASTIC_HEADER} ${GENERATED_ELASTIC_EXT} + DEPENDS generate_calculix_check + COMMENT "Generating CalculiX external material + standalone header" + VERBATIM) + +add_executable(calculix_check_driver + generated/calculix_check_driver.cpp + ${GENERATED_ELASTIC_HEADER} + ${GENERATED_ELASTIC_EXT}) +target_include_directories(calculix_check_driver + PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/generated + ${CMAKE_CURRENT_SOURCE_DIR}/support) +target_link_libraries(calculix_check_driver + PRIVATE numsim::codegen numsim_codegen_warnings_light + GTest::gtest GTest::gtest_main) +gtest_discover_tests(calculix_check_driver) + # ─── Header-purity guard ─────────────────────────────────────────────────── # # Public headers under include/ must never #include heavy template libraries diff --git a/tests/CalculiXTargetTest.cpp b/tests/CalculiXTargetTest.cpp new file mode 100644 index 0000000..27a915e --- /dev/null +++ b/tests/CalculiXTargetTest.cpp @@ -0,0 +1,188 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +namespace numsim::codegen { + +namespace { + +// Isotropic linear elasticity with its tangent: the supported stateless shape. +auto build_full_elastic(std::string name = "LinearElastic") -> ConstitutiveModel { + using namespace numsim::cas; + ConstitutiveModel m(std::move(name)); + auto lambda = m.add_parameter("lambda", 1.0); + auto mu = m.add_parameter("mu", 0.5); + auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain); + auto I = make_expression(std::size_t{3}, std::size_t{2}); + auto tr = make_expression(eps); + auto sigma = + lambda * make_expression(I, tr) + + 2 * mu * eps; + m.add_output("stress", sigma, roles::Stress); + m.add_algorithmic_tangent("dstress_deps", "stress", "eps"); + return m; +} + +// A non-symmetric F input: abq_std would silently truncate its skew part. +auto build_deformation_gradient_recipe() -> ConstitutiveModel { + using namespace numsim::cas; + ConstitutiveModel m("DefGrad"); + auto mu = m.add_parameter("mu", 0.5); + auto F = m.add_tensor_input("F", 3, 2, roles::DeformationGradient); + m.add_output("stress", 2 * mu * F, roles::Stress); + m.add_algorithmic_tangent("dstress_dF", "stress", "F"); + return m; +} + +} // namespace + +// ── Emission ───────────────────────────────────────────────────────────────── + +TEST(CalculiXTarget, EmitsSingleSourceFileNamedExt) { + CalculiXExternalTarget target; + auto files = target.emit(build_full_elastic()); + ASSERT_EQ(files.size(), 1u); + EXPECT_EQ(files[0].filename, "LinearElastic_ext.cpp"); + EXPECT_EQ(files[0].kind, EmittedFile::Kind::Source); + EXPECT_TRUE(files[0].install_subdir.empty()); +} + +TEST(CalculiXTarget, ExposesNcgUmatAndAbqStdBoundary) { + CalculiXExternalTarget target; + auto const files = target.emit(build_full_elastic()); + auto const &src = files[0].contents; + EXPECT_NE(src.find("extern \"C\" void NCG_UMAT("), std::string::npos) << src; + EXPECT_NE(src.find("LinearElastic_compute("), std::string::npos); + // abq_std ordering {11,22,33,12,13,23}; no plain voigt, no ×2 engineering shear. + EXPECT_NE(src.find("tmech::abq_std<3, false>"), std::string::npos); + EXPECT_EQ(src.find("tmech::voigt<"), std::string::npos); +} + +TEST(CalculiXTarget, PacksStiffColumnMajorUpper) { + CalculiXExternalTarget target; + auto const files = target.emit(build_full_elastic()); + auto const &src = files[0].contents; + EXPECT_NE(src.find("stiff[i + j * (j + 1) / 2]"), std::string::npos); + EXPECT_NE(src.find("!= 3"), std::string::npos); // icmd stress-only guard +} + +// Constants read from MPROPS per call, evaluator stateless: the cache-bug gate. +TEST(CalculiXTarget, ReadsConstantsPerCallNotCached) { + CalculiXExternalTarget target; + auto const files = target.emit(build_full_elastic()); + auto const &src = files[0].contents; + EXPECT_NE(src.find("double const lambda = mprops[0];"), std::string::npos); + EXPECT_NE(src.find("thread_local ncg_material const"), std::string::npos); + EXPECT_EQ(src.find(".emplace(mprops)"), std::string::npos) + << "constants must not be cached from the first call"; +} + +// ccx uppercases the deck name and drops '_': "J2_Plastic" -> libJ2PLASTIC.so. +TEST(CalculiXTarget, DeckNameIsUppercasedAlnum) { + CalculiXExternalTarget target; + auto const files = target.emit(build_full_elastic("J2_Plastic")); + auto const &src = files[0].contents; + EXPECT_NE(src.find("libJ2PLASTIC.so"), std::string::npos) << src; + EXPECT_NE(src.find("@J2PLASTIC_NCG_UMAT"), std::string::npos) << src; +} + +TEST(CalculiXTarget, FactorySelector) { + auto target = make_target("calculix"); + ASSERT_NE(target, nullptr); + EXPECT_EQ(target->target_name(), "CalculiXExternal"); +} + +// ── Scope guards (stateless linear-elastic first cut) ──────────────────────── + +TEST(CalculiXTarget, RejectsMissingConsistentTangent) { + using namespace numsim::cas; + CalculiXExternalTarget target; + ConstitutiveModel m("NoTangent"); + auto mu = m.add_parameter("mu", 0.5); + auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain); + m.add_output("stress", 2 * mu * eps, roles::Stress); // no add_algorithmic_tangent + EXPECT_THROW([[maybe_unused]] auto const discarded = target.emit(m), + std::runtime_error); +} + +TEST(CalculiXTarget, RejectsScalarInput) { + using namespace numsim::cas; + CalculiXExternalTarget target; + ConstitutiveModel m("WithScalarInput"); + auto mu = m.add_parameter("mu", 0.5); + auto T = m.add_scalar_input("T", roles::Temperature); + auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain); + m.add_output("stress", 2 * mu * (1 + T) * eps, roles::Stress); + m.add_algorithmic_tangent("dstress_deps", "stress", "eps"); + EXPECT_THROW([[maybe_unused]] auto const discarded = target.emit(m), + std::runtime_error); +} + +TEST(CalculiXTarget, RejectsMultipleTensorInputs) { + using namespace numsim::cas; + CalculiXExternalTarget target; + ConstitutiveModel m("TwoTensorInputs"); + auto mu = m.add_parameter("mu", 0.5); + auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain); + auto eps_p = m.add_tensor_input( + "eps_p", 3, 2, + Role{.name = "plastic_strain", .is_symmetric = true, .expected_rank = 2}); + m.add_output("stress", 2 * mu * (eps - eps_p), roles::Stress); + m.add_algorithmic_tangent("dstress_deps", "stress", "eps"); + EXPECT_THROW([[maybe_unused]] auto const discarded = target.emit(m), + std::runtime_error); +} + +TEST(CalculiXTarget, RejectsMultipleTensorOutputs) { + using namespace numsim::cas; + CalculiXExternalTarget target; + ConstitutiveModel m("TwoOutputs"); + auto mu = m.add_parameter("mu", 0.5); + auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain); + m.add_output("stress", 2 * mu * eps, roles::Stress); + m.add_output("extra", 3 * mu * eps, roles::Stress); // second rank-2 output + m.add_algorithmic_tangent("dstress_deps", "stress", "eps"); + EXPECT_THROW([[maybe_unused]] auto const discarded = target.emit(m), + std::runtime_error); +} + +// H3: a non-symmetric input must be rejected, not silently stripped of skew. +TEST(CalculiXTarget, RejectsNonSymmetricTensorInput) { + CalculiXExternalTarget target; + EXPECT_THROW([[maybe_unused]] auto const discarded = + target.emit(build_deformation_gradient_recipe()), + std::runtime_error); +} + +TEST(CalculiXTarget, RejectsStateVariable) { + using namespace numsim::cas; + CalculiXExternalTarget target; + ConstitutiveModel m("WithState"); + auto K = m.add_parameter("K", 1.0); + auto alpha = + m.add_scalar_state_variable("alpha", make_expression(0.0)); + m.add_scalar_evolution_equation(alpha, K * alpha.current); + m.enable_local_newton(); + EXPECT_THROW([[maybe_unused]] auto const discarded = target.emit(m), + std::runtime_error); +} + +// A name with no alphanumerics gives an empty library name: throw at emit, not +// inside ccx. +TEST(CalculiXTarget, RejectsEmptyLibraryName) { + CalculiXExternalTarget target; + EXPECT_THROW([[maybe_unused]] auto const discarded = + target.emit(build_full_elastic("_")), + std::runtime_error); +} + +} // namespace numsim::codegen diff --git a/tests/TargetFactoryTest.cpp b/tests/TargetFactoryTest.cpp index 8b247ae..cfe724c 100644 --- a/tests/TargetFactoryTest.cpp +++ b/tests/TargetFactoryTest.cpp @@ -16,6 +16,7 @@ TEST(TargetFactory, ConstructsEachKnownTarget) { EXPECT_EQ(make_target("numsim_material")->target_name(), "NumSimMaterial"); EXPECT_EQ(make_target("standalone")->target_name(), "StandaloneCxx"); EXPECT_EQ(make_target("moose")->target_name(), "MooseMaterial"); + EXPECT_EQ(make_target("calculix")->target_name(), "CalculiXExternal"); } TEST(TargetFactory, DefaultIsNumSimMaterial) { @@ -27,7 +28,7 @@ TEST(TargetFactory, NamesListDefaultFirstAndAreAllConstructible) { auto const &names = target_names(); // Pinned count: adding a target without updating the list (or vice-versa) // trips this — the cheap guard for the "three places, no enforcement" shape. - EXPECT_EQ(names.size(), 3u); + EXPECT_EQ(names.size(), 4u); ASSERT_FALSE(names.empty()); EXPECT_EQ(names.front(), default_target_name); // Every advertised name must actually construct. diff --git a/tests/generated/calculix_check_driver.cpp b/tests/generated/calculix_check_driver.cpp new file mode 100644 index 0000000..69f1664 --- /dev/null +++ b/tests/generated/calculix_check_driver.cpp @@ -0,0 +1,234 @@ +// CalculiX end-to-end gate driver. Two layers, no external dependencies: +// Recipe: FD-verify the emitted tangent through LinearElastic_compute and +// anchor the stress to the closed-form isotropic law. +// ABI: call the emitted NCG_UMAT as ccx would, and check stre(6)/stiff(21) +// against an INDEPENDENT oracle — pinning the Voigt boundary and the +// column-major packing before ccx is built. A negative control confirms the +// oracle discriminates order, not just values. + +#include "LinearElastic.h" // StandaloneCxx: LinearElastic_compute(...) + +#include "numerical_tangent_verifier.h" + +#include + +#include + +#include +#include +#include + +// ── The emitted external-behaviour ABI (LinearElastic_ext.cpp) — the exact +// `calculixptr` signature ccx dlsym's and calls. ───────────────────────────── +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, + double const *STATEV0, double *STATEV1, double *STRESS, double *DDSDDE, + int const *iorien, double const *pgauss, double const *orab, double *PNEWDT, + int const *ipkon, int size); + +namespace { + +constexpr double kLambda = 1.3; +constexpr double kMu = 0.7; + +using T2 = tmech::tensor; +using T4 = tmech::tensor; + +// Symmetric strain tensor from emec(6), abq_std order, tensorial shear. +auto strain_tensor_from_emec(std::array const &e) -> T2 { + T2 eps; + eps(0, 0) = e[0]; + eps(1, 1) = e[1]; + eps(2, 2) = e[2]; + eps(0, 1) = eps(1, 0) = e[3]; + eps(0, 2) = eps(2, 0) = e[4]; + eps(1, 2) = eps(2, 1) = e[5]; + return eps; +} + +// Closed-form isotropic stress σ = λ·tr(ε)·I + 2μ·ε. +auto isotropic_stress(T2 const &eps, double lambda, double mu) -> T2 { + double const tr = eps(0, 0) + eps(1, 1) + eps(2, 2); + T2 I = tmech::eye(); + return T2(lambda * tr * I + 2.0 * mu * eps); +} +auto isotropic_stress(T2 const &eps) -> T2 { + return isotropic_stress(eps, kLambda, kMu); +} + +// Independent isotropic tangent packed into stiff(21): column-major upper, +// k = i + j*(j+1)/2. Built WITHOUT touching the emitted code. +auto expected_stiff_column_major() -> std::array { + // Engineering 6×6 isotropic D matrix (abq_std order {11,22,33,12,13,23}). + double const a = kLambda + 2.0 * kMu; // diagonal normal + double const b = kLambda; // off-diagonal normal + double D[6][6] = {{a, b, b, 0, 0, 0}, {b, a, b, 0, 0, 0}, + {b, b, a, 0, 0, 0}, {0, 0, 0, kMu, 0, 0}, + {0, 0, 0, 0, kMu, 0}, {0, 0, 0, 0, 0, kMu}}; + std::array s{}; + for (int j = 0; j < 6; ++j) + for (int i = 0; i <= j; ++i) s[static_cast(i + j * (j + 1) / 2)] = D[i][j]; + return s; +} + +// Call NCG_UMAT for one integration point as ccx's dlopen path does. +void call_ext(std::array const &emec_in, std::array &stre, + std::array &stiff, double lambda = kLambda, + double mu = kMu, int icmd_val = 1) { + std::array mprops{lambda, mu}; + std::array stran1 = emec_in, stran0{}, beta{}; + int iel = 1, iint = 1, nprops = -102, icmd = icmd_val, ielas = 0, mi = 1, + nstatv = 0, iorien = 0, ipkon = 0, ithermal = 0; + double f0[9] = {}, voj = 1.0, f1[9] = {}, vj = 1.0, temp1 = 0.0, dtime = 1.0, + time_[2] = {}, ttime = 0.0, statev0 = 0.0, statev1 = 0.0, pgauss[3] = {}, + orab[7] = {}, pnewdt = 1.0; + char amat[81] = "@LINEARELASTIC_NCG_UMAT"; + stre.fill(0.0); + stiff.fill(0.0); + NCG_UMAT(amat, &iel, &iint, &nprops, mprops.data(), stran1.data(), + stran0.data(), beta.data(), f0, &voj, f1, &vj, &ithermal, &temp1, + &dtime, time_, &ttime, &icmd, &ielas, &mi, &nstatv, &statev0, + &statev1, stre.data(), stiff.data(), &iorien, pgauss, orab, &pnewdt, + &ipkon, 80); +} + +// A spread of strain states: uniaxial, pure shear, and a general symmetric one. +auto sample_strains() -> std::vector> { + return { + {{0.01, 0.0, 0.0, 0.0, 0.0, 0.0}}, // uniaxial ε11 + {{0.0, 0.0, 0.0, 0.02, 0.0, 0.0}}, // pure shear ε12 + {{0.005, -0.003, 0.002, 0.004, -0.001, 0.006}}, // general + }; +} + +} // namespace + +// ── Phase 0: recipe correctness through the standalone _compute ────────────── + +TEST(CalculiXGate, StandaloneStressMatchesIsotropicLaw) { + for (auto const &e : sample_strains()) { + T2 const eps = strain_tensor_from_emec(e); + T2 sigma; + T4 tangent; + LinearElastic_compute(kLambda, kMu, eps, sigma, tangent); + T2 const expected = isotropic_stress(eps); + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR(sigma(i, j), expected(i, j), 1e-12) + << "stress mismatch at (" << i << "," << j << ")"; + } +} + +TEST(CalculiXGate, ConsistentTangentMatchesFiniteDifference) { + auto stress_only = [&](auto const &e) -> T2 { + T2 s; + T4 t; + LinearElastic_compute(kLambda, kMu, e, s, t); + return s; + }; + numsim::codegen::verify::NumericalTangentVerifier<3> const verifier( + {.abs_tol = 1e-7, .rel_tol = 1e-6, .fd_step = 1e-6}); + for (auto const &e : sample_strains()) { + T2 const eps = strain_tensor_from_emec(e); + T2 s; + T4 emitted_tangent; + LinearElastic_compute(kLambda, kMu, eps, s, emitted_tangent); + auto const r = verifier.verify(stress_only, eps, emitted_tangent); + EXPECT_TRUE(r.passed) << "tangent FD mismatch: max_abs=" << r.max_abs_dev + << " max_rel=" << r.max_rel_dev; + } +} + +// ── External NCG_UMAT ABI boundary + stiff packing ─────────────────────────── + +TEST(CalculiXGate, ExternalUmatStressMatchesIsotropicOracle) { + for (auto const &e : sample_strains()) { + std::array stre{}; + std::array stiff{}; + call_ext(e, stre, stiff); + T2 const sigma = isotropic_stress(strain_tensor_from_emec(e)); + // stre is abq_std order {11,22,33,12,13,23}. + EXPECT_NEAR(stre[0], sigma(0, 0), 1e-12); + EXPECT_NEAR(stre[1], sigma(1, 1), 1e-12); + EXPECT_NEAR(stre[2], sigma(2, 2), 1e-12); + EXPECT_NEAR(stre[3], sigma(0, 1), 1e-12); + EXPECT_NEAR(stre[4], sigma(0, 2), 1e-12); + EXPECT_NEAR(stre[5], sigma(1, 2), 1e-12); + } +} + +TEST(CalculiXGate, ExternalUmatStiffMatchesColumnMajorPackedOracle) { + auto const expected = expected_stiff_column_major(); + // Strain-independent here; use a general state so every column participates. + std::array stre{}; + std::array stiff{}; + call_ext(sample_strains().back(), stre, stiff); + for (std::size_t k = 0; k < 21; ++k) + EXPECT_NEAR(stiff[k], expected[k], 1e-12) << "stiff mismatch at index " << k; +} + +// Negative control: packing ORDER is observable — index 2 holds D(0,2)=λ +// column-major vs D(1,1)=λ+2μ row-major. stiff must match column-major only. +TEST(CalculiXGate, StiffPackingIsColumnMajorNotRowMajor) { + double const a = kLambda + 2.0 * kMu, b = kLambda; + double D[6][6] = {{a, b, b, 0, 0, 0}, {b, a, b, 0, 0, 0}, + {b, b, a, 0, 0, 0}, {0, 0, 0, kMu, 0, 0}, + {0, 0, 0, 0, kMu, 0}, {0, 0, 0, 0, 0, kMu}}; + std::array row_major{}; + std::size_t k = 0; + for (int i = 0; i < 6; ++i) + for (int j = i; j < 6; ++j) row_major[k++] = D[i][j]; + + std::array stre{}; + std::array stiff{}; + call_ext(sample_strains().back(), stre, stiff); + + // stiff equals the column-major oracle but NOT the row-major one. + auto const col_major = expected_stiff_column_major(); + bool differs_from_row_major = false; + for (std::size_t idx = 0; idx < 21; ++idx) { + EXPECT_NEAR(stiff[idx], col_major[idx], 1e-12); + if (std::abs(stiff[idx] - row_major[idx]) > 1e-9) differs_from_row_major = true; + } + EXPECT_TRUE(differs_from_row_major) + << "column-major and row-major packings are indistinguishable here — " + "the negative control has no discriminating power"; +} + +// Regression: MPROPS must be read on EVERY call. The thread_local material once +// cached the first call's constants, so a second call with different λ/μ +// silently returned the first answer. Two back-to-back calls must both be +// correct. (CalculiX legitimately varies the constants by temperature.) +TEST(CalculiXGate, ExternalReadsConstantsEveryCall) { + std::array const e{{0.01, 0.0, 0.0, 0.0, 0.0, 0.0}}; + std::array s1{}, s2{}; + std::array k1{}, k2{}; + call_ext(e, s1, k1, /*lambda=*/1.3, /*mu=*/0.7); + call_ext(e, s2, k2, /*lambda=*/0.5, /*mu=*/0.2); // same thread, new constants + // S11 = (λ+2μ)·0.01: 0.027 then 0.009 — must differ and both be correct. + EXPECT_NEAR(s1[0], (1.3 + 2 * 0.7) * 0.01, 1e-12); + EXPECT_NEAR(s2[0], (0.5 + 2 * 0.2) * 0.01, 1e-12); + EXPECT_GT(std::abs(s1[0] - s2[0]), 1e-6) << "constants not re-read per call"; +} + +// icmd==3 requests stress only: stress must still be correct, and stiff must be +// left untouched. The call helpers pre-fill stiff with 0.0, so an untouched +// stiff stays all-zero — whereas the icmd=1 path writes nonzero packed entries +// (e.g. λ+2μ). Assert both to make "untouched" meaningful. +TEST(CalculiXGate, StressOnlyIcmd3LeavesStiffUntouched) { + auto const e = sample_strains().back(); // copy: don't bind a ref into a temporary + T2 const sigma = isotropic_stress(strain_tensor_from_emec(e)); + std::array stre{}; + std::array stiff3{}, stiff1{}; + call_ext(e, stre, stiff3, kLambda, kMu, /*icmd=*/3); + call_ext(e, stre, stiff1, kLambda, kMu, /*icmd=*/1); + EXPECT_NEAR(stre[0], sigma(0, 0), 1e-12); // stress still correct + for (std::size_t k = 0; k < 21; ++k) + EXPECT_DOUBLE_EQ(stiff3[k], 0.0) << "stiff written on icmd==3 at " << k; + EXPECT_GT(stiff1[0], 0.0) << "icmd==1 must fill stiff (control)"; +} diff --git a/tests/generated/generate_calculix_check.cpp b/tests/generated/generate_calculix_check.cpp new file mode 100644 index 0000000..88812c5 --- /dev/null +++ b/tests/generated/generate_calculix_check.cpp @@ -0,0 +1,84 @@ +// Generator for the CalculiX end-to-end gate; CMake runs it at build time: +// argv[1] = LinearElastic.h (StandaloneCxx) +// argv[2] = LinearElastic_ext.cpp (CalculiXExternal plugin) +// Both come from the SAME isotropic-elastic recipe. calculix_check_driver.cpp +// FD-verifies the tangent and checks NCG_UMAT against an independent oracle. + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + +// Isotropic linear elasticity, Lamé form. λ is declared before μ: that is the +// order ccx reads *USER MATERIAL, CONSTANTS=2 into MPROPS[0], MPROPS[1]. +auto build_linear_elastic() -> numsim::codegen::ConstitutiveModel { + using namespace numsim::cas; + using namespace numsim::codegen; + + ConstitutiveModel m("LinearElastic"); + auto lambda = m.add_parameter("lambda", 1.0, "First Lame parameter"); + auto mu = m.add_parameter("mu", 0.5, "Shear modulus (second Lame parameter)"); + // roles::Strain makes eps symmetric, so diff yields a minor-symmetric tangent, + // which is what the stiff packing assumes. + auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain); + + auto I = make_expression(std::size_t{3}, std::size_t{2}); + auto tr = make_expression(eps); // t2s + // λ·(tr(ε)·I): tensor_to_scalar_with_tensor_mul(I, tr) = tr(ε)·I, then ·λ. + auto vol = lambda * make_expression(I, tr); + auto sigma = vol + 2 * mu * eps; + + m.add_output("stress", sigma, roles::Stress); + m.add_algorithmic_tangent("dstress_deps", "stress", "eps"); + return m; +} + +auto write_single_file(std::vector const &files, + char const *out_path) -> int { + if (files.size() != 1) { + std::cerr << "expected single emitted file, got " << files.size() << "\n"; + return 1; + } + std::ofstream out(out_path); + if (!out) { + std::cerr << "could not open '" << out_path << "' for writing\n"; + return 1; + } + out << files[0].contents; + return out ? 0 : 1; +} + +} // namespace + +int main(int argc, char *argv[]) { + if (argc != 3) { + std::cerr << "usage: " << argv[0] + << " \n"; + return 1; + } + + auto const model = build_linear_elastic(); + + if (int rc = write_single_file( + numsim::codegen::StandaloneCxxTarget{}.emit(model), argv[1]); + rc != 0) { + return rc; + } + if (int rc = write_single_file( + numsim::codegen::CalculiXExternalTarget{}.emit(model), argv[2]); + rc != 0) { + return rc; + } + return 0; +}