From 0004506d0d477a20626be5d249186b9d8c9cca66 Mon Sep 17 00:00:00 2001 From: petlenz Date: Wed, 12 Aug 2026 10:07:53 +0200 Subject: [PATCH 01/11] =?UTF-8?q?materials:=20isotropic=5Ftangent=20+=20li?= =?UTF-8?q?near=5Fstress=20=E2=80=94=20a=20stiffness=20the=20graph=20can?= =?UTF-8?q?=20order?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit linear_elasticity owns its tangent and builds it once in the constructor, which is optimal while K and G are fixed but leaves no way to recompute it — and, more awkwardly, no way to recompute it CORRECTLY even if one were added. A material reading its own property creates no edge in the property graph. The engine builds edges from input_dependencies and from material-level input wiring, so `elastic::stress` and `elastic::tangent` both hang off the same single edge from the strain source and their relative order falls out of constructor registration order: strain_in::strain elastic::stress <- reads m_C elastic::tangent <- writes m_C, too late Attaching an update callback to the tangent would therefore have produced a stress lagging the constants by one call, silently. Splitting the stiffness into its own material turns that invisible intra-material read into a real Global edge, which the topological sort honours: stiffness::tangent <- produced first elastic::stress That is the change. Same physics, different graph shape, and the ordering is now guaranteed by construction rather than by registration order. isotropic_tangent binds its update callback ONLY when the "recompute" parameter is true. add_output ignores a null callback, so with recompute=false the property carries no callback at all and the engine skips it outright — the per-call cost is zero, not merely small, and the behaviour matches linear_elasticity exactly. The tangent is built once in the constructor either way, so it is valid before the first update() regardless. recompute = false constants fixed (Abaqus: PROPS cannot vary per name) recompute = true constants vary between calls (CalculiX interpolates *USER MATERIAL constants by temperature) It defaults to false, matching linear_elasticity. The failure mode of the wrong setting is worth naming: recompute=false while the constants do move leaves a stale tangent with the stress still correct, so it costs convergence rate rather than accuracy and produces no diagnostic. A debug assertion catches it; RecomputeFalseIgnoresALaterParameterWrite pins the behaviour so it is recorded rather than discovered. linear_stress is linear_elasticity minus tangent ownership: sigma = C : eps with C as a Global input. Any material producing a "tangent" property is a drop-in, so anisotropic, temperature-dependent or damage-degraded stiffnesses can be substituted without consumers knowing. This is the decomposition isotropic_damage already uses (state function + yield + damage law + assembler). Nothing existing is modified. linear_elasticity is untouched and remains the better choice whenever the moduli are fixed — one fewer material in the graph and the tangent computed once. j2_plasticity needs no change either: its elastic_source only requires some material producing "tangent". Depends on the parent commit for set_parameter, which the write-and-retrigger test uses. recompute=true has no purpose without a way to write parameters, so the two are only meaningful together. --- .../materials/isotropic_tangent.h | 120 ++++++++++ .../materials/linear_stress.h | 64 +++++ tests/CMakeLists.txt | 1 + tests/test_tangent_generator.cpp | 220 ++++++++++++++++++ 4 files changed, 405 insertions(+) create mode 100644 include/numsim-materials/materials/isotropic_tangent.h create mode 100644 include/numsim-materials/materials/linear_stress.h create mode 100644 tests/test_tangent_generator.cpp diff --git a/include/numsim-materials/materials/isotropic_tangent.h b/include/numsim-materials/materials/isotropic_tangent.h new file mode 100644 index 0000000..a13d6b6 --- /dev/null +++ b/include/numsim-materials/materials/isotropic_tangent.h @@ -0,0 +1,120 @@ +#ifndef NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H +#define NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H + +#include +#include "numsim-materials/core/material_base.h" + +namespace numsim::materials { + +/// Isotropic elastic stiffness as a material of its own. +/// +/// `linear_elasticity` owns its tangent and computes it once in the +/// constructor, which is optimal while K and G are fixed but leaves no way to +/// recompute it — and no way to order that recomputation correctly if there +/// were one. A material reading its OWN property creates no edge in the +/// property graph, so `elastic::stress` sorts before `elastic::tangent` and +/// would consume a stale stiffness. Splitting the stiffness into a separate +/// material turns that invisible dependency into a real Global edge, which the +/// topological sort then honours. +/// +/// It also makes the stiffness pluggable: anything producing a "tangent" +/// property — anisotropic, temperature-dependent, damage-degraded — is a drop-in +/// replacement, and consumers need not know which. +/// +/// ### recompute +/// +/// The update callback is bound ONLY when the "recompute" parameter is true. +/// With it false the tangent is built once, at construction, and the property +/// carries no callback at all — so `ctx.update()` skips it outright and the +/// per-call cost is exactly zero, matching `linear_elasticity`. +/// +/// recompute = false constants are fixed (Abaqus: PROPS cannot vary for a +/// given material name) +/// recompute = true constants change between calls (CalculiX interpolates +/// *USER MATERIAL constants by temperature) +/// +/// Setting it false while the constants do move yields a stale tangent, with +/// the stress still correct — the classic silently-wrong tangent, costing +/// convergence rate rather than accuracy. A debug assertion catches it. +template +class isotropic_tangent final + : public material_base, Traits> { +public: + using base = material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + using base::Dim; + using tensor4 = tmech::tensor; + + template + explicit isotropic_tangent(Args&&... args) + : base(std::forward(args)...), + m_K(base::template get_parameter("K")), + m_G(base::template get_parameter("G")), + m_recompute(base::template get_parameter("recompute")), + // Bind the callback only when asked. add_output ignores a null one, so + // with recompute=false the property has no callback and the engine + // never visits it. + m_C(base::template add_output( + "tangent", + base::template get_parameter("recompute") + ? &isotropic_tangent::update_tangent + : nullptr)) { + // Always compute once, so the tangent is valid before the first update() + // whether or not it will ever be recomputed. + update_tangent(); +#ifndef NDEBUG + m_K_at_construction = m_K; + m_G_at_construction = m_G; +#endif + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + para.template insert("K").template add(); + para.template insert("G").template add(); + // Defaults to the cheap behaviour, which matches linear_elasticity and is + // right whenever the host cannot vary the constants. + para.template insert("recompute").template add(false); + return para; + } + + void update_tangent() { + const auto I{tmech::eye()}; + const auto IIsym{(tmech::otimesu(I, I) + tmech::otimesl(I, I)) * 0.5}; + const auto IIvol{tmech::otimes(I, I) / Dim}; + const auto IIdev{IIsym - IIvol}; + m_C = 3 * m_K * IIvol + 2 * m_G * IIdev; + } + + /// True when this material will follow a change to K or G. + [[nodiscard]] bool recomputes() const noexcept { return m_recompute; } + +#ifndef NDEBUG + /// Debug-only guard against the one way this can be configured wrongly: + /// recompute=false while the constants actually move. Callers that write + /// parameters may invoke this after writing; it is a no-op in release. + void assert_constants_unchanged() const { + assert((m_recompute || (m_K == m_K_at_construction && + m_G == m_G_at_construction)) && + "isotropic_tangent: K or G changed but recompute=false, so the " + "tangent is stale"); + } +#else + void assert_constants_unchanged() const noexcept {} +#endif + +private: + const value_type& m_K; + const value_type& m_G; + const bool& m_recompute; + tensor4& m_C; +#ifndef NDEBUG + value_type m_K_at_construction{}; + value_type m_G_at_construction{}; +#endif +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H diff --git a/include/numsim-materials/materials/linear_stress.h b/include/numsim-materials/materials/linear_stress.h new file mode 100644 index 0000000..92b81da --- /dev/null +++ b/include/numsim-materials/materials/linear_stress.h @@ -0,0 +1,64 @@ +#ifndef NUMSIM_MATERIALS_LINEAR_STRESS_H +#define NUMSIM_MATERIALS_LINEAR_STRESS_H + +#include +#include "numsim-materials/core/material_base.h" + +namespace numsim::materials { + +/// sigma = C : eps, with C supplied by another material. +/// +/// The decomposed counterpart to `linear_elasticity`: same physics, but the +/// stiffness arrives as a Global input rather than being owned and computed +/// in the constructor. That single change is what makes the ordering correct — +/// consuming another material's property creates a graph edge, so the engine +/// guarantees the stiffness is produced before this material reads it. A +/// material reading its own property has no such edge (see isotropic_tangent). +/// +/// Pair it with any tangent generator. `linear_elasticity` remains the better +/// choice whenever the moduli are fixed, since it computes the tangent once and +/// costs one less material in the graph; reach for this pair when the stiffness +/// has to follow something that changes. +template +class linear_stress final + : public material_base, Traits> { +public: + using base = material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + using base::Dim; + using tensor2 = tmech::tensor; + using tensor4 = tmech::tensor; + + template + explicit linear_stress(Args&&... args) + : base(std::forward(args)...), + m_sig(base::template add_output( + "stress", &linear_stress::update_stress)), + m_C(base::template add_input( + base::template get_parameter("tangent_source"), + "tangent", EdgeKind::Global)), + m_eps(base::template add_input( + base::template get_parameter("strain_source"), + "strain", EdgeKind::Global)) {} + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + para.template insert("tangent_source") + .template add(); + para.template insert("strain_source") + .template add(); + return para; + } + + void update_stress() { m_sig = tmech::dcontract(m_C.get(), m_eps.get()); } + +private: + tensor2& m_sig; + const input_property& m_C; + const input_property& m_eps; +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_LINEAR_STRESS_H diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8a7433d..d03b294 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -23,6 +23,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_tangent_generator test_tangent_generator.cpp) target_link_libraries(test_umat_interface PRIVATE Threads::Threads) # Data dumper for plotting (not a test — standalone executable) diff --git a/tests/test_tangent_generator.cpp b/tests/test_tangent_generator.cpp new file mode 100644 index 0000000..3be8941 --- /dev/null +++ b/tests/test_tangent_generator.cpp @@ -0,0 +1,220 @@ +#include +#include +#include +#include +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/materials/isotropic_tangent.h" +#include "numsim-materials/materials/linear_elasticity.h" +#include "numsim-materials/materials/linear_stress.h" +#include "numsim-materials/umat/external_state_source.h" + +namespace { + +namespace nm = numsim::materials; + +using policy = nm::material_policy_default; +using T = policy::value_type; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; +using tensor2 = tmech::tensor; +using tensor4 = tmech::tensor; + +constexpr T K = 166.67; +constexpr T G = 76.92; + +/// strain source -> isotropic_tangent -> linear_stress +nm::external_strain_source& build_decomposed(ctx_type& ctx, T k, T g, + bool recompute) { + param_type p; + p.insert("name", "strain_in"); + auto& src = ctx.create>(p); + + p.clear(); + p.insert("name", "stiffness"); + p.insert("K", k); + p.insert("G", g); + p.insert("recompute", recompute); + ctx.create>(p); + + p.clear(); + p.insert("name", "elastic"); + p.insert("tangent_source", "stiffness"); + p.insert("strain_source", "strain_in"); + ctx.create>(p); + + ctx.finalize(); + return src; +} + +tensor2 uniaxial(T v) { + tensor2 e; + e.fill(0.0); + e(0, 0) = v; + return e; +} + +// --------------------------------------------------------------------------- +// The ordering this decomposition exists to fix +// --------------------------------------------------------------------------- + +/// The whole point. Within one material, `stress` sorts BEFORE `tangent` +/// because a material reading its own property creates no graph edge — so +/// linear_elasticity could never recompute its tangent safely. Consuming +/// another material's property creates a real Global edge, and the topological +/// sort then puts the producer first. +TEST(TangentGenerator, StiffnessIsOrderedBeforeTheStressThatConsumesIt) { + ctx_type ctx; + build_decomposed(ctx, K, G, /*recompute=*/true); + + std::size_t i_tangent = 0, i_stress = 0, n = 0; + for (const auto* prop : ctx.property_execution_order()) { + const auto& id = prop->traits().id; + if (id.owner == "stiffness" && id.name == "tangent") i_tangent = n; + if (id.owner == "elastic" && id.name == "stress") i_stress = n; + ++n; + } + EXPECT_LT(i_tangent, i_stress) + << "the stiffness must be produced before the stress reads it"; +} + +// --------------------------------------------------------------------------- +// Equivalence with the monolithic material +// --------------------------------------------------------------------------- + +/// Decomposed and monolithic must agree exactly — same physics, different +/// graph shape. +TEST(TangentGenerator, MatchesLinearElasticityExactly) { + ctx_type dec; + auto& dec_src = build_decomposed(dec, K, G, /*recompute=*/false); + + ctx_type mono; + nm::external_strain_source* mono_src = nullptr; + { + param_type p; + p.insert("name", "strain_in"); + mono_src = &mono.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("strain_producer_name", "strain_in"); + p.insert("K", K); + p.insert("G", G); + mono.create>(p); + mono.finalize(); + } + + for (int step = 1; step <= 5; ++step) { + const auto eps = uniaxial(0.001 * step); + dec_src.bind(eps, eps); + mono_src->bind(eps, eps); + dec.update(); + mono.update(); + + const auto& a = dec.get("elastic", "stress"); + const auto& b = mono.get("elastic", "stress"); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + EXPECT_DOUBLE_EQ(a(i, j), b(i, j)) << "step " << step; + } + + // ... and the stiffness itself. + const auto& Cd = dec.get("stiffness", "tangent"); + const auto& Cm = mono.get("elastic", "tangent"); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + for (int k = 0; k < 3; ++k) + for (int l = 0; l < 3; ++l) + EXPECT_DOUBLE_EQ(Cd(i, j, k, l), Cm(i, j, k, l)); +} + +// --------------------------------------------------------------------------- +// recompute: bound vs unbound callback +// --------------------------------------------------------------------------- + +/// With recompute=false the property carries NO callback, so the engine skips +/// it entirely — the per-call cost is zero rather than merely small. +TEST(TangentGenerator, RecomputeFalseLeavesThePropertyWithoutACallback) { + ctx_type ctx; + build_decomposed(ctx, K, G, /*recompute=*/false); + + const auto* prop = ctx.find_property("stiffness", "tangent"); + ASSERT_NE(prop, nullptr); + EXPECT_FALSE(static_cast(prop->traits().update)) + << "recompute=false must not bind an update callback"; +} + +TEST(TangentGenerator, RecomputeTrueBindsTheCallback) { + ctx_type ctx; + build_decomposed(ctx, K, G, /*recompute=*/true); + + const auto* prop = ctx.find_property("stiffness", "tangent"); + ASSERT_NE(prop, nullptr); + EXPECT_TRUE(static_cast(prop->traits().update)); +} + +/// Even with no callback the tangent must be valid: it is built once in the +/// constructor, before any update() runs. +TEST(TangentGenerator, TangentIsValidBeforeTheFirstUpdate) { + ctx_type ctx; + build_decomposed(ctx, K, G, /*recompute=*/false); + + const auto& C = ctx.get("stiffness", "tangent"); + // C_1111 = K + 4G/3 + EXPECT_NEAR(C(0, 0, 0, 0), K + 4.0 * G / 3.0, 1e-9); +} + +/// The behaviour the whole design is for: writing new constants in place and +/// having the graph pick them up, with correct ordering, on the next update. +/// +/// Writing goes through the non-const get(), which mutates the value inside +/// the std::any in place. insert() would replace the whole any and, for a type +/// past its small-buffer, relocate the object — invalidating the reference the +/// material bound at construction. +TEST(TangentGenerator, RecomputeTrueFollowsAParameterWrittenInPlace) { + ctx_type ctx; + auto& src = build_decomposed(ctx, K, G, /*recompute=*/true); + + const auto eps = uniaxial(0.001); + src.bind(eps, eps); + ctx.update(); + const T before = ctx.get("elastic", "stress")(0, 0); + EXPECT_NEAR(before, (K + 4.0 * G / 3.0) * 0.001, 1e-12); + + // Double the moduli in place, as a props writer would. + auto* stiffness = ctx.find("stiffness"); + ASSERT_NE(stiffness, nullptr); + auto* typed = dynamic_cast*>(stiffness); + ASSERT_NE(typed, nullptr); + ASSERT_TRUE(typed->recomputes()); + typed->template set_parameter("K", 2 * K); + typed->template set_parameter("G", 2 * G); + + ctx.update(); + const T after = ctx.get("elastic", "stress")(0, 0); + EXPECT_NEAR(after, 2 * (K + 4.0 * G / 3.0) * 0.001, 1e-12); + EXPECT_NEAR(after, 2 * before, 1e-12); +} + +/// The counterpart: with recompute=false the write is visible in the parameter +/// but the tangent does NOT follow it. This is the documented trade, and the +/// test exists so the behaviour is pinned rather than discovered. +TEST(TangentGenerator, RecomputeFalseIgnoresALaterParameterWrite) { + ctx_type ctx; + auto& src = build_decomposed(ctx, K, G, /*recompute=*/false); + + const auto eps = uniaxial(0.001); + src.bind(eps, eps); + ctx.update(); + const T before = ctx.get("elastic", "stress")(0, 0); + + auto* typed = + dynamic_cast*>(ctx.find("stiffness")); + ASSERT_NE(typed, nullptr); + EXPECT_FALSE(typed->recomputes()); + typed->template set_parameter("K", 2 * K); + + ctx.update(); + EXPECT_DOUBLE_EQ(ctx.get("elastic", "stress")(0, 0), before) + << "recompute=false must leave the tangent as built"; +} + +} // namespace From e1889d8a4276b305579b1f8d534ee9df7cf73f61 Mon Sep 17 00:00:00 2001 From: petlenz Date: Thu, 13 Aug 2026 21:48:59 +0200 Subject: [PATCH 02/11] =?UTF-8?q?materials:=20address=20review=20=E2=80=94?= =?UTF-8?q?=20flag=20honesty,=20composition=20with=20existing=20consumers?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Six findings from the review of #26. recomputes() could lie. "recompute" was held as a `const bool&` into the parameter store while the callback binding is frozen at construction, so flipping the flag with set_parameter made the accessor report that the tangent tracks K and G when no callback existed and it was permanently stale. It is now read once and stored BY VALUE, so the accessor describes what was actually bound. RecomputeIsReadOnceAndTheAccessorCannotLie pins it, and fails if the member is changed back to a reference. assert_constants_unchanged() had no callers anywhere, so the "a debug assertion catches it" promise in the header was false. Removed, along with the claim. The header now says plainly that nothing detects the stale-tangent misconfiguration automatically — the material cannot observe a write it has no callback for — and points at the test that pins what happens instead. Deleting it also removes the NDEBUG-conditional data members, which gave a header-only class template two layouts across mixed-config translation units, and the use of assert without including . The ordering test could pass on a graph missing the property it orders. Both indices used 0 as a not-found sentinel, and the tangent has no inputs so index 0 is exactly where it legitimately lands. Now std::optional with ASSERT_TRUE on both; verified by renaming the searched property, which the old version tolerated and the new one rejects. linear_stress could not substitute for linear_elasticity in existing consumers, which contradicted the header's drop-in claim. Measured, the blast radius was smaller than it looked: small_strain_plasticity and rk_plasticity already take their tangent from a separately NAMED source, so they work with the decomposed pair unchanged — a J2 chain with isotropic_tangent replacing linear_elasticity outright produces bit-identical results. Only three consumers resolved both stress and tangent from one name, and each now takes an optional "tangent_source" that defaults to the existing source: material_point_evaluator (config field), isotropic_damage (parameter), weighted_sum (a parallel list, one entry per term, since its sources are pairs). No existing configuration or behaviour changes; this widens an interface rather than refactoring one. DrivesJ2PlasticityIdenticallyToLinearElasticity is the test that should have existed from the start. The previous coverage compared two contexts built in the same file, which is exactly why a claim about pre-existing consumers went unchecked. Both materials are now registered in register_default_materials(), so they are reachable from the runtime factory rather than only from hand-written C++. bool is already in the JSON type registry, so "recompute" converts. update_tangent() reuses plasticity_detail::make_IIdev instead of spelling out the isotropic basis a fourth time. 193/193 tests pass. --- include/numsim-materials/default_materials.h | 4 + .../materials/isotropic_damage.h | 10 +- .../materials/isotropic_tangent.h | 45 +++---- .../numsim-materials/materials/weighted_sum.h | 20 ++- .../umat/material_point_evaluator.h | 13 +- tests/test_tangent_generator.cpp | 118 +++++++++++++++++- 6 files changed, 173 insertions(+), 37 deletions(-) diff --git a/include/numsim-materials/default_materials.h b/include/numsim-materials/default_materials.h index 45aed64..dbc56e9 100644 --- a/include/numsim-materials/default_materials.h +++ b/include/numsim-materials/default_materials.h @@ -6,7 +6,9 @@ #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/solvers/vector_newton.h" #include "numsim-materials/materials/scalar_stepper.h" +#include "numsim-materials/materials/isotropic_tangent.h" #include "numsim-materials/materials/linear_elasticity.h" +#include "numsim-materials/materials/linear_stress.h" #include "numsim-materials/materials/autocatalytic_reaction.h" #include "numsim-materials/materials/tensor_component_stepper.h" #include "numsim-materials/materials/scalar_identity_weight.h" @@ -69,6 +71,8 @@ void register_default_materials() { auto& factory = material_factory::instance(); factory.template register_type>("scalar_stepper"); factory.template register_type>("linear_elasticity"); + factory.template register_type>("isotropic_tangent"); + factory.template register_type>("linear_stress"); factory.template register_type>("autocatalytic_reaction"); factory.template register_type>("backward_euler"); factory.template register_type>("tensor_component_stepper_rank1"); diff --git a/include/numsim-materials/materials/isotropic_damage.h b/include/numsim-materials/materials/isotropic_damage.h index 889da36..1387c9f 100644 --- a/include/numsim-materials/materials/isotropic_damage.h +++ b/include/numsim-materials/materials/isotropic_damage.h @@ -50,8 +50,14 @@ class isotropic_damage final // inputs m_stress(base::template add_input( m_elastic_source, "stress", EdgeKind::Global)), + // Empty tangent_source means "same material as the stress", which is + // how a monolithic linear_elasticity is wired. Naming it separately + // allows a decomposed stiffness (isotropic_tangent + linear_stress). m_tangent(base::template add_input( - m_elastic_source, "tangent", EdgeKind::Global)), + base::template get_parameter("tangent_source").empty() + ? m_elastic_source + : base::template get_parameter("tangent_source"), + "tangent", EdgeKind::Global)), m_damage(base::template add_input( m_damage_source, "damage", EdgeKind::Global)), m_d_damage(base::template add_input( @@ -65,6 +71,8 @@ class isotropic_damage final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; para.template insert("elastic_source").template add(); + para.template insert("tangent_source") + .template add(std::string{}); para.template insert("damage_source").template add(); para.template insert("state_source").template add(); para.template insert("yield_source").template add(); diff --git a/include/numsim-materials/materials/isotropic_tangent.h b/include/numsim-materials/materials/isotropic_tangent.h index a13d6b6..cdcdb70 100644 --- a/include/numsim-materials/materials/isotropic_tangent.h +++ b/include/numsim-materials/materials/isotropic_tangent.h @@ -3,6 +3,7 @@ #include #include "numsim-materials/core/material_base.h" +#include "numsim-materials/materials/plasticity_utils.h" namespace numsim::materials { @@ -33,9 +34,19 @@ namespace numsim::materials { /// recompute = true constants change between calls (CalculiX interpolates /// *USER MATERIAL constants by temperature) /// +/// The flag is read ONCE, at construction, because that is when the callback is +/// bound — there is no way to attach one later. It is therefore stored by value +/// rather than as a reference into the parameter store: writing "recompute" +/// afterwards changes nothing, and an accessor reading the live parameter would +/// report that the tangent tracks K and G when in fact no callback exists. +/// `recomputes()` reports what was actually bound. +/// /// Setting it false while the constants do move yields a stale tangent, with /// the stress still correct — the classic silently-wrong tangent, costing -/// convergence rate rather than accuracy. A debug assertion catches it. +/// convergence rate rather than accuracy. Nothing detects that automatically: +/// the material cannot see a write it has no callback to observe. Choosing the +/// flag correctly is the caller's responsibility, and +/// RecomputeFalseIgnoresALaterParameterWrite pins what happens if they do not. template class isotropic_tangent final : public material_base, Traits> { @@ -51,6 +62,7 @@ class isotropic_tangent final : base(std::forward(args)...), m_K(base::template get_parameter("K")), m_G(base::template get_parameter("G")), + // By value, not by reference: see the class comment. m_recompute(base::template get_parameter("recompute")), // Bind the callback only when asked. add_output ignores a null one, so // with recompute=false the property has no callback and the engine @@ -63,10 +75,6 @@ class isotropic_tangent final // Always compute once, so the tangent is valid before the first update() // whether or not it will ever be recomputed. update_tangent(); -#ifndef NDEBUG - m_K_at_construction = m_K; - m_G_at_construction = m_G; -#endif } static input_parameter_controller parameters() { @@ -80,39 +88,22 @@ class isotropic_tangent final } void update_tangent() { + // IIdev comes from plasticity_utils rather than being spelled out again; + // the isotropic basis was already written three times in this repo. const auto I{tmech::eye()}; - const auto IIsym{(tmech::otimesu(I, I) + tmech::otimesl(I, I)) * 0.5}; const auto IIvol{tmech::otimes(I, I) / Dim}; - const auto IIdev{IIsym - IIvol}; - m_C = 3 * m_K * IIvol + 2 * m_G * IIdev; + m_C = 3 * m_K * IIvol + + 2 * m_G * plasticity_detail::make_IIdev(); } /// True when this material will follow a change to K or G. [[nodiscard]] bool recomputes() const noexcept { return m_recompute; } -#ifndef NDEBUG - /// Debug-only guard against the one way this can be configured wrongly: - /// recompute=false while the constants actually move. Callers that write - /// parameters may invoke this after writing; it is a no-op in release. - void assert_constants_unchanged() const { - assert((m_recompute || (m_K == m_K_at_construction && - m_G == m_G_at_construction)) && - "isotropic_tangent: K or G changed but recompute=false, so the " - "tangent is stale"); - } -#else - void assert_constants_unchanged() const noexcept {} -#endif - private: const value_type& m_K; const value_type& m_G; - const bool& m_recompute; + const bool m_recompute; tensor4& m_C; -#ifndef NDEBUG - value_type m_K_at_construction{}; - value_type m_G_at_construction{}; -#endif }; } // namespace numsim::materials diff --git a/include/numsim-materials/materials/weighted_sum.h b/include/numsim-materials/materials/weighted_sum.h index 3c19b29..93119d5 100644 --- a/include/numsim-materials/materials/weighted_sum.h +++ b/include/numsim-materials/materials/weighted_sum.h @@ -38,20 +38,33 @@ class weighted_sum final m_terms_param(base::template get_parameter("terms")), m_weight_property(base::template get_parameter("weight_property")), m_stress_property(base::template get_parameter("stress_property")), - m_tangent_property(base::template get_parameter("tangent_property")) + m_tangent_property(base::template get_parameter("tangent_property")), + m_tangent_sources(base::template get_parameter>("tangent_sources")) { // Dynamically create inputs for each term + std::size_t i = 0; for (const auto& [weight_name, mat_name] : m_terms_param) { + // A term's tangent normally comes from the same material as its stress. + // An entry in "tangent_sources" overrides that for one term, so a + // constituent whose stiffness lives in its own material (isotropic_tangent + // + linear_stress) can participate. A short or empty list leaves every + // unlisted term wired exactly as before. + const std::string& tangent_owner = + (i < m_tangent_sources.size() && !m_tangent_sources[i].empty()) + ? m_tangent_sources[i] + : mat_name; auto& w = base::template add_input(weight_name, m_weight_property, EdgeKind::Global); auto& s = base::template add_input(mat_name, m_stress_property, EdgeKind::Global); - auto& c = base::template add_input(mat_name, m_tangent_property, EdgeKind::Global); + auto& c = base::template add_input(tangent_owner, m_tangent_property, EdgeKind::Global); m_terms.push_back({&w, &s, &c}); + ++i; } - } static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; + para.template insert>("tangent_sources") + .template add(std::vector{}); para.template insert("terms").template add(); para.template insert("weight_property") .template add(std::string{"value"}); @@ -94,6 +107,7 @@ class weighted_sum final const std::string& m_weight_property; const std::string& m_stress_property; const std::string& m_tangent_property; + const std::vector& m_tangent_sources; std::vector m_terms; }; diff --git a/include/numsim-materials/umat/material_point_evaluator.h b/include/numsim-materials/umat/material_point_evaluator.h index bb377d3..6084ad7 100644 --- a/include/numsim-materials/umat/material_point_evaluator.h +++ b/include/numsim-materials/umat/material_point_evaluator.h @@ -43,8 +43,13 @@ class material_point_evaluator { struct config { /// Name of the external_strain_source material. std::string strain_source; - /// Material producing the stress and tangent the host wants back. + /// Material producing the stress the host wants back. std::string stress_source; + /// Material producing the tangent. Empty means "same as stress_source", + /// which is how every monolithic material (linear_elasticity, j2_plasticity) + /// is configured. Set it when the stiffness lives in its own material, as + /// with isotropic_tangent + linear_stress. + std::string tangent_source{}; std::string stress_property{"stress"}; std::string tangent_property{"tangent"}; /// Optional external_scalar_source carrying time; empty to omit. @@ -111,8 +116,10 @@ class material_point_evaluator { m_stress = resolve_property(m_cfg.stress_source, m_cfg.stress_property); - m_tangent = - resolve_property(m_cfg.stress_source, m_cfg.tangent_property); + const std::string& tangent_owner = m_cfg.tangent_source.empty() + ? m_cfg.stress_source + : m_cfg.tangent_source; + m_tangent = resolve_property(tangent_owner, m_cfg.tangent_property); if (!m_cfg.plastic_strain_property.empty()) { const auto src = connection_source::parse(m_cfg.plastic_strain_property); diff --git a/tests/test_tangent_generator.cpp b/tests/test_tangent_generator.cpp index 3be8941..be75114 100644 --- a/tests/test_tangent_generator.cpp +++ b/tests/test_tangent_generator.cpp @@ -1,11 +1,14 @@ #include -#include +#include #include #include #include "numsim-materials/core/material_context.h" #include "numsim-materials/materials/isotropic_tangent.h" #include "numsim-materials/materials/linear_elasticity.h" +#include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/linear_stress.h" +#include "numsim-materials/materials/small_strain_plasticity.h" +#include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/umat/external_state_source.h" namespace { @@ -66,14 +69,21 @@ TEST(TangentGenerator, StiffnessIsOrderedBeforeTheStressThatConsumesIt) { ctx_type ctx; build_decomposed(ctx, K, G, /*recompute=*/true); - std::size_t i_tangent = 0, i_stress = 0, n = 0; + // Not a sentinel of 0: the tangent has no inputs, so index 0 is exactly where + // it legitimately lands, and a not-found sentinel of 0 would be + // indistinguishable from the correct answer. The test would then pass on a + // graph that no longer contains the property being ordered. + std::optional i_tangent, i_stress; + std::size_t n = 0; for (const auto* prop : ctx.property_execution_order()) { const auto& id = prop->traits().id; if (id.owner == "stiffness" && id.name == "tangent") i_tangent = n; if (id.owner == "elastic" && id.name == "stress") i_stress = n; ++n; } - EXPECT_LT(i_tangent, i_stress) + ASSERT_TRUE(i_tangent.has_value()) << "stiffness::tangent is not in the graph"; + ASSERT_TRUE(i_stress.has_value()) << "elastic::stress is not in the graph"; + EXPECT_LT(*i_tangent, *i_stress) << "the stiffness must be produced before the stress reads it"; } @@ -217,4 +227,106 @@ TEST(TangentGenerator, RecomputeFalseIgnoresALaterParameterWrite) { << "recompute=false must leave the tangent as built"; } + +// --------------------------------------------------------------------------- +// The flag is construction-time only +// --------------------------------------------------------------------------- + +/// The callback can only be bound at construction, so "recompute" is read once +/// and stored BY VALUE. Reading the live parameter instead would let +/// recomputes() claim the tangent tracks K and G after someone flipped the flag +/// with set_parameter — while no callback exists and the tangent is permanently +/// stale. The accessor must describe what was actually bound. +TEST(TangentGenerator, RecomputeIsReadOnceAndTheAccessorCannotLie) { + ctx_type ctx; + build_decomposed(ctx, K, G, /*recompute=*/false); + + auto* typed = + dynamic_cast*>(ctx.find("stiffness")); + ASSERT_NE(typed, nullptr); + ASSERT_FALSE(typed->recomputes()); + + // Flipping the parameter afterwards cannot bind a callback... + typed->template set_parameter("recompute", true); + + // ... so the accessor must still report false, + EXPECT_FALSE(typed->recomputes()) + << "recomputes() must report what was bound, not the live parameter"; + // ... and the property must still carry no callback. + const auto* prop = ctx.find_property("stiffness", "tangent"); + ASSERT_NE(prop, nullptr); + EXPECT_FALSE(static_cast(prop->traits().update)); +} + +// --------------------------------------------------------------------------- +// Composition with a pre-existing consumer +// --------------------------------------------------------------------------- + +/// The header claims a tangent producer is a drop-in for any consumer. That has +/// to be tested against a consumer that already existed, not only against two +/// contexts built here. small_strain_plasticity takes its tangent from a NAMED +/// source, so it needs no change: point elastic_source at the generator. +TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { + auto drive = [](bool decomposed) { + ctx_type ctx; + param_type p; + p.insert("name", "strain_in"); + auto& src = ctx.create>(p); + + p.clear(); + p.insert("name", "stiffness"); + if (decomposed) { + p.insert("K", K); + p.insert("G", G); + p.insert("recompute", false); + ctx.create>(p); + } else { + p.insert("strain_producer_name", "strain_in"); + 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", T{1000}); + ctx.create>(p); + p.clear(); + p.insert("name", "j2"); + p.insert("elastic_source", "stiffness"); + p.insert("hardening_source", "hardening"); + p.insert("strain_source", "strain_in"); + p.insert("solver_source", "solver"); + p.insert("G", G); + p.insert("sigma_0", T{50}); + ctx.create>(p); + ctx.finalize(); + + std::vector out; + for (int step = 1; step <= 30; ++step) { + const auto eps = uniaxial(0.02 * step); // enough to pass sigma_0 = 50 + src.bind(eps, eps); + ctx.update(); + out.push_back(ctx.get("j2", "stress")(0, 0)); + out.push_back(ctx.get("j2", "equivalent_plastic_strain")); + ctx.commit(); + } + return out; + }; + + const auto with_generator = drive(true); + const auto with_monolith = drive(false); + ASSERT_EQ(with_generator.size(), with_monolith.size()); + bool went_plastic = false; + for (std::size_t i = 0; i < with_generator.size(); ++i) { + EXPECT_DOUBLE_EQ(with_generator[i], with_monolith[i]) << "sample " << i; + if (i % 2 == 1 && with_generator[i] > 1e-8) went_plastic = true; + } + EXPECT_TRUE(went_plastic) << "the path must yield for this to mean anything"; +} + } // namespace From 8a865fc485b7bb83f7f08905d902411324c18007 Mon Sep 17 00:00:00 2001 From: petlenz Date: Thu, 13 Aug 2026 21:58:27 +0200 Subject: [PATCH 03/11] materials: declare tangent_source optional rather than defaulting it to empty MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit There is no is_optional trait in this framework — the available checks are is_required, check_range, set_default, check_data_type and check_enum. But check_parameter only runs the checks that were registered, and the JSON visitor's accept() returns early for a key the input does not contain, so a parameter DECLARED WITH NO CHECK is already optional. That is how to spell it: para.template insert("tangent_source"); // optional paired with a contains() test at the use site, which is the idiom small_strain_plasticity already uses for "yield_function". This replaces set_default(std::string{}) plus an empty-string test. Two reasons it is better. Absent and empty become distinguishable: a deck that supplies tangent_source="" now fails to wire instead of silently falling back to the stress source, which is what a user typing an empty value would want to hear about. And unlike small_strain_plasticity's "yield_function" — which is probed with contains() but never declared — the key stays in the schema, so the JSON layer knows it rather than warning about an unrecognised key. material_point_evaluator's config is C++ rather than a parameter handler, so it uses std::optional for the same distinction. weighted_sum's per-term list follows the same rule: absent means every term takes its tangent from the material that produces its stress. Three tests pin the semantics, including the configuration the review showed was impossible before: the decomposed pair driving the UMAT evaluator, producing stress and DDSDDE identical to linear_elasticity over ten increments. The evaluator still throws when the tangent cannot be resolved at all, so omitting tangent_source with a stress-only material is a loud failure rather than a mystery. 193 -> 196 tests. --- .../materials/isotropic_damage.h | 20 ++-- .../numsim-materials/materials/weighted_sum.h | 17 ++-- .../umat/material_point_evaluator.h | 14 ++- tests/test_tangent_generator.cpp | 96 +++++++++++++++++++ 4 files changed, 128 insertions(+), 19 deletions(-) diff --git a/include/numsim-materials/materials/isotropic_damage.h b/include/numsim-materials/materials/isotropic_damage.h index 1387c9f..56609a3 100644 --- a/include/numsim-materials/materials/isotropic_damage.h +++ b/include/numsim-materials/materials/isotropic_damage.h @@ -50,13 +50,16 @@ class isotropic_damage final // inputs m_stress(base::template add_input( m_elastic_source, "stress", EdgeKind::Global)), - // Empty tangent_source means "same material as the stress", which is - // how a monolithic linear_elasticity is wired. Naming it separately + // An ABSENT tangent_source means "same material as the stress", which + // is how a monolithic linear_elasticity is wired. Naming it separately // allows a decomposed stiffness (isotropic_tangent + linear_stress). + // Tested with contains() rather than against an empty string: absent + // and empty are then distinct, so a deck that supplies "" gets a + // wiring error instead of silently falling back. m_tangent(base::template add_input( - base::template get_parameter("tangent_source").empty() - ? m_elastic_source - : base::template get_parameter("tangent_source"), + base::m_parameter_handler.contains("tangent_source") + ? base::template get_parameter("tangent_source") + : m_elastic_source, "tangent", EdgeKind::Global)), m_damage(base::template add_input( m_damage_source, "damage", EdgeKind::Global)), @@ -71,8 +74,11 @@ class isotropic_damage final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; para.template insert("elastic_source").template add(); - para.template insert("tangent_source") - .template add(std::string{}); + // Declared with no check, which is how this framework spells "optional": + // check_parameter only runs registered checks, and the JSON visitor skips + // keys the input does not contain. Declaring it anyway keeps it in the + // schema, so the JSON layer knows the key instead of warning about it. + para.template insert("tangent_source"); para.template insert("damage_source").template add(); para.template insert("state_source").template add(); para.template insert("yield_source").template add(); diff --git a/include/numsim-materials/materials/weighted_sum.h b/include/numsim-materials/materials/weighted_sum.h index 93119d5..d97730d 100644 --- a/include/numsim-materials/materials/weighted_sum.h +++ b/include/numsim-materials/materials/weighted_sum.h @@ -39,7 +39,10 @@ class weighted_sum final m_weight_property(base::template get_parameter("weight_property")), m_stress_property(base::template get_parameter("stress_property")), m_tangent_property(base::template get_parameter("tangent_property")), - m_tangent_sources(base::template get_parameter>("tangent_sources")) + m_has_tangent_sources(base::m_parameter_handler.contains("tangent_sources")), + m_tangent_sources(m_has_tangent_sources + ? base::template get_parameter>("tangent_sources") + : std::vector{}) { // Dynamically create inputs for each term std::size_t i = 0; @@ -50,9 +53,7 @@ class weighted_sum final // + linear_stress) can participate. A short or empty list leaves every // unlisted term wired exactly as before. const std::string& tangent_owner = - (i < m_tangent_sources.size() && !m_tangent_sources[i].empty()) - ? m_tangent_sources[i] - : mat_name; + (i < m_tangent_sources.size()) ? m_tangent_sources[i] : mat_name; auto& w = base::template add_input(weight_name, m_weight_property, EdgeKind::Global); auto& s = base::template add_input(mat_name, m_stress_property, EdgeKind::Global); auto& c = base::template add_input(tangent_owner, m_tangent_property, EdgeKind::Global); @@ -63,8 +64,9 @@ class weighted_sum final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert>("tangent_sources") - .template add(std::vector{}); + // Optional: declared with no check. Absent means every term takes its + // tangent from the same material as its stress. + para.template insert>("tangent_sources"); para.template insert("terms").template add(); para.template insert("weight_property") .template add(std::string{"value"}); @@ -107,7 +109,8 @@ class weighted_sum final const std::string& m_weight_property; const std::string& m_stress_property; const std::string& m_tangent_property; - const std::vector& m_tangent_sources; + const bool m_has_tangent_sources; + const std::vector m_tangent_sources; std::vector m_terms; }; diff --git a/include/numsim-materials/umat/material_point_evaluator.h b/include/numsim-materials/umat/material_point_evaluator.h index 6084ad7..45805e5 100644 --- a/include/numsim-materials/umat/material_point_evaluator.h +++ b/include/numsim-materials/umat/material_point_evaluator.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -45,11 +46,15 @@ class material_point_evaluator { std::string strain_source; /// Material producing the stress the host wants back. std::string stress_source; - /// Material producing the tangent. Empty means "same as stress_source", + /// Material producing the tangent. Unset means "same as stress_source", /// which is how every monolithic material (linear_elasticity, j2_plasticity) /// is configured. Set it when the stiffness lives in its own material, as /// with isotropic_tangent + linear_stress. - std::string tangent_source{}; + /// + /// std::optional rather than an empty string, so "not configured" and + /// "configured as empty" are distinct — the latter is a mistake and should + /// fail to resolve rather than silently fall back. + std::optional tangent_source{}; std::string stress_property{"stress"}; std::string tangent_property{"tangent"}; /// Optional external_scalar_source carrying time; empty to omit. @@ -116,9 +121,8 @@ class material_point_evaluator { m_stress = resolve_property(m_cfg.stress_source, m_cfg.stress_property); - const std::string& tangent_owner = m_cfg.tangent_source.empty() - ? m_cfg.stress_source - : m_cfg.tangent_source; + const std::string& tangent_owner = + m_cfg.tangent_source ? *m_cfg.tangent_source : m_cfg.stress_source; m_tangent = resolve_property(tangent_owner, m_cfg.tangent_property); if (!m_cfg.plastic_strain_property.empty()) { diff --git a/tests/test_tangent_generator.cpp b/tests/test_tangent_generator.cpp index be75114..e92f61b 100644 --- a/tests/test_tangent_generator.cpp +++ b/tests/test_tangent_generator.cpp @@ -10,10 +10,12 @@ #include "numsim-materials/materials/small_strain_plasticity.h" #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/umat/external_state_source.h" +#include "numsim-materials/umat/material_point_evaluator.h" namespace { namespace nm = numsim::materials; +namespace u = numsim::materials::umat; using policy = nm::material_policy_default; using T = policy::value_type; @@ -329,4 +331,98 @@ TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { EXPECT_TRUE(went_plastic) << "the path must yield for this to mean anything"; } + +// --------------------------------------------------------------------------- +// The optional tangent_source +// --------------------------------------------------------------------------- + +/// "Optional" in this framework means declared with no check: check_parameter +/// only runs registered checks, and the JSON visitor skips absent keys. These +/// pin that an ABSENT tangent_source falls back to the stress source, while a +/// SUPPLIED one is honoured — the two must be distinguishable, which an +/// empty-string sentinel could not express. +TEST(TangentSource, AbsentFallsBackToTheStressSource) { + ctx_type ctx; + param_type p; + p.insert("name", "strain_in"); + auto& src = ctx.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("strain_producer_name", "strain_in"); + p.insert("K", K); + p.insert("G", G); + ctx.create>(p); + ctx.finalize(); + + // No tangent_source configured: the monolithic material supplies both. + u::material_point_evaluator::config cfg; + cfg.strain_source = "strain_in"; + cfg.stress_source = "elastic"; + ASSERT_FALSE(cfg.tangent_source.has_value()); + EXPECT_NO_THROW(u::material_point_evaluator(ctx, cfg)); + + (void)src; +} + +TEST(TangentSource, SuppliedResolvesTheTangentElsewhere) { + ctx_type ctx; + build_decomposed(ctx, K, G, /*recompute=*/false); + + u::material_point_evaluator::config cfg; + cfg.strain_source = "strain_in"; + cfg.stress_source = "elastic"; // linear_stress: publishes only "stress" + // Without this the evaluator cannot find a tangent at all. + EXPECT_THROW(u::material_point_evaluator(ctx, cfg), u::fatal_error); + + cfg.tangent_source = "stiffness"; // the generator + EXPECT_NO_THROW(u::material_point_evaluator(ctx, cfg)); +} + +/// The decomposed pair driving a UMAT end to end — the configuration the review +/// showed was impossible before tangent_source existed. +TEST(TangentSource, DecomposedPairDrivesTheEvaluatorLikeLinearElasticity) { + ctx_type dec; + build_decomposed(dec, K, G, /*recompute=*/false); + u::material_point_evaluator::config dcfg; + dcfg.strain_source = "strain_in"; + dcfg.stress_source = "elastic"; + dcfg.tangent_source = "stiffness"; + u::material_point_evaluator deval(dec, dcfg); + + ctx_type mono; + { + param_type p; + p.insert("name", "strain_in"); + mono.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("strain_producer_name", "strain_in"); + p.insert("K", K); + p.insert("G", G); + mono.create>(p); + mono.finalize(); + } + u::material_point_evaluator::config mcfg; + mcfg.strain_source = "strain_in"; + mcfg.stress_source = "elastic"; + u::material_point_evaluator meval(mono, mcfg); + + std::vector dsv(deval.nstatv(), 0.0), msv(meval.nstatv(), 0.0); + T stran[6] = {0, 0, 0, 0, 0, 0}; + const T dstran[6] = {0.002, -0.0005, 0.0, 0.001, 0.0, 0.0}; + T ds[6], dd[36], ms[6], md[36]; + + for (int step = 0; step < 10; ++step) { + deval.evaluate({.stran = stran, .dstran = dstran, .stress = ds, + .ddsdde = dd, .statev = dsv}); + meval.evaluate({.stran = stran, .dstran = dstran, .stress = ms, + .ddsdde = md, .statev = msv}); + for (std::size_t i = 0; i < 6; ++i) + EXPECT_DOUBLE_EQ(ds[i], ms[i]) << "step " << step << " stress " << i; + for (std::size_t i = 0; i < 36; ++i) + EXPECT_DOUBLE_EQ(dd[i], md[i]) << "step " << step << " ddsdde " << i; + for (std::size_t i = 0; i < 6; ++i) stran[i] += dstran[i]; + } +} + } // namespace From e49790eb72dff7598a439ebae4109df0e803345b Mon Sep 17 00:00:00 2001 From: petlenz Date: Thu, 13 Aug 2026 22:01:52 +0200 Subject: [PATCH 04/11] materials: shorten the comments on the tangent generator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Roughly halved: 132 comment lines down to 70 across the new headers and tests. Kept the facts that took work to establish — no graph edge for an intra-material read, recompute is read once because that is when the callback is bound, the optional-by-no-check idiom, the 0-sentinel hazard in the ordering test — and cut the narrative around them. --- .../materials/isotropic_damage.h | 13 +--- .../materials/isotropic_tangent.h | 62 +++++------------- .../materials/linear_stress.h | 15 ++--- .../numsim-materials/materials/weighted_sum.h | 9 +-- .../umat/material_point_evaluator.h | 11 +--- tests/test_tangent_generator.cpp | 65 +++++-------------- 6 files changed, 48 insertions(+), 127 deletions(-) diff --git a/include/numsim-materials/materials/isotropic_damage.h b/include/numsim-materials/materials/isotropic_damage.h index 56609a3..d551ec5 100644 --- a/include/numsim-materials/materials/isotropic_damage.h +++ b/include/numsim-materials/materials/isotropic_damage.h @@ -50,12 +50,8 @@ class isotropic_damage final // inputs m_stress(base::template add_input( m_elastic_source, "stress", EdgeKind::Global)), - // An ABSENT tangent_source means "same material as the stress", which - // is how a monolithic linear_elasticity is wired. Naming it separately - // allows a decomposed stiffness (isotropic_tangent + linear_stress). - // Tested with contains() rather than against an empty string: absent - // and empty are then distinct, so a deck that supplies "" gets a - // wiring error instead of silently falling back. + // Absent tangent_source: same material as the stress. contains() + // rather than an empty-string test, so "" is an error, not a fallback. m_tangent(base::template add_input( base::m_parameter_handler.contains("tangent_source") ? base::template get_parameter("tangent_source") @@ -74,10 +70,7 @@ class isotropic_damage final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; para.template insert("elastic_source").template add(); - // Declared with no check, which is how this framework spells "optional": - // check_parameter only runs registered checks, and the JSON visitor skips - // keys the input does not contain. Declaring it anyway keeps it in the - // schema, so the JSON layer knows the key instead of warning about it. + // No check == optional; declared so the JSON schema still knows the key. para.template insert("tangent_source"); para.template insert("damage_source").template add(); para.template insert("state_source").template add(); diff --git a/include/numsim-materials/materials/isotropic_tangent.h b/include/numsim-materials/materials/isotropic_tangent.h index cdcdb70..82c7965 100644 --- a/include/numsim-materials/materials/isotropic_tangent.h +++ b/include/numsim-materials/materials/isotropic_tangent.h @@ -9,44 +9,21 @@ namespace numsim::materials { /// Isotropic elastic stiffness as a material of its own. /// -/// `linear_elasticity` owns its tangent and computes it once in the -/// constructor, which is optimal while K and G are fixed but leaves no way to -/// recompute it — and no way to order that recomputation correctly if there -/// were one. A material reading its OWN property creates no edge in the -/// property graph, so `elastic::stress` sorts before `elastic::tangent` and -/// would consume a stale stiffness. Splitting the stiffness into a separate -/// material turns that invisible dependency into a real Global edge, which the -/// topological sort then honours. +/// linear_elasticity owns its tangent, and a material reading its own property +/// creates no graph edge — so elastic::stress sorts before elastic::tangent and +/// would consume a stale stiffness. Consuming ANOTHER material's property is a +/// real edge, which the topological sort honours. Also makes the stiffness +/// pluggable: anything producing "tangent" is a drop-in. /// -/// It also makes the stiffness pluggable: anything producing a "tangent" -/// property — anisotropic, temperature-dependent, damage-degraded — is a drop-in -/// replacement, and consumers need not know which. +/// "recompute" binds the update callback. False (the default, matching +/// linear_elasticity) leaves the property with no callback, so the engine skips +/// it and per-call cost is zero — right when constants are fixed (Abaqus PROPS). +/// True is for constants that vary per call (CalculiX interpolates them by +/// temperature). Read once, at construction: that is when the callback is bound. /// -/// ### recompute -/// -/// The update callback is bound ONLY when the "recompute" parameter is true. -/// With it false the tangent is built once, at construction, and the property -/// carries no callback at all — so `ctx.update()` skips it outright and the -/// per-call cost is exactly zero, matching `linear_elasticity`. -/// -/// recompute = false constants are fixed (Abaqus: PROPS cannot vary for a -/// given material name) -/// recompute = true constants change between calls (CalculiX interpolates -/// *USER MATERIAL constants by temperature) -/// -/// The flag is read ONCE, at construction, because that is when the callback is -/// bound — there is no way to attach one later. It is therefore stored by value -/// rather than as a reference into the parameter store: writing "recompute" -/// afterwards changes nothing, and an accessor reading the live parameter would -/// report that the tangent tracks K and G when in fact no callback exists. -/// `recomputes()` reports what was actually bound. -/// -/// Setting it false while the constants do move yields a stale tangent, with -/// the stress still correct — the classic silently-wrong tangent, costing -/// convergence rate rather than accuracy. Nothing detects that automatically: -/// the material cannot see a write it has no callback to observe. Choosing the -/// flag correctly is the caller's responsibility, and -/// RecomputeFalseIgnoresALaterParameterWrite pins what happens if they do not. +/// recompute=false while constants move leaves a stale tangent with the stress +/// still correct — costs convergence rate, not accuracy, and nothing detects it. +/// See RecomputeFalseIgnoresALaterParameterWrite. template class isotropic_tangent final : public material_base, Traits> { @@ -62,11 +39,10 @@ class isotropic_tangent final : base(std::forward(args)...), m_K(base::template get_parameter("K")), m_G(base::template get_parameter("G")), - // By value, not by reference: see the class comment. + // By value: the callback is bound once, so a live reference would let + // recomputes() claim tracking that does not exist. m_recompute(base::template get_parameter("recompute")), - // Bind the callback only when asked. add_output ignores a null one, so - // with recompute=false the property has no callback and the engine - // never visits it. + // add_output ignores a null callback. m_C(base::template add_output( "tangent", base::template get_parameter("recompute") @@ -81,22 +57,18 @@ class isotropic_tangent final input_parameter_controller para{base::parameters()}; para.template insert("K").template add(); para.template insert("G").template add(); - // Defaults to the cheap behaviour, which matches linear_elasticity and is - // right whenever the host cannot vary the constants. para.template insert("recompute").template add(false); return para; } void update_tangent() { - // IIdev comes from plasticity_utils rather than being spelled out again; - // the isotropic basis was already written three times in this repo. const auto I{tmech::eye()}; const auto IIvol{tmech::otimes(I, I) / Dim}; m_C = 3 * m_K * IIvol + 2 * m_G * plasticity_detail::make_IIdev(); } - /// True when this material will follow a change to K or G. + /// Whether a callback was bound, i.e. whether K/G changes are followed. [[nodiscard]] bool recomputes() const noexcept { return m_recompute; } private: diff --git a/include/numsim-materials/materials/linear_stress.h b/include/numsim-materials/materials/linear_stress.h index 92b81da..0538d35 100644 --- a/include/numsim-materials/materials/linear_stress.h +++ b/include/numsim-materials/materials/linear_stress.h @@ -8,17 +8,12 @@ namespace numsim::materials { /// sigma = C : eps, with C supplied by another material. /// -/// The decomposed counterpart to `linear_elasticity`: same physics, but the -/// stiffness arrives as a Global input rather than being owned and computed -/// in the constructor. That single change is what makes the ordering correct — -/// consuming another material's property creates a graph edge, so the engine -/// guarantees the stiffness is produced before this material reads it. A -/// material reading its own property has no such edge (see isotropic_tangent). +/// linear_elasticity with the tangent as a Global input instead of an owned +/// output. That is what makes the ordering correct: consuming another +/// material's property creates a graph edge; reading your own does not. /// -/// Pair it with any tangent generator. `linear_elasticity` remains the better -/// choice whenever the moduli are fixed, since it computes the tangent once and -/// costs one less material in the graph; reach for this pair when the stiffness -/// has to follow something that changes. +/// Pair with any tangent generator. linear_elasticity stays the better choice +/// when the moduli are fixed — one fewer material, tangent computed once. template class linear_stress final : public material_base, Traits> { diff --git a/include/numsim-materials/materials/weighted_sum.h b/include/numsim-materials/materials/weighted_sum.h index d97730d..bd0154d 100644 --- a/include/numsim-materials/materials/weighted_sum.h +++ b/include/numsim-materials/materials/weighted_sum.h @@ -47,11 +47,7 @@ class weighted_sum final // Dynamically create inputs for each term std::size_t i = 0; for (const auto& [weight_name, mat_name] : m_terms_param) { - // A term's tangent normally comes from the same material as its stress. - // An entry in "tangent_sources" overrides that for one term, so a - // constituent whose stiffness lives in its own material (isotropic_tangent - // + linear_stress) can participate. A short or empty list leaves every - // unlisted term wired exactly as before. + // Unlisted terms take the tangent from the material producing the stress. const std::string& tangent_owner = (i < m_tangent_sources.size()) ? m_tangent_sources[i] : mat_name; auto& w = base::template add_input(weight_name, m_weight_property, EdgeKind::Global); @@ -64,8 +60,7 @@ class weighted_sum final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - // Optional: declared with no check. Absent means every term takes its - // tangent from the same material as its stress. + // Optional (no check): absent means every term uses its stress material. para.template insert>("tangent_sources"); para.template insert("terms").template add(); para.template insert("weight_property") diff --git a/include/numsim-materials/umat/material_point_evaluator.h b/include/numsim-materials/umat/material_point_evaluator.h index 45805e5..d2c15f1 100644 --- a/include/numsim-materials/umat/material_point_evaluator.h +++ b/include/numsim-materials/umat/material_point_evaluator.h @@ -46,14 +46,9 @@ class material_point_evaluator { std::string strain_source; /// Material producing the stress the host wants back. std::string stress_source; - /// Material producing the tangent. Unset means "same as stress_source", - /// which is how every monolithic material (linear_elasticity, j2_plasticity) - /// is configured. Set it when the stiffness lives in its own material, as - /// with isotropic_tangent + linear_stress. - /// - /// std::optional rather than an empty string, so "not configured" and - /// "configured as empty" are distinct — the latter is a mistake and should - /// fail to resolve rather than silently fall back. + /// Unset means "same as stress_source", which is how monolithic materials + /// are configured. Set it when the stiffness is its own material. + /// optional, not an empty string, so unset and "" stay distinct. std::optional tangent_source{}; std::string stress_property{"stress"}; std::string tangent_property{"tangent"}; diff --git a/tests/test_tangent_generator.cpp b/tests/test_tangent_generator.cpp index e92f61b..c08f391 100644 --- a/tests/test_tangent_generator.cpp +++ b/tests/test_tangent_generator.cpp @@ -62,19 +62,14 @@ tensor2 uniaxial(T v) { // The ordering this decomposition exists to fix // --------------------------------------------------------------------------- -/// The whole point. Within one material, `stress` sorts BEFORE `tangent` -/// because a material reading its own property creates no graph edge — so -/// linear_elasticity could never recompute its tangent safely. Consuming -/// another material's property creates a real Global edge, and the topological -/// sort then puts the producer first. +/// The whole point: a cross-material tangent is ordered before its consumer, +/// where an intra-material one is not. TEST(TangentGenerator, StiffnessIsOrderedBeforeTheStressThatConsumesIt) { ctx_type ctx; build_decomposed(ctx, K, G, /*recompute=*/true); - // Not a sentinel of 0: the tangent has no inputs, so index 0 is exactly where - // it legitimately lands, and a not-found sentinel of 0 would be - // indistinguishable from the correct answer. The test would then pass on a - // graph that no longer contains the property being ordered. + // optional, not 0: the tangent legitimately lands at index 0, so a 0 sentinel + // would pass even with the property missing from the graph. std::optional i_tangent, i_stress; std::size_t n = 0; for (const auto* prop : ctx.property_execution_order()) { @@ -93,8 +88,7 @@ TEST(TangentGenerator, StiffnessIsOrderedBeforeTheStressThatConsumesIt) { // Equivalence with the monolithic material // --------------------------------------------------------------------------- -/// Decomposed and monolithic must agree exactly — same physics, different -/// graph shape. +/// Same physics, different graph shape: must agree exactly. TEST(TangentGenerator, MatchesLinearElasticityExactly) { ctx_type dec; auto& dec_src = build_decomposed(dec, K, G, /*recompute=*/false); @@ -142,8 +136,7 @@ TEST(TangentGenerator, MatchesLinearElasticityExactly) { // recompute: bound vs unbound callback // --------------------------------------------------------------------------- -/// With recompute=false the property carries NO callback, so the engine skips -/// it entirely — the per-call cost is zero rather than merely small. +/// recompute=false must leave no callback at all, not merely a cheap one. TEST(TangentGenerator, RecomputeFalseLeavesThePropertyWithoutACallback) { ctx_type ctx; build_decomposed(ctx, K, G, /*recompute=*/false); @@ -163,8 +156,7 @@ TEST(TangentGenerator, RecomputeTrueBindsTheCallback) { EXPECT_TRUE(static_cast(prop->traits().update)); } -/// Even with no callback the tangent must be valid: it is built once in the -/// constructor, before any update() runs. +/// Built in the constructor, so valid before the first update(). TEST(TangentGenerator, TangentIsValidBeforeTheFirstUpdate) { ctx_type ctx; build_decomposed(ctx, K, G, /*recompute=*/false); @@ -174,13 +166,7 @@ TEST(TangentGenerator, TangentIsValidBeforeTheFirstUpdate) { EXPECT_NEAR(C(0, 0, 0, 0), K + 4.0 * G / 3.0, 1e-9); } -/// The behaviour the whole design is for: writing new constants in place and -/// having the graph pick them up, with correct ordering, on the next update. -/// -/// Writing goes through the non-const get(), which mutates the value inside -/// the std::any in place. insert() would replace the whole any and, for a type -/// past its small-buffer, relocate the object — invalidating the reference the -/// material bound at construction. +/// What the design is for: write constants, graph picks them up next update. TEST(TangentGenerator, RecomputeTrueFollowsAParameterWrittenInPlace) { ctx_type ctx; auto& src = build_decomposed(ctx, K, G, /*recompute=*/true); @@ -206,9 +192,7 @@ TEST(TangentGenerator, RecomputeTrueFollowsAParameterWrittenInPlace) { EXPECT_NEAR(after, 2 * before, 1e-12); } -/// The counterpart: with recompute=false the write is visible in the parameter -/// but the tangent does NOT follow it. This is the documented trade, and the -/// test exists so the behaviour is pinned rather than discovered. +/// The documented trade, pinned: recompute=false ignores a later write. TEST(TangentGenerator, RecomputeFalseIgnoresALaterParameterWrite) { ctx_type ctx; auto& src = build_decomposed(ctx, K, G, /*recompute=*/false); @@ -234,11 +218,8 @@ TEST(TangentGenerator, RecomputeFalseIgnoresALaterParameterWrite) { // The flag is construction-time only // --------------------------------------------------------------------------- -/// The callback can only be bound at construction, so "recompute" is read once -/// and stored BY VALUE. Reading the live parameter instead would let -/// recomputes() claim the tangent tracks K and G after someone flipped the flag -/// with set_parameter — while no callback exists and the tangent is permanently -/// stale. The accessor must describe what was actually bound. +/// recomputes() must report what was bound, not the live parameter — flipping +/// the flag afterwards cannot bind a callback. TEST(TangentGenerator, RecomputeIsReadOnceAndTheAccessorCannotLie) { ctx_type ctx; build_decomposed(ctx, K, G, /*recompute=*/false); @@ -248,13 +229,10 @@ TEST(TangentGenerator, RecomputeIsReadOnceAndTheAccessorCannotLie) { ASSERT_NE(typed, nullptr); ASSERT_FALSE(typed->recomputes()); - // Flipping the parameter afterwards cannot bind a callback... typed->template set_parameter("recompute", true); - // ... so the accessor must still report false, EXPECT_FALSE(typed->recomputes()) << "recomputes() must report what was bound, not the live parameter"; - // ... and the property must still carry no callback. const auto* prop = ctx.find_property("stiffness", "tangent"); ASSERT_NE(prop, nullptr); EXPECT_FALSE(static_cast(prop->traits().update)); @@ -264,10 +242,9 @@ TEST(TangentGenerator, RecomputeIsReadOnceAndTheAccessorCannotLie) { // Composition with a pre-existing consumer // --------------------------------------------------------------------------- -/// The header claims a tangent producer is a drop-in for any consumer. That has -/// to be tested against a consumer that already existed, not only against two -/// contexts built here. small_strain_plasticity takes its tangent from a NAMED -/// source, so it needs no change: point elastic_source at the generator. +/// The drop-in claim, tested against a PRE-EXISTING consumer. +/// small_strain_plasticity already names its tangent source, so it needs no +/// change: point elastic_source at the generator. TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { auto drive = [](bool decomposed) { ctx_type ctx; @@ -336,11 +313,8 @@ TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { // The optional tangent_source // --------------------------------------------------------------------------- -/// "Optional" in this framework means declared with no check: check_parameter -/// only runs registered checks, and the JSON visitor skips absent keys. These -/// pin that an ABSENT tangent_source falls back to the stress source, while a -/// SUPPLIED one is honoured — the two must be distinguishable, which an -/// empty-string sentinel could not express. +/// Absent tangent_source falls back to the stress source; supplied is honoured. +/// The two must stay distinguishable. TEST(TangentSource, AbsentFallsBackToTheStressSource) { ctx_type ctx; param_type p; @@ -354,7 +328,6 @@ TEST(TangentSource, AbsentFallsBackToTheStressSource) { ctx.create>(p); ctx.finalize(); - // No tangent_source configured: the monolithic material supplies both. u::material_point_evaluator::config cfg; cfg.strain_source = "strain_in"; cfg.stress_source = "elastic"; @@ -370,16 +343,14 @@ TEST(TangentSource, SuppliedResolvesTheTangentElsewhere) { u::material_point_evaluator::config cfg; cfg.strain_source = "strain_in"; - cfg.stress_source = "elastic"; // linear_stress: publishes only "stress" - // Without this the evaluator cannot find a tangent at all. + cfg.stress_source = "elastic"; // linear_stress publishes only "stress" EXPECT_THROW(u::material_point_evaluator(ctx, cfg), u::fatal_error); cfg.tangent_source = "stiffness"; // the generator EXPECT_NO_THROW(u::material_point_evaluator(ctx, cfg)); } -/// The decomposed pair driving a UMAT end to end — the configuration the review -/// showed was impossible before tangent_source existed. +/// The decomposed pair driving a UMAT end to end. TEST(TangentSource, DecomposedPairDrivesTheEvaluatorLikeLinearElasticity) { ctx_type dec; build_decomposed(dec, K, G, /*recompute=*/false); From 2b2d38dbec10361c674cdccb90509eca2d1cac5e Mon Sep 17 00:00:00 2001 From: petlenz Date: Fri, 14 Aug 2026 21:55:00 +0200 Subject: [PATCH 05/11] materials: moduli become graph inputs; remove the recompute flag MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The recompute flag encoded a choice — does this stiffness follow its moduli — as a boolean, and a boolean can disagree with how the graph was actually built. It did: the flag was a live reference into the parameter store while the callback was bound once at construction, so flipping it made recomputes() claim tracking that did not exist. That was a symptom. The cause is that a parameter has no edge to the material reading it, so nothing orders a change to it against the values derived from it. K and G are now Global inputs. Which material you wire IS the choice, and it cannot desynchronise from anything: constant_scalar fixed moduli, a plain (non-history) property any "value" producer moduli that vary, e.g. with temperature constant_scalar publishes a PLAIN property deliberately. statev_map enumerates history properties, so a history-valued constant would consume one STATEV slot per integration point, per constant, and would have to be named in the exclusion list forever. Measured: with a plain property, nstatv is 0 with the moduli not mentioned at all. An earlier prototype used external_scalar_source and did cost those slots — that material is history-based because it exists for time, where the old/new pair is the point. The callback is now always bound, because inputs are not wired until finalize() and nothing can be computed in a constructor that reads them. It self-guards on the moduli, so fixed constants cost two comparisons rather than a rank-4 rebuild. Measured per ctx.update(): 28.1 ns guarded against 308.5 ns unguarded, and 31.7 ns for the old parameter-with-no-callback path — so the wired version is not slower than what it replaces, and 11x faster than rebuilding blindly. FixedModuliAreRebuiltExactlyOnce pins the guard via a recomputations() counter. One behaviour genuinely changes: the tangent is built on the first update rather than at construction, since the inputs do not exist before finalize(). Every consumer goes through ctx.update() first, but TangentIsBuiltOnTheFirstUpdateNotAtConstruction records it. EveryProducerIsOrderedBeforeItsConsumer now checks the whole chain — K and G before the stiffness, the stiffness before the stress — rather than one pair. Note this branch no longer uses set_parameter at all: the test writes through ctx.get_mutable on the constant's published property instead. #26 could be rebased off #25 and reviewed independently. 194 tests pass. --- include/numsim-materials/default_materials.h | 2 + .../materials/constant_scalar.h | 48 ++++ .../materials/isotropic_tangent.h | 77 ++++--- tests/test_tangent_generator.cpp | 206 +++++++++--------- 4 files changed, 191 insertions(+), 142 deletions(-) create mode 100644 include/numsim-materials/materials/constant_scalar.h diff --git a/include/numsim-materials/default_materials.h b/include/numsim-materials/default_materials.h index dbc56e9..1ba4672 100644 --- a/include/numsim-materials/default_materials.h +++ b/include/numsim-materials/default_materials.h @@ -6,6 +6,7 @@ #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/solvers/vector_newton.h" #include "numsim-materials/materials/scalar_stepper.h" +#include "numsim-materials/materials/constant_scalar.h" #include "numsim-materials/materials/isotropic_tangent.h" #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_stress.h" @@ -71,6 +72,7 @@ void register_default_materials() { auto& factory = material_factory::instance(); factory.template register_type>("scalar_stepper"); factory.template register_type>("linear_elasticity"); + factory.template register_type>("constant_scalar"); factory.template register_type>("isotropic_tangent"); factory.template register_type>("linear_stress"); factory.template register_type>("autocatalytic_reaction"); diff --git a/include/numsim-materials/materials/constant_scalar.h b/include/numsim-materials/materials/constant_scalar.h new file mode 100644 index 0000000..bfe384a --- /dev/null +++ b/include/numsim-materials/materials/constant_scalar.h @@ -0,0 +1,48 @@ +#ifndef NUMSIM_MATERIALS_CONSTANT_SCALAR_H +#define NUMSIM_MATERIALS_CONSTANT_SCALAR_H + +#include "numsim-materials/core/material_base.h" + +namespace numsim::materials { + +/// A fixed scalar, published as a graph property. +/// +/// Lets a quantity that some material derives from — a modulus, a yield stress — +/// be a real graph edge instead of a parameter. The consumer is then ordered +/// after it by the topological sort, which is not true of a parameter: nothing +/// connects a parameter to the material that reads it. +/// +/// The property is PLAIN, not history, so statev_map never sees it and it costs +/// no STATEV slot. Use external_scalar_source instead when the value genuinely +/// changes per call and a consumer needs its old/new pair. +/// +/// No update callback: the value is set once, at construction, so the engine +/// skips this property entirely. +template +class constant_scalar final + : public material_base, Traits> { +public: + using base = material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit constant_scalar(Args&&... args) + : base(std::forward(args)...), + m_value(base::template add_output("value")) { + m_value = base::template get_parameter("value"); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + para.template insert("value").template add(); + return para; + } + +private: + value_type& m_value; +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_CONSTANT_SCALAR_H diff --git a/include/numsim-materials/materials/isotropic_tangent.h b/include/numsim-materials/materials/isotropic_tangent.h index 82c7965..3d399fd 100644 --- a/include/numsim-materials/materials/isotropic_tangent.h +++ b/include/numsim-materials/materials/isotropic_tangent.h @@ -1,29 +1,30 @@ #ifndef NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H #define NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H +#include #include #include "numsim-materials/core/material_base.h" #include "numsim-materials/materials/plasticity_utils.h" namespace numsim::materials { -/// Isotropic elastic stiffness as a material of its own. +/// Isotropic elastic stiffness as a material, with K and G as graph inputs. /// /// linear_elasticity owns its tangent, and a material reading its own property /// creates no graph edge — so elastic::stress sorts before elastic::tangent and /// would consume a stale stiffness. Consuming ANOTHER material's property is a -/// real edge, which the topological sort honours. Also makes the stiffness -/// pluggable: anything producing "tangent" is a drop-in. +/// real edge, which the topological sort honours. /// -/// "recompute" binds the update callback. False (the default, matching -/// linear_elasticity) leaves the property with no callback, so the engine skips -/// it and per-call cost is zero — right when constants are fixed (Abaqus PROPS). -/// True is for constants that vary per call (CalculiX interpolates them by -/// temperature). Read once, at construction: that is when the callback is bound. +/// The moduli are inputs rather than parameters for the same reason: a +/// parameter has no edge to the material that reads it, so nothing orders a +/// change to it against the values derived from it. Wire from constant_scalar +/// when they are fixed, or from any material publishing "value" when they vary +/// (temperature dependence). Which one you wire IS the choice — there is no +/// flag that can disagree with how the graph was built. /// -/// recompute=false while constants move leaves a stale tangent with the stress -/// still correct — costs convergence rate, not accuracy, and nothing detects it. -/// See RecomputeFalseIgnoresALaterParameterWrite. +/// The callback is always bound: inputs are not wired until finalize(), so +/// nothing can be computed in the constructor. It self-guards on the moduli, so +/// the fixed case costs two comparisons rather than a rank-4 rebuild. template class isotropic_tangent final : public material_base, Traits> { @@ -37,45 +38,51 @@ class isotropic_tangent final template explicit isotropic_tangent(Args&&... args) : base(std::forward(args)...), - m_K(base::template get_parameter("K")), - m_G(base::template get_parameter("G")), - // By value: the callback is bound once, so a live reference would let - // recomputes() claim tracking that does not exist. - m_recompute(base::template get_parameter("recompute")), - // add_output ignores a null callback. m_C(base::template add_output( - "tangent", - base::template get_parameter("recompute") - ? &isotropic_tangent::update_tangent - : nullptr)) { - // Always compute once, so the tangent is valid before the first update() - // whether or not it will ever be recomputed. - update_tangent(); - } + "tangent", &isotropic_tangent::update_tangent)), + m_K(base::template add_input( + base::template get_parameter("K_source"), "value", + EdgeKind::Global)), + m_G(base::template add_input( + base::template get_parameter("G_source"), "value", + EdgeKind::Global)) {} static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("K").template add(); - para.template insert("G").template add(); - para.template insert("recompute").template add(false); + // Both sources must publish a scalar property called "value", the + // convention scalar_identity_weight and constant_scalar already follow. + para.template insert("K_source").template add(); + para.template insert("G_source").template add(); return para; } void update_tangent() { + if (m_valid && m_K.get() == m_K_cached && m_G.get() == m_G_cached) return; + m_K_cached = m_K.get(); + m_G_cached = m_G.get(); + m_valid = true; + const auto I{tmech::eye()}; const auto IIvol{tmech::otimes(I, I) / Dim}; - m_C = 3 * m_K * IIvol + - 2 * m_G * plasticity_detail::make_IIdev(); + m_C = 3 * m_K_cached * IIvol + + 2 * m_G_cached * plasticity_detail::make_IIdev(); + ++m_recomputations; } - /// Whether a callback was bound, i.e. whether K/G changes are followed. - [[nodiscard]] bool recomputes() const noexcept { return m_recompute; } + /// How often the stiffness was actually rebuilt. Diagnostics for the guard: + /// with fixed moduli this stays at 1 however many updates run. + [[nodiscard]] std::size_t recomputations() const noexcept { + return m_recomputations; + } private: - const value_type& m_K; - const value_type& m_G; - const bool m_recompute; tensor4& m_C; + const input_property& m_K; + const input_property& m_G; + value_type m_K_cached{}; + value_type m_G_cached{}; + bool m_valid{false}; + std::size_t m_recomputations{0}; }; } // namespace numsim::materials diff --git a/tests/test_tangent_generator.cpp b/tests/test_tangent_generator.cpp index c08f391..c7556e5 100644 --- a/tests/test_tangent_generator.cpp +++ b/tests/test_tangent_generator.cpp @@ -1,8 +1,10 @@ #include #include #include +#include #include #include "numsim-materials/core/material_context.h" +#include "numsim-materials/materials/constant_scalar.h" #include "numsim-materials/materials/isotropic_tangent.h" #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_isotropic_hardening.h" @@ -11,6 +13,7 @@ #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/umat/external_state_source.h" #include "numsim-materials/umat/material_point_evaluator.h" +#include "numsim-materials/umat/statev_map.h" namespace { @@ -27,18 +30,26 @@ using tensor4 = tmech::tensor; constexpr T K = 166.67; constexpr T G = 76.92; -/// strain source -> isotropic_tangent -> linear_stress -nm::external_strain_source& build_decomposed(ctx_type& ctx, T k, T g, - bool recompute) { +/// strain source + constant moduli -> isotropic_tangent -> linear_stress +nm::external_strain_source& build_decomposed(ctx_type& ctx, T k, T g) { param_type p; p.insert("name", "strain_in"); auto& src = ctx.create>(p); + p.clear(); + p.insert("name", "K"); + p.insert("value", k); + ctx.create>(p); + + p.clear(); + p.insert("name", "G"); + p.insert("value", g); + ctx.create>(p); + p.clear(); p.insert("name", "stiffness"); - p.insert("K", k); - p.insert("G", g); - p.insert("recompute", recompute); + p.insert("K_source", "K"); + p.insert("G_source", "G"); ctx.create>(p); p.clear(); @@ -59,29 +70,34 @@ tensor2 uniaxial(T v) { } // --------------------------------------------------------------------------- -// The ordering this decomposition exists to fix +// Ordering — what the decomposition exists for // --------------------------------------------------------------------------- -/// The whole point: a cross-material tangent is ordered before its consumer, -/// where an intra-material one is not. -TEST(TangentGenerator, StiffnessIsOrderedBeforeTheStressThatConsumesIt) { +/// A cross-material property is ordered before its consumer; an intra-material +/// one is not. Both the moduli and the stiffness are edges here. +TEST(TangentGenerator, EveryProducerIsOrderedBeforeItsConsumer) { ctx_type ctx; - build_decomposed(ctx, K, G, /*recompute=*/true); + build_decomposed(ctx, K, G); - // optional, not 0: the tangent legitimately lands at index 0, so a 0 sentinel + // optional, not 0: a producer legitimately lands at index 0, so a 0 sentinel // would pass even with the property missing from the graph. - std::optional i_tangent, i_stress; + std::optional i_k, i_g, i_tangent, i_stress; std::size_t n = 0; for (const auto* prop : ctx.property_execution_order()) { const auto& id = prop->traits().id; + if (id.owner == "K" && id.name == "value") i_k = n; + if (id.owner == "G" && id.name == "value") i_g = n; if (id.owner == "stiffness" && id.name == "tangent") i_tangent = n; if (id.owner == "elastic" && id.name == "stress") i_stress = n; ++n; } + ASSERT_TRUE(i_k.has_value() && i_g.has_value()); ASSERT_TRUE(i_tangent.has_value()) << "stiffness::tangent is not in the graph"; ASSERT_TRUE(i_stress.has_value()) << "elastic::stress is not in the graph"; - EXPECT_LT(*i_tangent, *i_stress) - << "the stiffness must be produced before the stress reads it"; + + EXPECT_LT(*i_k, *i_tangent); + EXPECT_LT(*i_g, *i_tangent); + EXPECT_LT(*i_tangent, *i_stress); } // --------------------------------------------------------------------------- @@ -91,7 +107,7 @@ TEST(TangentGenerator, StiffnessIsOrderedBeforeTheStressThatConsumesIt) { /// Same physics, different graph shape: must agree exactly. TEST(TangentGenerator, MatchesLinearElasticityExactly) { ctx_type dec; - auto& dec_src = build_decomposed(dec, K, G, /*recompute=*/false); + auto& dec_src = build_decomposed(dec, K, G); ctx_type mono; nm::external_strain_source* mono_src = nullptr; @@ -122,7 +138,6 @@ TEST(TangentGenerator, MatchesLinearElasticityExactly) { EXPECT_DOUBLE_EQ(a(i, j), b(i, j)) << "step " << step; } - // ... and the stiffness itself. const auto& Cd = dec.get("stiffness", "tangent"); const auto& Cm = mono.get("elastic", "tangent"); for (int i = 0; i < 3; ++i) @@ -133,109 +148,81 @@ TEST(TangentGenerator, MatchesLinearElasticityExactly) { } // --------------------------------------------------------------------------- -// recompute: bound vs unbound callback +// Constants as materials // --------------------------------------------------------------------------- -/// recompute=false must leave no callback at all, not merely a cheap one. -TEST(TangentGenerator, RecomputeFalseLeavesThePropertyWithoutACallback) { +/// constant_scalar publishes a PLAIN property, so it costs no STATEV slot and +/// needs no exclusion. A history property here would be one wasted slot per +/// integration point, per constant. +TEST(TangentGenerator, ConstantsCostNoStatevSlot) { ctx_type ctx; - build_decomposed(ctx, K, G, /*recompute=*/false); - - const auto* prop = ctx.find_property("stiffness", "tangent"); - ASSERT_NE(prop, nullptr); - EXPECT_FALSE(static_cast(prop->traits().update)) - << "recompute=false must not bind an update callback"; + build_decomposed(ctx, K, G); + // Only the host-driven strain is excluded; the moduli are not mentioned. + const u::statev_map map(ctx, {{"strain_in", "strain"}}); + EXPECT_EQ(map.nstatv(), 0u); } -TEST(TangentGenerator, RecomputeTrueBindsTheCallback) { +/// Inputs are not wired until finalize(), so nothing can be computed in the +/// constructor: the stiffness is built on the first update, not before it. +/// Everything that reads it goes through ctx.update() first, but the ordering +/// is worth pinning because it differs from the parameter-based version. +TEST(TangentGenerator, TangentIsBuiltOnTheFirstUpdateNotAtConstruction) { ctx_type ctx; - build_decomposed(ctx, K, G, /*recompute=*/true); + auto& src = build_decomposed(ctx, K, G); - const auto* prop = ctx.find_property("stiffness", "tangent"); - ASSERT_NE(prop, nullptr); - EXPECT_TRUE(static_cast(prop->traits().update)); -} - -/// Built in the constructor, so valid before the first update(). -TEST(TangentGenerator, TangentIsValidBeforeTheFirstUpdate) { - ctx_type ctx; - build_decomposed(ctx, K, G, /*recompute=*/false); - - const auto& C = ctx.get("stiffness", "tangent"); - // C_1111 = K + 4G/3 - EXPECT_NEAR(C(0, 0, 0, 0), K + 4.0 * G / 3.0, 1e-9); -} - -/// What the design is for: write constants, graph picks them up next update. -TEST(TangentGenerator, RecomputeTrueFollowsAParameterWrittenInPlace) { - ctx_type ctx; - auto& src = build_decomposed(ctx, K, G, /*recompute=*/true); + auto* typed = + dynamic_cast*>(ctx.find("stiffness")); + ASSERT_NE(typed, nullptr); + EXPECT_EQ(typed->recomputations(), 0u); const auto eps = uniaxial(0.001); src.bind(eps, eps); ctx.update(); - const T before = ctx.get("elastic", "stress")(0, 0); - EXPECT_NEAR(before, (K + 4.0 * G / 3.0) * 0.001, 1e-12); - - // Double the moduli in place, as a props writer would. - auto* stiffness = ctx.find("stiffness"); - ASSERT_NE(stiffness, nullptr); - auto* typed = dynamic_cast*>(stiffness); - ASSERT_NE(typed, nullptr); - ASSERT_TRUE(typed->recomputes()); - typed->template set_parameter("K", 2 * K); - typed->template set_parameter("G", 2 * G); - ctx.update(); - const T after = ctx.get("elastic", "stress")(0, 0); - EXPECT_NEAR(after, 2 * (K + 4.0 * G / 3.0) * 0.001, 1e-12); - EXPECT_NEAR(after, 2 * before, 1e-12); + EXPECT_EQ(typed->recomputations(), 1u); + EXPECT_NEAR(ctx.get("stiffness", "tangent")(0, 0, 0, 0), + K + 4.0 * G / 3.0, 1e-9); } -/// The documented trade, pinned: recompute=false ignores a later write. -TEST(TangentGenerator, RecomputeFalseIgnoresALaterParameterWrite) { +/// The self-guard is what makes fixed moduli free: the callback is always bound, +/// but it rebuilds only when a modulus actually moves. +TEST(TangentGenerator, FixedModuliAreRebuiltExactlyOnce) { ctx_type ctx; - auto& src = build_decomposed(ctx, K, G, /*recompute=*/false); - - const auto eps = uniaxial(0.001); - src.bind(eps, eps); - ctx.update(); - const T before = ctx.get("elastic", "stress")(0, 0); - + auto& src = build_decomposed(ctx, K, G); auto* typed = dynamic_cast*>(ctx.find("stiffness")); ASSERT_NE(typed, nullptr); - EXPECT_FALSE(typed->recomputes()); - typed->template set_parameter("K", 2 * K); - - ctx.update(); - EXPECT_DOUBLE_EQ(ctx.get("elastic", "stress")(0, 0), before) - << "recompute=false must leave the tangent as built"; -} + const auto eps = uniaxial(0.001); + src.bind(eps, eps); + for (int i = 0; i < 50; ++i) ctx.update(); -// --------------------------------------------------------------------------- -// The flag is construction-time only -// --------------------------------------------------------------------------- + EXPECT_EQ(typed->recomputations(), 1u) + << "the guard must skip every update after the first"; +} -/// recomputes() must report what was bound, not the live parameter — flipping -/// the flag afterwards cannot bind a callback. -TEST(TangentGenerator, RecomputeIsReadOnceAndTheAccessorCannotLie) { +/// And when a modulus does move, the stiffness follows on the next update — +/// with the ordering guaranteed by the edge, not by registration order. +TEST(TangentGenerator, StiffnessFollowsAChangedModulus) { ctx_type ctx; - build_decomposed(ctx, K, G, /*recompute=*/false); - + auto& src = build_decomposed(ctx, K, G); auto* typed = dynamic_cast*>(ctx.find("stiffness")); ASSERT_NE(typed, nullptr); - ASSERT_FALSE(typed->recomputes()); - typed->template set_parameter("recompute", true); + const auto eps = uniaxial(0.001); + src.bind(eps, eps); + ctx.update(); + const T before = ctx.get("elastic", "stress")(0, 0); + EXPECT_NEAR(before, (K + 4.0 * G / 3.0) * 0.001, 1e-12); + + // Write the constant material's published value directly. + ctx.get_mutable("K", "value") = 2 * K; + ctx.update(); - EXPECT_FALSE(typed->recomputes()) - << "recomputes() must report what was bound, not the live parameter"; - const auto* prop = ctx.find_property("stiffness", "tangent"); - ASSERT_NE(prop, nullptr); - EXPECT_FALSE(static_cast(prop->traits().update)); + EXPECT_EQ(typed->recomputations(), 2u); + EXPECT_NEAR(ctx.get("elastic", "stress")(0, 0), + (2 * K + 4.0 * G / 3.0) * 0.001, 1e-12); } // --------------------------------------------------------------------------- @@ -252,14 +239,23 @@ TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { p.insert("name", "strain_in"); auto& src = ctx.create>(p); - p.clear(); - p.insert("name", "stiffness"); if (decomposed) { - p.insert("K", K); - p.insert("G", G); - p.insert("recompute", false); + p.clear(); + p.insert("name", "K"); + p.insert("value", K); + ctx.create>(p); + p.clear(); + p.insert("name", "G"); + p.insert("value", G); + ctx.create>(p); + p.clear(); + p.insert("name", "stiffness"); + p.insert("K_source", "K"); + p.insert("G_source", "G"); ctx.create>(p); } else { + p.clear(); + p.insert("name", "stiffness"); p.insert("strain_producer_name", "strain_in"); p.insert("K", K); p.insert("G", G); @@ -287,7 +283,7 @@ TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { std::vector out; for (int step = 1; step <= 30; ++step) { - const auto eps = uniaxial(0.02 * step); // enough to pass sigma_0 = 50 + const auto eps = uniaxial(0.02 * step); src.bind(eps, eps); ctx.update(); out.push_back(ctx.get("j2", "stress")(0, 0)); @@ -308,18 +304,16 @@ TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { EXPECT_TRUE(went_plastic) << "the path must yield for this to mean anything"; } - // --------------------------------------------------------------------------- // The optional tangent_source // --------------------------------------------------------------------------- /// Absent tangent_source falls back to the stress source; supplied is honoured. -/// The two must stay distinguishable. TEST(TangentSource, AbsentFallsBackToTheStressSource) { ctx_type ctx; param_type p; p.insert("name", "strain_in"); - auto& src = ctx.create>(p); + ctx.create>(p); p.clear(); p.insert("name", "elastic"); p.insert("strain_producer_name", "strain_in"); @@ -333,27 +327,25 @@ TEST(TangentSource, AbsentFallsBackToTheStressSource) { cfg.stress_source = "elastic"; ASSERT_FALSE(cfg.tangent_source.has_value()); EXPECT_NO_THROW(u::material_point_evaluator(ctx, cfg)); - - (void)src; } TEST(TangentSource, SuppliedResolvesTheTangentElsewhere) { ctx_type ctx; - build_decomposed(ctx, K, G, /*recompute=*/false); + build_decomposed(ctx, K, G); u::material_point_evaluator::config cfg; cfg.strain_source = "strain_in"; cfg.stress_source = "elastic"; // linear_stress publishes only "stress" EXPECT_THROW(u::material_point_evaluator(ctx, cfg), u::fatal_error); - cfg.tangent_source = "stiffness"; // the generator + cfg.tangent_source = "stiffness"; EXPECT_NO_THROW(u::material_point_evaluator(ctx, cfg)); } /// The decomposed pair driving a UMAT end to end. TEST(TangentSource, DecomposedPairDrivesTheEvaluatorLikeLinearElasticity) { ctx_type dec; - build_decomposed(dec, K, G, /*recompute=*/false); + build_decomposed(dec, K, G); u::material_point_evaluator::config dcfg; dcfg.strain_source = "strain_in"; dcfg.stress_source = "elastic"; From ffc29d7216d214a90c15cacf14923f231729e529 Mon Sep 17 00:00:00 2001 From: petlenz Date: Fri, 14 Aug 2026 23:24:34 +0200 Subject: [PATCH 06/11] materials: fix the tangent_sources index trap; give weighted_sum a test file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Five findings from the review of #26. tangent_sources was matched to terms by index with no length check, so only TRAILING terms could be left unlisted. A user with two constituents writing ["stiff_B"] for the second one silently got it applied to the FIRST: both names resolve, wire_inputs() succeeds, and the only symptom is a wrong summed tangent — degraded or failed global Newton while the stresses still converge correctly, which is close to the hardest failure to trace. A longer list was truncated in silence. It now requires one entry per term, and "" keeps that term's own tangent so a non-leading term can be overridden alone. I removed that escape myself in an earlier commit when switching to contains(), which is what made leading-only overrides the only possibility. weighted_sum had no test file at all — the one material in the repo with zero coverage, and I had just added a parameter to it. It now has five tests: the mixture rule itself, absent/overridden/empty-entry tangent ownership, and the length check. The last two are mutation-verified. isotropic_tangent's doc told users to wire varying moduli from "any material publishing value", but the property name was hardcoded and the only host-driven scalar in the repo publishes "state" — following the comment gave wire_inputs(): property 'temperature::value' not found at finalize(). Rather than correct the comment, K_property/G_property now default to "value" and can be overridden, so the documented path works. isotropic_tangent gains invalidate(). The memo guard keys on K and G, which assumes this material is the sole writer of its "tangent" property; material_context::get_mutable() hands out a mutable reference to exactly that, so anything writing it without writing it back would wedge the memo permanently with no way out. material_point_evaluator::config's tangent_source moved from the middle of the aggregate to the end. std::string converts implicitly to std::optional, so a downstream aggregate initialiser would still compile and silently re-bind its trailing arguments — a confusing upgrade error in a header meant to be embedded in third-party UMATs. isotropic_damage's class doc now mentions the optional tangent_source, which the change was about and the doc still denied. 189 -> 194 tests. --- .../materials/isotropic_damage.h | 2 + .../materials/isotropic_tangent.h | 29 ++- .../numsim-materials/materials/weighted_sum.h | 26 ++- .../umat/material_point_evaluator.h | 12 +- tests/CMakeLists.txt | 1 + tests/test_weighted_sum.cpp | 215 ++++++++++++++++++ 6 files changed, 269 insertions(+), 16 deletions(-) create mode 100644 tests/test_weighted_sum.cpp diff --git a/include/numsim-materials/materials/isotropic_damage.h b/include/numsim-materials/materials/isotropic_damage.h index d551ec5..3660aa7 100644 --- a/include/numsim-materials/materials/isotropic_damage.h +++ b/include/numsim-materials/materials/isotropic_damage.h @@ -24,6 +24,8 @@ namespace numsim::materials { /// damage_source::damage, damage_source::d_damage — from propagation law /// state_source::d_equivalent_strain — from state function /// yield_source::is_yielding — from yield function +/// The tangent may come from a different material than the stress: set the +/// optional "tangent_source". Absent, both come from "elastic_source". template class isotropic_damage final : public material_base, Traits> { diff --git a/include/numsim-materials/materials/isotropic_tangent.h b/include/numsim-materials/materials/isotropic_tangent.h index 3d399fd..587a974 100644 --- a/include/numsim-materials/materials/isotropic_tangent.h +++ b/include/numsim-materials/materials/isotropic_tangent.h @@ -18,9 +18,9 @@ namespace numsim::materials { /// The moduli are inputs rather than parameters for the same reason: a /// parameter has no edge to the material that reads it, so nothing orders a /// change to it against the values derived from it. Wire from constant_scalar -/// when they are fixed, or from any material publishing "value" when they vary -/// (temperature dependence). Which one you wire IS the choice — there is no -/// flag that can disagree with how the graph was built. +/// when they are fixed, or from any scalar producer when they vary (temperature +/// dependence) — set K_property/G_property if it does not publish under "value". +/// Which one you wire IS the choice, and no flag can disagree with it. /// /// The callback is always bound: inputs are not wired until finalize(), so /// nothing can be computed in the constructor. It self-guards on the moduli, so @@ -41,18 +41,25 @@ class isotropic_tangent final m_C(base::template add_output( "tangent", &isotropic_tangent::update_tangent)), m_K(base::template add_input( - base::template get_parameter("K_source"), "value", + base::template get_parameter("K_source"), + base::template get_parameter("K_property"), EdgeKind::Global)), m_G(base::template add_input( - base::template get_parameter("G_source"), "value", + base::template get_parameter("G_source"), + base::template get_parameter("G_property"), EdgeKind::Global)) {} static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - // Both sources must publish a scalar property called "value", the - // convention scalar_identity_weight and constant_scalar already follow. para.template insert("K_source").template add(); para.template insert("G_source").template add(); + // Defaults to the scalar-output convention constant_scalar and + // scalar_identity_weight follow. Override for a source that publishes under + // another name — external_scalar_source publishes "state", for instance. + para.template insert("K_property") + .template add(std::string{"value"}); + para.template insert("G_property") + .template add(std::string{"value"}); return para; } @@ -69,6 +76,14 @@ class isotropic_tangent final ++m_recomputations; } + /// Force a rebuild on the next update. + /// + /// The guard keys on K and G, which assumes this material is the only writer + /// of its "tangent" property. material_context::get_mutable() hands out a + /// mutable reference to exactly that, so anything writing it without writing + /// it back would otherwise wedge the memo permanently. + void invalidate() noexcept { m_valid = false; } + /// How often the stiffness was actually rebuilt. Diagnostics for the guard: /// with fixed moduli this stays at 1 however many updates run. [[nodiscard]] std::size_t recomputations() const noexcept { diff --git a/include/numsim-materials/materials/weighted_sum.h b/include/numsim-materials/materials/weighted_sum.h index bd0154d..3eaeb1e 100644 --- a/include/numsim-materials/materials/weighted_sum.h +++ b/include/numsim-materials/materials/weighted_sum.h @@ -1,6 +1,7 @@ #ifndef WEIGHTED_SUM_H #define WEIGHTED_SUM_H +#include #include #include "numsim-materials/core/material_base.h" @@ -18,6 +19,8 @@ namespace numsim::materials { /// Parameters: /// "name": material name /// "terms": vector> — (weight_name, constituent_name) pairs +/// "tangent_sources": optional vector — one per term, naming a +/// different material for that term's tangent; "" keeps the term's own. template class weighted_sum final : public material_base, Traits> { @@ -39,17 +42,30 @@ class weighted_sum final m_weight_property(base::template get_parameter("weight_property")), m_stress_property(base::template get_parameter("stress_property")), m_tangent_property(base::template get_parameter("tangent_property")), - m_has_tangent_sources(base::m_parameter_handler.contains("tangent_sources")), - m_tangent_sources(m_has_tangent_sources + m_tangent_sources(base::m_parameter_handler.contains("tangent_sources") ? base::template get_parameter>("tangent_sources") : std::vector{}) { + // Positional, so a shorter list would silently shift every override onto + // the wrong constituent — both names resolve and the only symptom is a + // wrong summed tangent. Require one entry per term, and let "" mean "this + // term keeps its own", so a non-leading term can be overridden alone. + if (!m_tangent_sources.empty() && + m_tangent_sources.size() != m_terms_param.size()) + throw std::runtime_error( + "weighted_sum '" + base::name() + "': tangent_sources has " + + std::to_string(m_tangent_sources.size()) + " entries but there are " + + std::to_string(m_terms_param.size()) + + " terms — supply one per term (\"\" keeps a term's own tangent) or " + "omit it entirely"); + // Dynamically create inputs for each term std::size_t i = 0; for (const auto& [weight_name, mat_name] : m_terms_param) { - // Unlisted terms take the tangent from the material producing the stress. + const bool overridden = + i < m_tangent_sources.size() && !m_tangent_sources[i].empty(); const std::string& tangent_owner = - (i < m_tangent_sources.size()) ? m_tangent_sources[i] : mat_name; + overridden ? m_tangent_sources[i] : mat_name; auto& w = base::template add_input(weight_name, m_weight_property, EdgeKind::Global); auto& s = base::template add_input(mat_name, m_stress_property, EdgeKind::Global); auto& c = base::template add_input(tangent_owner, m_tangent_property, EdgeKind::Global); @@ -61,6 +77,7 @@ class weighted_sum final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; // Optional (no check): absent means every term uses its stress material. + // If given, one entry per term; "" keeps that term's own tangent. para.template insert>("tangent_sources"); para.template insert("terms").template add(); para.template insert("weight_property") @@ -104,7 +121,6 @@ class weighted_sum final const std::string& m_weight_property; const std::string& m_stress_property; const std::string& m_tangent_property; - const bool m_has_tangent_sources; const std::vector m_tangent_sources; std::vector m_terms; }; diff --git a/include/numsim-materials/umat/material_point_evaluator.h b/include/numsim-materials/umat/material_point_evaluator.h index d2c15f1..614ea26 100644 --- a/include/numsim-materials/umat/material_point_evaluator.h +++ b/include/numsim-materials/umat/material_point_evaluator.h @@ -46,10 +46,6 @@ class material_point_evaluator { std::string strain_source; /// Material producing the stress the host wants back. std::string stress_source; - /// Unset means "same as stress_source", which is how monolithic materials - /// are configured. Set it when the stiffness is its own material. - /// optional, not an empty string, so unset and "" stay distinct. - std::optional tangent_source{}; std::string stress_property{"stress"}; std::string tangent_property{"tangent"}; /// Optional external_scalar_source carrying time; empty to omit. @@ -64,6 +60,14 @@ class material_point_evaluator { /// Additional host-owned history to keep out of STATEV, beyond the strain /// and time sources (which are excluded automatically). std::vector extra_exclusions{}; + /// Unset means "same as stress_source", which is how monolithic materials + /// are configured. Set it when the stiffness is its own material. + /// optional, not an empty string, so unset and "" stay distinct. + /// + /// Appended rather than inserted: this header is embedded in third-party + /// UMATs, and a new field in the middle would silently re-bind the trailing + /// arguments of an existing aggregate initialiser. + std::optional tangent_source{}; }; /// One host call's arguments. diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d03b294..1f94db2 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -24,6 +24,7 @@ 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_tangent_generator test_tangent_generator.cpp) +add_numsim_test(test_weighted_sum test_weighted_sum.cpp) target_link_libraries(test_umat_interface PRIVATE Threads::Threads) # Data dumper for plotting (not a test — standalone executable) diff --git a/tests/test_weighted_sum.cpp b/tests/test_weighted_sum.cpp new file mode 100644 index 0000000..4562a60 --- /dev/null +++ b/tests/test_weighted_sum.cpp @@ -0,0 +1,215 @@ +#include +#include +#include +#include +#include +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/materials/constant_scalar.h" +#include "numsim-materials/materials/isotropic_tangent.h" +#include "numsim-materials/materials/linear_elasticity.h" +#include "numsim-materials/materials/linear_stress.h" +#include "numsim-materials/materials/scalar_identity_weight.h" +#include "numsim-materials/materials/scalar_stepper.h" +#include "numsim-materials/materials/weighted_sum.h" +#include "numsim-materials/umat/external_state_source.h" + +namespace { + +namespace nm = numsim::materials; + +using policy = nm::material_policy_default; +using T = policy::value_type; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; +using tensor2 = tmech::tensor; +using tensor4 = tmech::tensor; +using terms_type = std::vector>; + +constexpr T KA = 100.0, GA = 40.0; +constexpr T KB = 300.0, GB = 140.0; + +/// A weight in [0,1] published as "value", via scalar_identity_weight over a +/// scalar_stepper's history. +void add_weight(ctx_type& ctx, const std::string& name, T increment) { + param_type p; + p.insert("name", name + "_drv"); + p.insert("increment", increment); + ctx.create>(p); + p.clear(); + p.insert("name", name); + p.insert("source", name + "_drv::state"); + ctx.create>(p); +} + +void add_monolithic(ctx_type& ctx, const std::string& name, T k, T g) { + param_type p; + p.insert("name", name); + p.insert("strain_producer_name", "strain_in"); + p.insert("K", k); + p.insert("G", g); + ctx.create>(p); +} + +/// A constituent whose stiffness lives in its OWN material, so its stress and +/// tangent come from two different names. +void add_decomposed(ctx_type& ctx, const std::string& name, T k, T g) { + param_type p; + p.insert("name", name + "_K"); + p.insert("value", k); + ctx.create>(p); + p.clear(); + p.insert("name", name + "_G"); + p.insert("value", g); + ctx.create>(p); + p.clear(); + p.insert("name", name + "_stiff"); + p.insert("K_source", name + "_K"); + p.insert("G_source", name + "_G"); + ctx.create>(p); + p.clear(); + p.insert("name", name); + p.insert("tangent_source", name + "_stiff"); + p.insert("strain_source", "strain_in"); + ctx.create>(p); +} + +nm::external_strain_source& add_strain(ctx_type& ctx) { + param_type p; + p.insert("name", "strain_in"); + return ctx.create>(p); +} + +tensor2 uniaxial(T v) { + tensor2 e; + e.fill(0.0); + e(0, 0) = v; + return e; +} + +// --------------------------------------------------------------------------- +// Baseline: the mixture rule itself +// --------------------------------------------------------------------------- + +/// Two monolithic constituents at weights 0.25 and 0.5. Both stress and tangent +/// must be the weighted sum. +TEST(WeightedSum, SumsStressAndTangentByWeight) { + ctx_type ctx; + auto& src = add_strain(ctx); + add_weight(ctx, "wA", 0.25); + add_weight(ctx, "wB", 0.5); + add_monolithic(ctx, "matA", KA, GA); + add_monolithic(ctx, "matB", KB, GB); + + param_type p; + p.insert("name", "mix"); + p.insert("terms", {{"wA", "matA"}, {"wB", "matB"}}); + ctx.create>(p); + ctx.finalize(); + + src.bind(uniaxial(0.001), uniaxial(0.001)); + ctx.update(); + + const T cA = KA + 4.0 * GA / 3.0; + const T cB = KB + 4.0 * GB / 3.0; + EXPECT_NEAR(ctx.get("mix", "stress")(0, 0), + (0.25 * cA + 0.5 * cB) * 0.001, 1e-10); + EXPECT_NEAR(ctx.get("mix", "tangent")(0, 0, 0, 0), + 0.25 * cA + 0.5 * cB, 1e-9); +} + +// --------------------------------------------------------------------------- +// tangent_sources +// --------------------------------------------------------------------------- + +/// Absent: every term takes its tangent from the material producing its stress. +TEST(WeightedSum, AbsentTangentSourcesUsesEachTermsOwnMaterial) { + ctx_type ctx; + auto& src = add_strain(ctx); + add_weight(ctx, "wA", 1.0); + add_monolithic(ctx, "matA", KA, GA); + + param_type p; + p.insert("name", "mix"); + p.insert("terms", {{"wA", "matA"}}); + ctx.create>(p); + ctx.finalize(); + + src.bind(uniaxial(0.001), uniaxial(0.001)); + ctx.update(); + EXPECT_NEAR(ctx.get("mix", "tangent")(0, 0, 0, 0), + KA + 4.0 * GA / 3.0, 1e-9); +} + +/// The case tangent_sources exists for: a constituent whose stiffness is its own +/// material, so its stress and tangent have different owners. +TEST(WeightedSum, OverridesOneTermsTangentOwner) { + ctx_type ctx; + auto& src = add_strain(ctx); + add_weight(ctx, "wA", 1.0); + add_decomposed(ctx, "matA", KA, GA); + + param_type p; + p.insert("name", "mix"); + p.insert("terms", {{"wA", "matA"}}); + p.insert>("tangent_sources", {"matA_stiff"}); + ctx.create>(p); + ctx.finalize(); + + src.bind(uniaxial(0.001), uniaxial(0.001)); + ctx.update(); + EXPECT_NEAR(ctx.get("mix", "tangent")(0, 0, 0, 0), + KA + 4.0 * GA / 3.0, 1e-9); +} + +/// An empty entry keeps that term's own tangent, so a NON-LEADING term can be +/// overridden alone. Without it, a positional list could only ever override a +/// prefix, and a short list would silently shift every override left. +TEST(WeightedSum, EmptyEntryKeepsATermsOwnTangent) { + ctx_type ctx; + auto& src = add_strain(ctx); + add_weight(ctx, "wA", 0.5); + add_weight(ctx, "wB", 0.5); + add_monolithic(ctx, "matA", KA, GA); // keeps its own tangent + add_decomposed(ctx, "matB", KB, GB); // tangent lives elsewhere + + param_type p; + p.insert("name", "mix"); + p.insert("terms", {{"wA", "matA"}, {"wB", "matB"}}); + p.insert>("tangent_sources", {"", "matB_stiff"}); + ctx.create>(p); + ctx.finalize(); + + src.bind(uniaxial(0.001), uniaxial(0.001)); + ctx.update(); + + const T cA = KA + 4.0 * GA / 3.0; + const T cB = KB + 4.0 * GB / 3.0; + EXPECT_NEAR(ctx.get("mix", "tangent")(0, 0, 0, 0), + 0.5 * cA + 0.5 * cB, 1e-9); +} + +/// A shorter list is rejected. Positional matching means it would otherwise +/// apply the override to the WRONG constituent: both names resolve, +/// wire_inputs() succeeds, and the only symptom is a wrong summed tangent — +/// degraded Newton convergence while the stresses still converge correctly. +TEST(WeightedSum, RejectsATangentSourcesListThatDoesNotMatchTheTermCount) { + auto build = [](std::vector sources) { + ctx_type ctx; + add_strain(ctx); + add_weight(ctx, "wA", 0.5); + add_weight(ctx, "wB", 0.5); + add_monolithic(ctx, "matA", KA, GA); + add_decomposed(ctx, "matB", KB, GB); + param_type p; + p.insert("name", "mix"); + p.insert("terms", {{"wA", "matA"}, {"wB", "matB"}}); + p.insert>("tangent_sources", std::move(sources)); + ctx.create>(p); + }; + + EXPECT_THROW(build({"matB_stiff"}), std::runtime_error); // too short + EXPECT_THROW(build({"", "matB_stiff", "extra"}), std::runtime_error); // long + EXPECT_NO_THROW(build({"", "matB_stiff"})); // exact +} + +} // namespace From 27b52930aad952b6a64bb4e4e3f53310e3b08533 Mon Sep 17 00:00:00 2001 From: petlenz Date: Fri, 14 Aug 2026 23:49:24 +0200 Subject: [PATCH 07/11] materials: drop recomputations() and invalidate() from isotropic_tangent MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Neither had a caller outside the tests that existed to exercise them, which is the same reason #25 was closed: do not ship public API nothing uses. recomputations() was a counter for testing the memo guard. invalidate() was added to let a caller recover after writing the "tangent" property directly through material_context::get_mutable(), which would otherwise leave the guard satisfied and the stiffness never rebuilt — but nothing in the repo does that, and get_mutable already documents that the caller is responsible for restoring consistency. The three tests that read the counter now check observable behaviour instead: the tangent is zero before the first update and correct after, repeated updates with fixed moduli stay stable, and a changed modulus still propagates. Two consequences worth stating rather than discovering later. The self-guard is now untested. It is a pure optimisation — 28.1 ns per update against 308.5 ns without it — with no behavioural signature, so removing it would be a silent 11x regression on this path, not a test failure. The get_mutable wedge is reopened. Writing stiffness::tangent directly and not writing it back leaves m_valid true with matching cached moduli, so the guard returns early forever. Recovering means rebuilding the context. 194 tests pass. --- .../materials/isotropic_tangent.h | 17 ----------- tests/test_tangent_generator.cpp | 30 +++++++------------ 2 files changed, 10 insertions(+), 37 deletions(-) diff --git a/include/numsim-materials/materials/isotropic_tangent.h b/include/numsim-materials/materials/isotropic_tangent.h index 587a974..8c1aadc 100644 --- a/include/numsim-materials/materials/isotropic_tangent.h +++ b/include/numsim-materials/materials/isotropic_tangent.h @@ -1,7 +1,6 @@ #ifndef NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H #define NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H -#include #include #include "numsim-materials/core/material_base.h" #include "numsim-materials/materials/plasticity_utils.h" @@ -73,21 +72,6 @@ class isotropic_tangent final const auto IIvol{tmech::otimes(I, I) / Dim}; m_C = 3 * m_K_cached * IIvol + 2 * m_G_cached * plasticity_detail::make_IIdev(); - ++m_recomputations; - } - - /// Force a rebuild on the next update. - /// - /// The guard keys on K and G, which assumes this material is the only writer - /// of its "tangent" property. material_context::get_mutable() hands out a - /// mutable reference to exactly that, so anything writing it without writing - /// it back would otherwise wedge the memo permanently. - void invalidate() noexcept { m_valid = false; } - - /// How often the stiffness was actually rebuilt. Diagnostics for the guard: - /// with fixed moduli this stays at 1 however many updates run. - [[nodiscard]] std::size_t recomputations() const noexcept { - return m_recomputations; } private: @@ -97,7 +81,6 @@ class isotropic_tangent final value_type m_K_cached{}; value_type m_G_cached{}; bool m_valid{false}; - std::size_t m_recomputations{0}; }; } // namespace numsim::materials diff --git a/tests/test_tangent_generator.cpp b/tests/test_tangent_generator.cpp index c7556e5..f725a16 100644 --- a/tests/test_tangent_generator.cpp +++ b/tests/test_tangent_generator.cpp @@ -170,35 +170,29 @@ TEST(TangentGenerator, TangentIsBuiltOnTheFirstUpdateNotAtConstruction) { ctx_type ctx; auto& src = build_decomposed(ctx, K, G); - auto* typed = - dynamic_cast*>(ctx.find("stiffness")); - ASSERT_NE(typed, nullptr); - EXPECT_EQ(typed->recomputations(), 0u); + // Default-constructed until something runs the callback. + EXPECT_DOUBLE_EQ(ctx.get("stiffness", "tangent")(0, 0, 0, 0), 0.0); const auto eps = uniaxial(0.001); src.bind(eps, eps); ctx.update(); - EXPECT_EQ(typed->recomputations(), 1u); EXPECT_NEAR(ctx.get("stiffness", "tangent")(0, 0, 0, 0), K + 4.0 * G / 3.0, 1e-9); } -/// The self-guard is what makes fixed moduli free: the callback is always bound, -/// but it rebuilds only when a modulus actually moves. -TEST(TangentGenerator, FixedModuliAreRebuiltExactlyOnce) { +/// Repeated updates with fixed moduli must keep giving the same answer. The +/// self-guard makes that cheap, but it is an optimisation with no observable +/// behaviour — this checks the observable part. +TEST(TangentGenerator, RepeatedUpdatesWithFixedModuliAreStable) { ctx_type ctx; auto& src = build_decomposed(ctx, K, G); - auto* typed = - dynamic_cast*>(ctx.find("stiffness")); - ASSERT_NE(typed, nullptr); - const auto eps = uniaxial(0.001); src.bind(eps, eps); - for (int i = 0; i < 50; ++i) ctx.update(); - - EXPECT_EQ(typed->recomputations(), 1u) - << "the guard must skip every update after the first"; + ctx.update(); + const T first = ctx.get("elastic", "stress")(0, 0); + for (int i = 0; i < 20; ++i) ctx.update(); + EXPECT_DOUBLE_EQ(ctx.get("elastic", "stress")(0, 0), first); } /// And when a modulus does move, the stiffness follows on the next update — @@ -206,9 +200,6 @@ TEST(TangentGenerator, FixedModuliAreRebuiltExactlyOnce) { TEST(TangentGenerator, StiffnessFollowsAChangedModulus) { ctx_type ctx; auto& src = build_decomposed(ctx, K, G); - auto* typed = - dynamic_cast*>(ctx.find("stiffness")); - ASSERT_NE(typed, nullptr); const auto eps = uniaxial(0.001); src.bind(eps, eps); @@ -220,7 +211,6 @@ TEST(TangentGenerator, StiffnessFollowsAChangedModulus) { ctx.get_mutable("K", "value") = 2 * K; ctx.update(); - EXPECT_EQ(typed->recomputations(), 2u); EXPECT_NEAR(ctx.get("elastic", "stress")(0, 0), (2 * K + 4.0 * G / 3.0) * 0.001, 1e-12); } From 5c27bed3166ac882252c9c2bd7f7068632295ed1 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 15 Aug 2026 11:57:06 +0200 Subject: [PATCH 08/11] materials: isotropic_tangent rebuilds unconditionally; drop the memo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit One job: recompute the stiffness from its inputs on every update. No cached moduli, no validity flag, nothing to invalidate. What that removes, beyond the code: the guard was an optimisation with no behavioural signature, so it could only be tested through a counter that existed for the test — and the counter was itself unused API. It also assumed this material was the sole writer of its "tangent" property, so anything writing that through material_context::get_mutable() and not writing it back left the guard permanently satisfied and the stiffness never rebuilt. With no memo there is no wedge, and that review finding closes by construction rather than by adding an escape hatch for it. The cost is real and worth stating plainly. Measured per ctx.update(): isotropic_tangent + linear_stress 308.2 ns linear_elasticity (monolithic) 27.0 ns So this pair is ~11x the monolithic material on the elastic path, and the header now says so and points at linear_elasticity for fixed moduli. That is the right default anyway: the decomposition exists for moduli that CHANGE, where a memo would be missing most of the time, and for pluggable stiffness, where correctness of ordering is the point rather than throughput. No test changed; behaviour is identical, only slower. 194 tests pass. --- .../materials/isotropic_tangent.h | 23 ++++++++----------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/include/numsim-materials/materials/isotropic_tangent.h b/include/numsim-materials/materials/isotropic_tangent.h index 8c1aadc..0222837 100644 --- a/include/numsim-materials/materials/isotropic_tangent.h +++ b/include/numsim-materials/materials/isotropic_tangent.h @@ -21,9 +21,14 @@ namespace numsim::materials { /// dependence) — set K_property/G_property if it does not publish under "value". /// Which one you wire IS the choice, and no flag can disagree with it. /// -/// The callback is always bound: inputs are not wired until finalize(), so -/// nothing can be computed in the constructor. It self-guards on the moduli, so -/// the fixed case costs two comparisons rather than a rank-4 rebuild. +/// One job: rebuild the stiffness from its inputs on every update. No memo, so +/// no cached state to go stale and nothing to invalidate. The callback is always +/// bound, since inputs are not wired until finalize() and nothing can be +/// computed in a constructor that reads them. +/// +/// Costs a rank-4 rebuild per update (measured: ~309 ns against ~28 ns for a +/// memoised version). Where the moduli are fixed and that matters, +/// linear_elasticity computes its tangent once and is the cheaper choice. template class isotropic_tangent final : public material_base, Traits> { @@ -63,24 +68,16 @@ class isotropic_tangent final } void update_tangent() { - if (m_valid && m_K.get() == m_K_cached && m_G.get() == m_G_cached) return; - m_K_cached = m_K.get(); - m_G_cached = m_G.get(); - m_valid = true; - const auto I{tmech::eye()}; const auto IIvol{tmech::otimes(I, I) / Dim}; - m_C = 3 * m_K_cached * IIvol + - 2 * m_G_cached * plasticity_detail::make_IIdev(); + m_C = 3 * m_K.get() * IIvol + + 2 * m_G.get() * plasticity_detail::make_IIdev(); } private: tensor4& m_C; const input_property& m_K; const input_property& m_G; - value_type m_K_cached{}; - value_type m_G_cached{}; - bool m_valid{false}; }; } // namespace numsim::materials From 5b7036d95bb89adf0e4c823279a7fa04af9db7e6 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 15 Aug 2026 12:19:16 +0200 Subject: [PATCH 09/11] tests: compare tensors with tmech instead of component loops MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit test_tangent_generator compared tensors by nested index loops and by picking out single components — C(0,0,0,0) against a hand-derived K + 4G/3. That is the style 41d217f already replaced elsewhere in this repo with tmech expressions. Comparisons now go through tmech: TensorsIdentical(a, b) norm(a - b) == 0, for the exactness claims tmech::almost_equal(a,b,e) for the tolerant ones, as test_vector_newton does and expected values are built as tensors rather than scalars — isotropic(K, G) constructs the reference stiffness independently, so the tests check the whole tensor rather than one component that happens to be easy to derive by hand. Fixtures are tmech expressions too: uniaxial() is v * (e1 (x) e1). The evaluator tests compared the host's flat buffers slot by slot. They now convert both sides back to tensors first, which is both closer to what the host consumes and immune to a slot permutation applied on both sides — which a componentwise check would cancel out. Net: 6 index loops replaced by 12 tensor comparisons, and the J2 comparison carries stress tensors through instead of flattening to one component per step. Verified the tests still bite: perturbing the volumetric coefficient from 3 to 3.0001 fails 5 of them. (My first attempt at that check reported zero failures — the grep pattern missed gtest's "[ FAILED ]" spacing, not a coverage hole.) 194 tests pass. --- tests/test_tangent_generator.cpp | 101 ++++++++++++++++++++----------- 1 file changed, 67 insertions(+), 34 deletions(-) diff --git a/tests/test_tangent_generator.cpp b/tests/test_tangent_generator.cpp index f725a16..a901a2d 100644 --- a/tests/test_tangent_generator.cpp +++ b/tests/test_tangent_generator.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include "numsim-materials/core/material_context.h" @@ -14,6 +15,7 @@ #include "numsim-materials/umat/external_state_source.h" #include "numsim-materials/umat/material_point_evaluator.h" #include "numsim-materials/umat/statev_map.h" +#include "numsim-materials/umat/tensor_conversion.h" namespace { @@ -62,11 +64,33 @@ nm::external_strain_source& build_decomposed(ctx_type& ctx, T k, T g) { return src; } +/// v * (e1 (x) e1) — a uniaxial strain, built as a tensor expression. tensor2 uniaxial(T v) { - tensor2 e; - e.fill(0.0); - e(0, 0) = v; - return e; + tmech::tensor e1; + e1.fill(0.0); + e1(0) = 1.0; + tensor2 out; + out = v * tmech::otimes(e1, e1); + return out; +} + +/// The isotropic stiffness the generator should produce, built independently. +tensor4 isotropic(T k, T g) { + const auto I = tmech::eye(); + const auto IIsym = (tmech::otimesu(I, I) + tmech::otimesl(I, I)) * 0.5; + const auto IIvol = tmech::otimes(I, I) / 3.0; + tensor4 C; + C = 3.0 * k * IIvol + 2.0 * g * (IIsym - IIvol); + return C; +} + +/// Exact tensor equality, expressed through tmech rather than a component loop. +/// norm(a - b) is identically zero only when every component matches. +template +::testing::AssertionResult TensorsIdentical(const A& a, const B& b) { + const auto d = tmech::norm(a - b); + if (d == T{0}) return ::testing::AssertionSuccess(); + return ::testing::AssertionFailure() << "norm(a - b) = " << d; } // --------------------------------------------------------------------------- @@ -131,20 +155,13 @@ TEST(TangentGenerator, MatchesLinearElasticityExactly) { dec.update(); mono.update(); - const auto& a = dec.get("elastic", "stress"); - const auto& b = mono.get("elastic", "stress"); - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - EXPECT_DOUBLE_EQ(a(i, j), b(i, j)) << "step " << step; + EXPECT_TRUE(TensorsIdentical(dec.get("elastic", "stress"), + mono.get("elastic", "stress"))) + << "step " << step; } - const auto& Cd = dec.get("stiffness", "tangent"); - const auto& Cm = mono.get("elastic", "tangent"); - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - for (int k = 0; k < 3; ++k) - for (int l = 0; l < 3; ++l) - EXPECT_DOUBLE_EQ(Cd(i, j, k, l), Cm(i, j, k, l)); + EXPECT_TRUE(TensorsIdentical(dec.get("stiffness", "tangent"), + mono.get("elastic", "tangent"))); } // --------------------------------------------------------------------------- @@ -171,14 +188,16 @@ TEST(TangentGenerator, TangentIsBuiltOnTheFirstUpdateNotAtConstruction) { auto& src = build_decomposed(ctx, K, G); // Default-constructed until something runs the callback. - EXPECT_DOUBLE_EQ(ctx.get("stiffness", "tangent")(0, 0, 0, 0), 0.0); + tensor4 zero; + zero.fill(0.0); + EXPECT_TRUE(TensorsIdentical(ctx.get("stiffness", "tangent"), zero)); const auto eps = uniaxial(0.001); src.bind(eps, eps); ctx.update(); - EXPECT_NEAR(ctx.get("stiffness", "tangent")(0, 0, 0, 0), - K + 4.0 * G / 3.0, 1e-9); + EXPECT_TRUE(tmech::almost_equal(ctx.get("stiffness", "tangent"), + isotropic(K, G), 1e-12)); } /// Repeated updates with fixed moduli must keep giving the same answer. The @@ -190,9 +209,9 @@ TEST(TangentGenerator, RepeatedUpdatesWithFixedModuliAreStable) { const auto eps = uniaxial(0.001); src.bind(eps, eps); ctx.update(); - const T first = ctx.get("elastic", "stress")(0, 0); + const tensor2 first = ctx.get("elastic", "stress"); for (int i = 0; i < 20; ++i) ctx.update(); - EXPECT_DOUBLE_EQ(ctx.get("elastic", "stress")(0, 0), first); + EXPECT_TRUE(TensorsIdentical(ctx.get("elastic", "stress"), first)); } /// And when a modulus does move, the stiffness follows on the next update — @@ -204,15 +223,18 @@ TEST(TangentGenerator, StiffnessFollowsAChangedModulus) { const auto eps = uniaxial(0.001); src.bind(eps, eps); ctx.update(); - const T before = ctx.get("elastic", "stress")(0, 0); - EXPECT_NEAR(before, (K + 4.0 * G / 3.0) * 0.001, 1e-12); + tensor2 expected; + expected = tmech::dcontract(isotropic(K, G), eps); + EXPECT_TRUE(tmech::almost_equal(ctx.get("elastic", "stress"), + expected, 1e-12)); // Write the constant material's published value directly. ctx.get_mutable("K", "value") = 2 * K; ctx.update(); - EXPECT_NEAR(ctx.get("elastic", "stress")(0, 0), - (2 * K + 4.0 * G / 3.0) * 0.001, 1e-12); + expected = tmech::dcontract(isotropic(2 * K, G), eps); + EXPECT_TRUE(tmech::almost_equal(ctx.get("elastic", "stress"), + expected, 1e-12)); } // --------------------------------------------------------------------------- @@ -271,13 +293,13 @@ TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { ctx.create>(p); ctx.finalize(); - std::vector out; + std::vector> out; for (int step = 1; step <= 30; ++step) { const auto eps = uniaxial(0.02 * step); src.bind(eps, eps); ctx.update(); - out.push_back(ctx.get("j2", "stress")(0, 0)); - out.push_back(ctx.get("j2", "equivalent_plastic_strain")); + out.emplace_back(ctx.get("j2", "stress"), + ctx.get("j2", "equivalent_plastic_strain")); ctx.commit(); } return out; @@ -288,8 +310,12 @@ TEST(TangentGenerator, DrivesJ2PlasticityIdenticallyToLinearElasticity) { ASSERT_EQ(with_generator.size(), with_monolith.size()); bool went_plastic = false; for (std::size_t i = 0; i < with_generator.size(); ++i) { - EXPECT_DOUBLE_EQ(with_generator[i], with_monolith[i]) << "sample " << i; - if (i % 2 == 1 && with_generator[i] > 1e-8) went_plastic = true; + EXPECT_TRUE(TensorsIdentical(with_generator[i].first, + with_monolith[i].first)) + << "stress at step " << i; + EXPECT_DOUBLE_EQ(with_generator[i].second, with_monolith[i].second) + << "equivalent plastic strain at step " << i; + if (with_generator[i].second > 1e-8) went_plastic = true; } EXPECT_TRUE(went_plastic) << "the path must yield for this to mean anything"; } @@ -370,10 +396,17 @@ TEST(TangentSource, DecomposedPairDrivesTheEvaluatorLikeLinearElasticity) { .ddsdde = dd, .statev = dsv}); meval.evaluate({.stran = stran, .dstran = dstran, .stress = ms, .ddsdde = md, .statev = msv}); - for (std::size_t i = 0; i < 6; ++i) - EXPECT_DOUBLE_EQ(ds[i], ms[i]) << "step " << step << " stress " << i; - for (std::size_t i = 0; i < 36; ++i) - EXPECT_DOUBLE_EQ(dd[i], md[i]) << "step " << step << " ddsdde " << i; + + // Compare in tensor space rather than slot by slot: a slot permutation on + // both sides would cancel in a componentwise check, and the tensors are + // what the host actually consumes. + EXPECT_TRUE(TensorsIdentical(u::stress_from_buffer(ds), + u::stress_from_buffer(ms))) + << "stress at step " << step; + EXPECT_TRUE(TensorsIdentical(u::tangent_from_buffer(dd), + u::tangent_from_buffer(md))) + << "tangent at step " << step; + for (std::size_t i = 0; i < 6; ++i) stran[i] += dstran[i]; } } From f888e77276a1d617490b5dccc169e89923c66ca0 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 15 Aug 2026 23:35:50 +0200 Subject: [PATCH 10/11] umat: define models in JSON, with the deck's constants bound into the graph MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Models were registered as std::function builders, so a new material meant recompiling the shared library — which defeats the point of the deck driving the model. The rejection of JSON earlier rested on a performance premise I later measured to be wrong: the builder runs once per THREAD, not per integration point, so a parse costs roughly N_threads, not N_gauss_points. make_json_builder turns a document into the builder the registry already takes: { "materials": [ {"type": "external_strain_source", "name": "strain_in"}, {"type": "constant_scalar", "name": "K", "value": 0}, {"type": "constant_scalar", "name": "G", "value": 0}, {"type": "isotropic_tangent", "name": "stiffness", "K_source": "K", "G_source": "G"}, {"type": "linear_stress", "name": "elastic", "tangent_source": "stiffness", "strain_source": "strain_in"} ], "props": ["K.value", "G.value"] } "props" binds *USER MATERIAL constants positionally to "material.parameter" targets, substituted before creation. Pairing that with constant_scalar is what makes it worth doing: a deck constant enters as a graph PROPERTY, so consumers are ordered after it and follow it, rather than being a parameter with no edge to anything. nlohmann/json is now a fetched dependency rather than reached through __has_include against whatever is installed system-wide. A configuration mechanism cannot be optional; previously the JSON tests silently vanished on a machine without it, and this repo only had them because /usr/local happened to carry a copy. The host-driven source materials are now registered with the runtime factory — without that no JSON document could name them. They are registered from the umat layer rather than from register_default_materials(), so the core defaults keep no dependency on the UMAT code. The PROPS-consistency guard now compares VALUES, not just the count. Writing the tests found the hole: the per-thread cache is keyed on CMNAME, so a second call with a same-length but different PROPS array returned the FIRST call's stiffness with no error. That cannot arise from a well-formed deck — one *MATERIAL name carries one constants array — but the failure was silent, and NPROPS doubles is nothing next to an evaluation. My first test was itself unrealistic (one name, two constant sets); it is now split into the real shape, two names, plus a negative test pinning that reusing a name with different constants is fatal. Both the substitution and the value guard are mutation-verified. 194 -> 201 tests. --- CMakeLists.txt | 20 +- include/numsim-materials/umat/json_model.h | 161 +++++++++++++++ .../numsim-materials/umat/umat_interface.h | 31 +-- tests/CMakeLists.txt | 1 + tests/test_json_model.cpp | 190 ++++++++++++++++++ 5 files changed, 389 insertions(+), 14 deletions(-) create mode 100644 include/numsim-materials/umat/json_model.h create mode 100644 tests/test_json_model.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index dfe3cd7..6556295 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,6 +51,24 @@ else() FetchContent_MakeAvailable(tmech) endif() +# --- Dependencies (nlohmann/json for the configuration layer) --- +# Fetched rather than left optional. Model configuration is meant to be JSON +# driven, so a build without it is missing the primary way to define a material, +# not an extra. It was previously reached only via __has_include against whatever +# happened to be installed system-wide, which meant the JSON tests silently +# vanished on a machine without it. +find_package(nlohmann_json 3.11 QUIET) +if(NOT nlohmann_json_FOUND) + FetchContent_Declare( + nlohmann_json + GIT_REPOSITORY https://github.com/nlohmann/json + GIT_TAG v3.11.3 + GIT_SHALLOW TRUE + ) + set(JSON_BuildTests OFF CACHE INTERNAL "") + FetchContent_MakeAvailable(nlohmann_json) +endif() + # --- Dependencies (Eigen for linear algebra) --- find_package(Eigen3 QUIET) if(NOT Eigen3_FOUND) @@ -76,7 +94,7 @@ target_include_directories(${PROJECT_NAME} ) target_compile_features(${PROJECT_NAME} INTERFACE cxx_std_23) -target_link_libraries(${PROJECT_NAME} INTERFACE numsim-core) +target_link_libraries(${PROJECT_NAME} INTERFACE numsim-core nlohmann_json::nlohmann_json) # tmech is header-only — add its include path without linking a target # (linking would pull it into the install export set) diff --git a/include/numsim-materials/umat/json_model.h b/include/numsim-materials/umat/json_model.h new file mode 100644 index 0000000..9b10516 --- /dev/null +++ b/include/numsim-materials/umat/json_model.h @@ -0,0 +1,161 @@ +#ifndef NUMSIM_MATERIALS_UMAT_JSON_MODEL_H +#define NUMSIM_MATERIALS_UMAT_JSON_MODEL_H + +#include +#include +#include +#include +#include + +#include +#include "numsim-materials/default_materials.h" +#include "numsim-materials/io/json_material_factory.h" +#include "numsim-materials/umat/errors.h" +#include "numsim-materials/umat/external_state_source.h" +#include "numsim-materials/umat/umat_interface.h" + +/// Define a UMAT model from JSON rather than from compiled C++. +/// +/// A builder written as a lambda forces a rebuild of the shared library for +/// every new material, which defeats the point of the deck driving the model. +/// This turns a JSON document into the same builder the registry already takes, +/// so a new material means editing a config file. +/// +/// The document is the one io/json_material_factory already understands, plus +/// an optional "props" array binding the deck's *USER MATERIAL constants to +/// named parameters: +/// +/// { +/// "materials": [ +/// {"type": "external_strain_source", "name": "strain_in"}, +/// {"type": "constant_scalar", "name": "K", "value": 0}, +/// {"type": "constant_scalar", "name": "G", "value": 0}, +/// {"type": "isotropic_tangent", "name": "stiffness", +/// "K_source": "K", "G_source": "G"}, +/// {"type": "linear_stress", "name": "elastic", +/// "tangent_source": "stiffness", "strain_source": "strain_in"} +/// ], +/// "props": ["K.value", "G.value"] +/// } +/// +/// PROPS[i] replaces the parameter named by props[i], written "material.param", +/// before that material is created. Values in the document are therefore +/// placeholders for anything listed there. Pairing this with constant_scalar +/// means a deck constant enters as a graph property, so consumers are ordered +/// after it and follow it — see materials/isotropic_tangent.h. +namespace numsim::materials::umat { + +/// Register the host-driven source materials with the runtime factory. +/// +/// Kept here rather than in register_default_materials() so the core defaults +/// stay free of any dependency on the UMAT layer; these materials only mean +/// something when a host is driving the graph. +template +void register_umat_materials() { + auto& factory = material_factory::instance(); + factory.template register_type>( + "external_strain_source"); + factory.template register_type>( + "external_scalar_source"); +} + +/// Split "material.parameter". Both halves must be non-empty. +inline std::pair split_props_target( + const std::string& target) { + const auto dot = target.find('.'); + if (dot == std::string::npos || dot == 0 || dot + 1 == target.size()) + throw fatal_error( + "json_model: props entry '" + target + + "' must be written \"material.parameter\""); + return {target.substr(0, dot), target.substr(dot + 1)}; +} + +/// Build a registry builder from a JSON document. +/// +/// Parsing happens once, here; the returned builder only substitutes PROPS and +/// creates. Any error in the document surfaces on the first UMAT call for the +/// material, as a fatal_error — a malformed config is a setup fault, not +/// something a smaller increment fixes. +template +typename umat_registry::builder make_json_builder( + const std::string& document) { + nlohmann::json parsed; + try { + parsed = nlohmann::json::parse(document); + } catch (const std::exception& e) { + throw fatal_error(std::string("json_model: cannot parse the model " + "document: ") + + e.what()); + } + if (!parsed.contains("materials") || !parsed["materials"].is_array()) + throw fatal_error("json_model: the document needs a \"materials\" array"); + + // Validate the props bindings now rather than on first use, so a typo is + // reported when the model is registered rather than mid-analysis. + std::vector> bindings; + if (parsed.contains("props")) { + if (!parsed["props"].is_array()) + throw fatal_error("json_model: \"props\" must be an array of " + "\"material.parameter\" strings"); + for (const auto& entry : parsed["props"]) { + if (!entry.is_string()) + throw fatal_error("json_model: every \"props\" entry must be a string"); + auto binding = split_props_target(entry.get()); + const bool known = std::any_of( + parsed["materials"].begin(), parsed["materials"].end(), + [&](const nlohmann::json& m) { + return m.contains("name") && + m["name"].get() == binding.first; + }); + if (!known) + throw fatal_error("json_model: props entry targets material '" + + binding.first + + "', which the document does not define"); + bindings.push_back(std::move(binding)); + } + } + + return [parsed, bindings](material_context& ctx, + std::span props) { + static std::once_flag once; + std::call_once(once, [] { + register_default_materials(); + register_umat_materials(); + }); + + if (props.size() < bindings.size()) + throw fatal_error( + "json_model: the document binds " + std::to_string(bindings.size()) + + " material constants but the deck supplied " + + std::to_string(props.size()) + + " — check *USER MATERIAL, CONSTANTS="); + + // Substitute into a copy, so the registered document stays a template and + // a second thread building the same model is unaffected. + nlohmann::json doc = parsed; + for (std::size_t i = 0; i < bindings.size(); ++i) + for (auto& material : doc["materials"]) + if (material.contains("name") && + material["name"].get() == bindings[i].first) + material[bindings[i].second] = props[i]; + + for (const auto& material : doc["materials"]) + create_from_json(ctx, material); + ctx.finalize(); + }; +} + +/// Register a model defined by a JSON document. +template +void register_json_model( + std::string cmname, const std::string& document, + typename umat_registry::config cfg, + typename plane_stress_evaluator::options ps_opts = {}) { + umat_registry::instance().register_model( + std::move(cmname), make_json_builder(document), std::move(cfg), + ps_opts); +} + +} // namespace numsim::materials::umat + +#endif // NUMSIM_MATERIALS_UMAT_JSON_MODEL_H diff --git a/include/numsim-materials/umat/umat_interface.h b/include/numsim-materials/umat/umat_interface.h index b10ce03..b908a83 100644 --- a/include/numsim-materials/umat/umat_interface.h +++ b/include/numsim-materials/umat/umat_interface.h @@ -207,10 +207,10 @@ class umat_registry { std::unique_ptr ctx; std::unique_ptr solid; std::unique_ptr ps; - /// How many constants the context was built from. The graph is built once - /// and reused, so a later call arriving with a different count would mean - /// the cached parameters no longer describe this material. - std::size_t nprops{0}; + /// The constants the context was built from. The graph is built once and + /// reused, so a later call arriving with different ones would mean the + /// cached parameters no longer describe this material. + std::vector props; }; static std::unordered_mapsecond.nprops != props.size()) + // have distinct names — so anything different means the deck contradicts + // the cached graph, and the constants baked into it would be silently + // wrong for every subsequent call. Comparing the VALUES, not just the + // count: a same-length array with different numbers is the case that + // actually reaches a material, and NPROPS doubles is nothing next to an + // evaluation. + if (!std::equal(it->second.props.begin(), it->second.props.end(), + props.begin(), props.end())) throw fatal_error( "numsim UMAT: material '" + std::string(key) + - "' was built from " + std::to_string(it->second.nprops) + - " constants but this call supplies " + - std::to_string(props.size()) + - " — PROPS must be constant for a given material name"); + "' was built from a different set of " + + std::to_string(it->second.props.size()) + + " constants than this call supplies — PROPS must be constant for a " + "given material name; use distinct *MATERIAL names for distinct " + "constants"); return it->second; } @@ -269,7 +274,7 @@ class umat_registry { thread_state ts; ts.ctx = std::make_unique(); m.build(*ts.ctx, props); - ts.nprops = props.size(); + ts.props.assign(props.begin(), props.end()); if (!ts.ctx->is_finalized()) throw fatal_error( "the builder returned without calling finalize() on the context"); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1f94db2..d65d8c3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -25,6 +25,7 @@ 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_tangent_generator test_tangent_generator.cpp) add_numsim_test(test_weighted_sum test_weighted_sum.cpp) +add_numsim_test(test_json_model test_json_model.cpp) target_link_libraries(test_umat_interface PRIVATE Threads::Threads) # Data dumper for plotting (not a test — standalone executable) diff --git a/tests/test_json_model.cpp b/tests/test_json_model.cpp new file mode 100644 index 0000000..9f6a0b0 --- /dev/null +++ b/tests/test_json_model.cpp @@ -0,0 +1,190 @@ +#include +#include +#include +#include +#include +#include "numsim-materials/umat/json_model.h" + +// The Fortran-callable symbol, so the JSON path is exercised through the real +// ABI rather than only through the C++ evaluator. +NUMSIM_MATERIALS_DEFINE_UMAT(numsim::materials::material_policy_default) + +namespace { + +namespace nm = numsim::materials; +namespace u = numsim::materials::umat; + +using policy = nm::material_policy_default; +using T = policy::value_type; +using registry = u::umat_registry; + +/// Elastic model with the moduli bound to the deck's constants. Nothing here is +/// compiled: adding a material means editing this string. +const char* kElastic = R"({ + "materials": [ + {"type": "external_strain_source", "name": "strain_in"}, + {"type": "constant_scalar", "name": "K", "value": 0}, + {"type": "constant_scalar", "name": "G", "value": 0}, + {"type": "isotropic_tangent", "name": "stiffness", + "K_source": "K", "G_source": "G"}, + {"type": "linear_stress", "name": "elastic", + "tangent_source": "stiffness", "strain_source": "strain_in"} + ], + "props": ["K.value", "G.value"] +})"; + +registry::config elastic_config() { + registry::config cfg; + cfg.strain_source = "strain_in"; + cfg.stress_source = "elastic"; + cfg.tangent_source = "stiffness"; + return cfg; +} + +struct fortran_name { + char buf[80]; + explicit fortran_name(const std::string& s) { + for (auto& c : buf) c = ' '; + for (std::size_t i = 0; i < s.size() && i < 80; ++i) buf[i] = s[i]; + } +}; + +/// DDSDDE(1,1) for a uniaxial increment, through the real umat_ entry point. +T uniaxial_tangent(const std::string& name, const T* props, int nprops) { + const fortran_name cm(name); + T statev[1] = {0}; + const T stran[6] = {0, 0, 0, 0, 0, 0}; + const T dstran[6] = {0.001, 0, 0, 0, 0, 0}; + T stress[6] = {0}, ddsdde[36] = {0}, pnewdt = 1.0; + T sse = 0, spd = 0, scd = 0, rpl = 0, ddsddt[6] = {0}, drplde[6] = {0}; + T drpldt = 0; + const T time[2] = {0, 0}; + T dtime = 0.1; + const T temp = 0, dtemp = 0, predef = 0, dpred = 0, celent = 1; + const T coords[3] = {0}, drot[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + const T dfg[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + int noel = 1, npt = 1, layer = 1, kspt = 1, jstep = 1, kinc = 1; + int ndi = 3, nshr = 3, ntens = 6, nstatv = 0; + + 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, props, &nprops, coords, drot, &pnewdt, + &celent, dfg, dfg, &noel, &npt, &layer, &kspt, &jstep, &kinc, 80); + + EXPECT_DOUBLE_EQ(pnewdt, 1.0); + return ddsdde[0]; +} + +struct Registration { + Registration() { + // One registered document, two deck materials — which is how a real deck + // expresses two parameter sets: distinct *MATERIAL names. + u::register_json_model("JSONSOFT", kElastic, elastic_config()); + u::register_json_model("JSONSTIFF", kElastic, elastic_config()); + u::register_json_model("JSONELASTIC", kElastic, elastic_config()); + } +}; +const Registration registration_{}; + +// --------------------------------------------------------------------------- + +/// The point of the whole exercise: the model is a document, the constants come +/// from the deck, and neither requires recompiling the UMAT. +TEST(JsonModel, DeckConstantsDriveAModelDefinedEntirelyInJson) { + const T soft[2] = {100.0, 40.0}; + const T stiff[2] = {300.0, 140.0}; + + EXPECT_NEAR(uniaxial_tangent("JSONSOFT", soft, 2), + 100.0 + 4.0 * 40.0 / 3.0, 1e-9); + EXPECT_NEAR(uniaxial_tangent("JSONSTIFF", stiff, 2), + 300.0 + 4.0 * 140.0 / 3.0, 1e-9); +} + +/// The values written in the document are placeholders for anything listed in +/// "props" — the deck wins. +TEST(JsonModel, DocumentValuesArePlaceholdersForBoundConstants) { + // The document says 0 for both; if substitution failed the tangent would be + // zero rather than wrong-but-plausible. + const T props[2] = {250.0, 90.0}; + EXPECT_NEAR(uniaxial_tangent("JSONELASTIC", props, 2), + 250.0 + 4.0 * 90.0 / 3.0, 1e-9); +} + +// --------------------------------------------------------------------------- +// Validation, at registration rather than mid-analysis +// --------------------------------------------------------------------------- + +TEST(JsonModel, RejectsAMalformedDocument) { + EXPECT_THROW(u::make_json_builder("{not json"), u::fatal_error); + EXPECT_THROW(u::make_json_builder(R"({"nope": 1})"), u::fatal_error); +} + +TEST(JsonModel, RejectsAPropsEntryThatIsNotMaterialDotParameter) { + const char* doc = R"({ + "materials": [{"type": "constant_scalar", "name": "K", "value": 0}], + "props": ["Kvalue"] + })"; + EXPECT_THROW(u::make_json_builder(doc), u::fatal_error); +} + +/// A props entry naming a material the document does not define is a typo that +/// would otherwise substitute nothing and leave the placeholder in place — a +/// wrong-but-plausible modulus rather than an error. +TEST(JsonModel, RejectsAPropsEntryTargetingAnUndefinedMaterial) { + const char* doc = R"({ + "materials": [{"type": "constant_scalar", "name": "K", "value": 0}], + "props": ["Gee.value"] + })"; + EXPECT_THROW(u::make_json_builder(doc), u::fatal_error); +} + +/// One material NAME carries one PROPS array. Supplying different constants for +/// a name whose graph is already built would otherwise return the first call's +/// stiffness forever, silently — the cached context is keyed on the name. +TEST(JsonModel, ChangingConstantsForOneMaterialNameIsFatal) { + const T first[2] = {100.0, 40.0}; + const T second[2] = {300.0, 140.0}; + + EXPECT_NEAR(uniaxial_tangent("JSONELASTIC", first, 2), + 100.0 + 4.0 * 40.0 / 3.0, 1e-9); + + int fatal_count = 0; + static int* counter = &fatal_count; + u::set_fatal_handler([](const char*) { ++*counter; }); + uniaxial_tangent("JSONELASTIC", second, 2); + u::set_fatal_handler(nullptr); + + EXPECT_EQ(fatal_count, 1) + << "same name, different constants must be reported, not ignored"; +} + +/// Fewer constants than the document binds is a *DEPVAR-style setup error, and +/// must be fatal rather than a cutback. +TEST(JsonModel, TooFewDeckConstantsIsFatal) { + const T only_one[1] = {100.0}; + const fortran_name cm("JSONELASTIC"); + T statev[1] = {0}; + const T stran[6] = {0}, dstran[6] = {0.001, 0, 0, 0, 0, 0}; + T stress[6] = {0}, ddsdde[36] = {0}, pnewdt = 1.0; + T sse = 0, spd = 0, scd = 0, rpl = 0, ddsddt[6] = {0}, drplde[6] = {0}; + T drpldt = 0; + const T time[2] = {0, 0}; + T dtime = 0.1; + const T temp = 0, dtemp = 0, predef = 0, dpred = 0, celent = 1; + const T coords[3] = {0}, drot[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + const T dfg[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + int noel = 1, npt = 1, layer = 1, kspt = 1, jstep = 1, kinc = 1; + int ndi = 3, nshr = 3, ntens = 6, nstatv = 0, nprops = 1; + + bool fatal = false; + u::set_fatal_handler([](const char*) {}); + 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, only_one, &nprops, coords, drot, &pnewdt, + &celent, dfg, dfg, &noel, &npt, &layer, &kspt, &jstep, &kinc, 80); + fatal = (pnewdt == 1.0); // fatal path leaves PNEWDT alone + u::set_fatal_handler(nullptr); + EXPECT_TRUE(fatal) << "a wrong CONSTANTS= count must not request a cutback"; +} + +} // namespace From 8948786f1ca635c8bcc1bba4b1b58c573d5cce15 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 16 Aug 2026 02:26:04 +0200 Subject: [PATCH 11/11] umat: name the deck-constants binding "constants", using the library's :: syntax MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two naming problems in the JSON model layer, both of which invited a reader to guess wrong. "props" collided with the framework's central concept. In this library a PROPERTY is a graph node; naming the deck's numbers after it means two unrelated things share a word in the one document where both appear. It is now "constants", which is also what the deck calls them (*USER MATERIAL, CONSTANTS=). "K.value" invented a second syntax for a material-qualified name. The library already has one — "time::state", "stepper::strain", parsed by connection_source::parse — so the targets are now "K::value" and go through that same parser instead of a private splitter. Note the right-hand side is a PARAMETER here, not a property, which the doc now says outright since the syntax does not distinguish them. Renaming exposed a silent failure the old spelling had all along: an unrecognised top-level key was ignored, so a document still saying "props" was accepted with every constant UNBOUND — the placeholders written in the file became the material's moduli, wrong but plausible and with no diagnostic. The document now rejects any top-level key it does not recognise, and names "constants" specifically when it sees "props". json_to_parameters already warns about unknown keys per material; this is the same check one level up. { "materials": [ ... ], "constants": ["K::value", "G::value"] } 201 -> 202 tests. --- include/numsim-materials/umat/json_model.h | 86 ++++++++++++++-------- tests/test_json_model.cpp | 33 +++++++-- 2 files changed, 82 insertions(+), 37 deletions(-) diff --git a/include/numsim-materials/umat/json_model.h b/include/numsim-materials/umat/json_model.h index 9b10516..8d1c578 100644 --- a/include/numsim-materials/umat/json_model.h +++ b/include/numsim-materials/umat/json_model.h @@ -10,6 +10,7 @@ #include #include "numsim-materials/default_materials.h" #include "numsim-materials/io/json_material_factory.h" +#include "numsim-materials/core/input_types.h" #include "numsim-materials/umat/errors.h" #include "numsim-materials/umat/external_state_source.h" #include "numsim-materials/umat/umat_interface.h" @@ -22,8 +23,11 @@ /// so a new material means editing a config file. /// /// The document is the one io/json_material_factory already understands, plus -/// an optional "props" array binding the deck's *USER MATERIAL constants to -/// named parameters: +/// an optional "constants" array binding the deck's *USER MATERIAL constants to +/// named parameters. It is spelled "constants" rather than "props" because in +/// this library a PROPERTY is a graph node — reusing that word for the deck's +/// numbers would name two unrelated things the same. "constants" is also what +/// the deck itself calls them (*USER MATERIAL, CONSTANTS=). /// /// { /// "materials": [ @@ -35,14 +39,19 @@ /// {"type": "linear_stress", "name": "elastic", /// "tangent_source": "stiffness", "strain_source": "strain_in"} /// ], -/// "props": ["K.value", "G.value"] +/// "constants": ["K::value", "G::value"] /// } /// -/// PROPS[i] replaces the parameter named by props[i], written "material.param", -/// before that material is created. Values in the document are therefore -/// placeholders for anything listed there. Pairing this with constant_scalar -/// means a deck constant enters as a graph property, so consumers are ordered -/// after it and follow it — see materials/isotropic_tangent.h. +/// PROPS[i] replaces the parameter named by constants[i], written +/// "material::parameter" — the same qualified-name syntax the rest of the +/// library uses for wiring ("time::state"), parsed by the same +/// connection_source::parse. Note the right-hand side is a PARAMETER here, not +/// a property. +/// +/// Values written in the document are placeholders for anything listed there. +/// Pairing this with constant_scalar means a deck constant enters as a graph +/// property, so consumers are ordered after it and follow it — see +/// materials/isotropic_tangent.h. namespace numsim::materials::umat { /// Register the host-driven source materials with the runtime factory. @@ -59,15 +68,18 @@ void register_umat_materials() { "external_scalar_source"); } -/// Split "material.parameter". Both halves must be non-empty. -inline std::pair split_props_target( - const std::string& target) { - const auto dot = target.find('.'); - if (dot == std::string::npos || dot == 0 || dot + 1 == target.size()) +/// Parse a "material::parameter" target with the library's existing splitter, +/// so this does not invent a second syntax for a qualified name. +inline connection_source parse_constant_target(const std::string& target) { + try { + auto src = connection_source::parse(target); + if (src.material.empty() || src.property.empty()) throw std::invalid_argument(""); + return src; + } catch (const std::invalid_argument&) { throw fatal_error( - "json_model: props entry '" + target + - "' must be written \"material.parameter\""); - return {target.substr(0, dot), target.substr(dot + 1)}; + "json_model: constants entry '" + target + + "' must be written \"material::parameter\""); + } } /// Build a registry builder from a JSON document. @@ -90,26 +102,42 @@ typename umat_registry::builder make_json_builder( if (!parsed.contains("materials") || !parsed["materials"].is_array()) throw fatal_error("json_model: the document needs a \"materials\" array"); + // An unrecognised top-level key is a setup fault, not something to ignore. + // A document still spelling the binding array "props" would otherwise be + // accepted with every constant silently unbound, leaving the placeholders in + // the document as the material's moduli. json_to_parameters already warns + // about unknown keys per material; this is the same check one level up. + for (const auto& [key, value] : parsed.items()) { + if (key == "materials" || key == "constants") continue; + throw fatal_error( + "json_model: unrecognised top-level key \"" + key + + "\"; the document takes \"materials\" and \"constants\"" + + (key == "props" ? " (the binding array is named \"constants\", since " + "\"property\" already means a graph node here)" + : "")); + } + // Validate the props bindings now rather than on first use, so a typo is // reported when the model is registered rather than mid-analysis. - std::vector> bindings; - if (parsed.contains("props")) { - if (!parsed["props"].is_array()) - throw fatal_error("json_model: \"props\" must be an array of " - "\"material.parameter\" strings"); - for (const auto& entry : parsed["props"]) { + std::vector bindings; + if (parsed.contains("constants")) { + if (!parsed["constants"].is_array()) + throw fatal_error("json_model: \"constants\" must be an array of " + "\"material::parameter\" strings"); + for (const auto& entry : parsed["constants"]) { if (!entry.is_string()) - throw fatal_error("json_model: every \"props\" entry must be a string"); - auto binding = split_props_target(entry.get()); + throw fatal_error( + "json_model: every \"constants\" entry must be a string"); + auto binding = parse_constant_target(entry.get()); const bool known = std::any_of( parsed["materials"].begin(), parsed["materials"].end(), [&](const nlohmann::json& m) { return m.contains("name") && - m["name"].get() == binding.first; + m["name"].get() == binding.material; }); if (!known) - throw fatal_error("json_model: props entry targets material '" + - binding.first + + throw fatal_error("json_model: constants entry targets material '" + + binding.material + "', which the document does not define"); bindings.push_back(std::move(binding)); } @@ -136,8 +164,8 @@ typename umat_registry::builder make_json_builder( for (std::size_t i = 0; i < bindings.size(); ++i) for (auto& material : doc["materials"]) if (material.contains("name") && - material["name"].get() == bindings[i].first) - material[bindings[i].second] = props[i]; + material["name"].get() == bindings[i].material) + material[bindings[i].property] = props[i]; for (const auto& material : doc["materials"]) create_from_json(ctx, material); diff --git a/tests/test_json_model.cpp b/tests/test_json_model.cpp index 9f6a0b0..2af52bb 100644 --- a/tests/test_json_model.cpp +++ b/tests/test_json_model.cpp @@ -30,7 +30,7 @@ const char* kElastic = R"({ {"type": "linear_stress", "name": "elastic", "tangent_source": "stiffness", "strain_source": "strain_in"} ], - "props": ["K.value", "G.value"] + "constants": ["K::value", "G::value"] })"; registry::config elastic_config() { @@ -100,8 +100,8 @@ TEST(JsonModel, DeckConstantsDriveAModelDefinedEntirelyInJson) { 300.0 + 4.0 * 140.0 / 3.0, 1e-9); } -/// The values written in the document are placeholders for anything listed in -/// "props" — the deck wins. +/// Values in the document are placeholders for anything listed in "constants" +/// — the deck wins. TEST(JsonModel, DocumentValuesArePlaceholdersForBoundConstants) { // The document says 0 for both; if substitution failed the tangent would be // zero rather than wrong-but-plausible. @@ -119,21 +119,38 @@ TEST(JsonModel, RejectsAMalformedDocument) { EXPECT_THROW(u::make_json_builder(R"({"nope": 1})"), u::fatal_error); } -TEST(JsonModel, RejectsAPropsEntryThatIsNotMaterialDotParameter) { +/// A document using the old "props" spelling would otherwise be accepted with +/// every constant unbound, leaving the placeholders as the material's moduli — +/// wrong but plausible, and completely silent. +TEST(JsonModel, RejectsAnUnrecognisedTopLevelKey) { + const char* old_spelling = R"({ + "materials": [{"type": "constant_scalar", "name": "K", "value": 0}], + "props": ["K::value"] + })"; + EXPECT_THROW(u::make_json_builder(old_spelling), u::fatal_error); + + const char* typo = R"({ + "materials": [{"type": "constant_scalar", "name": "K", "value": 0}], + "constant": ["K::value"] + })"; + EXPECT_THROW(u::make_json_builder(typo), u::fatal_error); +} + +TEST(JsonModel, RejectsAConstantsEntryThatIsNotQualified) { const char* doc = R"({ "materials": [{"type": "constant_scalar", "name": "K", "value": 0}], - "props": ["Kvalue"] + "constants": ["Kvalue"] })"; EXPECT_THROW(u::make_json_builder(doc), u::fatal_error); } -/// A props entry naming a material the document does not define is a typo that +/// A constants entry naming a material the document does not define is a typo that /// would otherwise substitute nothing and leave the placeholder in place — a /// wrong-but-plausible modulus rather than an error. -TEST(JsonModel, RejectsAPropsEntryTargetingAnUndefinedMaterial) { +TEST(JsonModel, RejectsAConstantsEntryTargetingAnUndefinedMaterial) { const char* doc = R"({ "materials": [{"type": "constant_scalar", "name": "K", "value": 0}], - "props": ["Gee.value"] + "constants": ["Gee::value"] })"; EXPECT_THROW(u::make_json_builder(doc), u::fatal_error); }