diff --git a/src/targets/numsim_material.cpp b/src/targets/numsim_material.cpp index 1d2d4f5..bdb7ae1 100644 --- a/src/targets/numsim_material.cpp +++ b/src/targets/numsim_material.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include @@ -89,13 +90,28 @@ void check_scope(ConstitutiveModel const &model) { // properties with their own update callbacks. Only tensor INTERNAL STATE is // still blocked (Mandel + vector solver) — a tensor stress from a scalar state // needs neither, so it is supported here. - if (!model.tangents().empty()) { - // Tensor stress outputs ARE emitted now, so a tangent's `of_output` resolves - // — but the consistent tangent dσ/dε itself is Phase D (it emits the - // tangent-block properties the solver assembles), so reject it here. - throw std::runtime_error( - "NumSimMaterialTarget: algorithmic tangents are out of scope for the " - "rate material (Phase D emits tangent-block properties)."); + // Algorithmic/consistent tangents dσ/dε ARE emitted (Phase D) as rank-4 tensor + // properties via cas::diff. In THIS target's scope the rate cannot depend on + // strain (rate-leaf guard), so the integrated state is independent of the + // current strain (dx/dε ≡ 0) and dσ/dε = ∂σ/∂ε exactly — no J⁻¹ solve / no + // rank-4 inverse, just diff(stress, strain). Validate each spec resolves. + for (auto const &t : model.tangents()) { + bool of_ok = false; + for (auto const &o : model.outputs()) + if (o.name == t.of_output && o.kind == OutputDecl::Kind::Tensor) of_ok = true; + if (!of_ok) { + throw std::runtime_error( + "NumSimMaterialTarget: tangent '" + t.name + "' differentiates '" + + t.of_output + "', which is not a declared tensor output (stress)."); + } + bool wrt_ok = false; + for (auto const &in : model.inputs()) + if (in.name == t.wrt_input && in.kind == SymbolDecl::Kind::Tensor) wrt_ok = true; + if (!wrt_ok) { + throw std::runtime_error( + "NumSimMaterialTarget: tangent '" + t.name + "' differentiates w.r.t. '" + + t.wrt_input + "', which is not a declared tensor input (strain)."); + } } // Tensor inputs (e.g. strain) ARE wired (Global-edge input_property); scalar // inputs are a separate small follow-up — reject those loudly for now. @@ -303,6 +319,44 @@ auto NumSimMaterialTarget::emit(ConstitutiveModel const &model) const } } + // Consistent tangents dσ/dε: a rank-4 tensor property = cas::diff(stress, strain). + // Folded into `outputs` as rank-4 tensor outputs, so the ctor/callback/member + // emission below handles them uniformly. The strain's roles::Strain symmetric + // space makes diff return the minor-symmetric rank-4 identity (P_sym), so the + // tangent is minor-symmetric as a stress-strain tangent must be. + for (auto const &t : model.tangents()) { + cas::expression_holder sigma; + std::size_t tdim = 0; + for (auto const &o : model.outputs()) { + if (o.name == t.of_output && o.kind == OutputDecl::Kind::Tensor) { + sigma = std::get>(o.expr); + tdim = o.dim; + } + } + cas::expression_holder eps; + for (auto const &[name, h] : model.tensor_symbol_map()) { + if (name == t.wrt_input) eps = h; + } + // check_scope already verified both resolve; guard belt-and-braces so a + // future check_scope/emit drift surfaces as a clear error, not cas::diff UB. + if (!sigma.is_valid() || !eps.is_valid()) { + throw std::runtime_error( + "NumSimMaterialTarget: tangent '" + t.name + + "' could not resolve its stress output / strain input handle."); + } + auto const tangent_expr = cas::diff(sigma, eps); // rank-4 ∂σ/∂ε + + CodeGenContext tc; + CodeEmitPipeline tp(tc); + register_scalars(tc); + for (auto const &[name, expr] : model.tensor_symbol_map()) { + tc.register_symbol_tensor(expr, name); + } + tc.reset(); + auto const trhs = tp.tensor().apply(tangent_expr); + outputs.push_back({t.name, tc.render_statements(" "), trhs, true, tdim, 4}); + } + bool has_tensor = !tensor_inputs.empty(); for (auto const &o : outputs) if (o.is_tensor) has_tensor = true; diff --git a/tests/NumSimMaterialTargetTest.cpp b/tests/NumSimMaterialTargetTest.cpp index c813f25..4740e66 100644 --- a/tests/NumSimMaterialTargetTest.cpp +++ b/tests/NumSimMaterialTargetTest.cpp @@ -253,18 +253,40 @@ TEST(NumSimMaterialTarget, RejectsReservedStateName) { std::string::npos); } -// The stress TENSOR output is now emitted, so a tangent-bearing recipe is no -// longer caught by an output guard — the consistent tangent dσ/dε itself is -// Phase D, so the tangents() guard is now the load-bearing rejection. -TEST(NumSimMaterialTarget, RejectsTangentBearingRecipe) { +// Phase D: the consistent tangent dσ/dε is emitted as a rank-4 tensor property +// via cas::diff(stress, strain). With a strain-INDEPENDENT rate (dα/dε=0) this is +// the exact consistent tangent. For σ=2μ·ε the tangent is 2μ·P_sym, so the +// emitted rank-4 RHS is minor-symmetric (otimesu + otimesl). +TEST(NumSimMaterialTarget, EmitsConsistentTangent) { ConstitutiveModel m("WithTangent"); auto mu = m.add_parameter("mu", 0.5); auto a = m.add_scalar_state_variable("a", make_expression(0.0)); m.add_scalar_evolution_equation(a, mu * a.current); - auto eps = m.add_tensor_input("eps", 3, 2, roles::Strain); + auto eps = m.add_tensor_input("strain", 3, 2, roles::Strain); m.add_output("stress", 2 * mu * eps, roles::Stress); - m.add_algorithmic_tangent("dstress_deps", "stress", "eps"); - EXPECT_NE(emit_throw_message(m).find("algorithmic tangents"), + m.add_algorithmic_tangent("dstress_dstrain", "stress", "strain"); + auto const h = header_of(NumSimMaterialTarget{}.emit(m)); + // rank-4 tensor property + its own callback + EXPECT_NE(h.find("add_output>(\n" + " \"dstress_dstrain\", " + "&WithTangent::update_dstress_dstrain"), + std::string::npos) + << h; + EXPECT_NE(h.find("void update_dstress_dstrain() {"), std::string::npos) << h; + // minor-symmetric: both otimesu AND otimesl present in the tangent RHS + EXPECT_NE(h.find("tmech::otimesu"), std::string::npos) << h; + EXPECT_NE(h.find("tmech::otimesl"), std::string::npos) << h; +} + +// A tangent whose `of_output` is not a declared tensor output is rejected. +TEST(NumSimMaterialTarget, RejectsTangentOfMissingStress) { + ConstitutiveModel m("BadTangent"); + auto mu = m.add_parameter("mu", 0.5); + auto a = m.add_scalar_state_variable("a", make_expression(0.0)); + m.add_scalar_evolution_equation(a, mu * a.current); + m.add_tensor_input("strain", 3, 2, roles::Strain); + m.add_algorithmic_tangent("dstress_dstrain", "stress", "strain"); // no "stress" + EXPECT_NE(emit_throw_message(m).find("not a declared tensor output"), std::string::npos); } diff --git a/tests/generated/generate_numsim_material_check.cpp b/tests/generated/generate_numsim_material_check.cpp index 5c534fd..3825432 100644 --- a/tests/generated/generate_numsim_material_check.cpp +++ b/tests/generated/generate_numsim_material_check.cpp @@ -90,6 +90,9 @@ int main(int argc, char** argv) { viscoelastic.add_scalar_evolution_equation(alpha, K * alpha.current); auto eps = viscoelastic.add_tensor_input("strain", 3, 2, roles::Strain); viscoelastic.add_output("stress", alpha.current * eps); + // Phase D: consistent tangent dσ/dε. Since the rate (K·α) is strain- + // independent, dα/dε=0, so dσ/dε = ∂σ/∂ε = α·P_sym (minor-symmetric). + viscoelastic.add_algorithmic_tangent("dstress_dstrain", "stress", "strain"); } if (!write_header(linear, argv[1])) return 1; diff --git a/tests/generated/numsim_material_check_driver.cpp b/tests/generated/numsim_material_check_driver.cpp index 4463476..b37e435 100644 --- a/tests/generated/numsim_material_check_driver.cpp +++ b/tests/generated/numsim_material_check_driver.cpp @@ -260,6 +260,17 @@ TEST(NumSimMaterialEndToEnd, TensorStressFromScalarStateAndStrain) { const auto sig = read_tensor(ctx, "Viscoelastic", "stress"); EXPECT_NEAR(sig(0, 0), alpha * eps(0, 0), 1e-12); // σ = α·ε on the driven comp EXPECT_NEAR(sig(0, 1), T{0}, 1e-12); // off-diagonal strain is 0 + + // Phase D: consistent tangent dσ/dε = α·P_sym (rank-4, minor-symmetric). + using tensor4 = tmech::tensor; + auto* cp = dynamic_cast< + numsim_core::property*>( + ctx.find_property("Viscoelastic", "dstress_dstrain")); + ASSERT_NE(cp, nullptr); + const auto C = cp->get(); + EXPECT_NEAR(C(0, 0, 0, 0), alpha, 1e-12); // P_sym(0,0,0,0)=1 + EXPECT_NEAR(C(0, 1, 0, 1), alpha * 0.5, 1e-12); // minor-symmetric ½ + EXPECT_NEAR(C(0, 0, 1, 1), T{0}, 1e-12); // off-block zero } #endif // NCG_TENSOR_E2E