diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 203690a..399589e 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -21,3 +21,6 @@ target_link_libraries(stvenant_kirchhoff PRIVATE numsim::codegen) add_executable(first_piola_svk first_piola_svk.cpp) target_link_libraries(first_piola_svk PRIVATE numsim::codegen) + +add_executable(svk_from_energy svk_from_energy.cpp) +target_link_libraries(svk_from_energy PRIVATE numsim::codegen) diff --git a/examples/svk_from_energy.cpp b/examples/svk_from_energy.cpp new file mode 100644 index 0000000..782e844 --- /dev/null +++ b/examples/svk_from_energy.cpp @@ -0,0 +1,12 @@ +// Emits the energy-derived St. Venant–Kirchhoff material (#108 front-end) via the +// StandaloneCxx target. See svk_from_energy_recipe.h. +#include "svk_from_energy_recipe.h" +#include +#include +int main() { + using namespace numsim::codegen; + auto const model = examples::make_svk_from_energy(); + StandaloneCxxTarget target; + for (auto const &f : target.emit(model)) std::cout << f.contents; + return 0; +} diff --git a/examples/svk_from_energy_recipe.h b/examples/svk_from_energy_recipe.h new file mode 100644 index 0000000..96e2603 --- /dev/null +++ b/examples/svk_from_energy_recipe.h @@ -0,0 +1,113 @@ +#ifndef NUMSIM_CODEGEN_EXAMPLES_SVK_FROM_ENERGY_RECIPE_H +#define NUMSIM_CODEGEN_EXAMPLES_SVK_FROM_ENERGY_RECIPE_H + +// St. Venant–Kirchhoff DERIVED FROM ITS ENERGY POTENTIAL — the #108 +// material-compiler front-end. The human writes only the strain-energy density +// +// ψ(E) = ½λ (tr E)² + μ (E : E) +// +// and `add_hyperelastic_potential` derives the 2nd-Piola–Kirchhoff stress and +// the consistent tangent by differentiation: +// +// S = ∂ψ/∂E = λ tr(E) I + 2μ E (NOT written by hand) +// dS/dE = ∂²ψ/∂E² (NOT written by hand) +// +// This is the same material as the hand-written St. Venant–Kirchhoff example — +// the e2e checks the derived S against that closed form — but here the stress +// and tangent fall out of `cas::diff` of ψ, the AceGen "differentiate the +// potential" paradigm. +// +// STRESS/STRAIN MEASURE: (S, Green–Lagrange E), 2nd Piola–Kirchhoff conjugate +// pair — same as the hand-written SVK example. +// CONSTRAINTS: `E` must be a declared INPUT leaf (add_hyperelastic_potential +// recovers its name from the handle and rejects a derived measure), and it is +// roles::Strain (symmetric), which is what makes the derived tangent +// minor-symmetric — a non-symmetric leaf would not. +// VERIFICATION BOUNDARY: the e2e compiles the StandaloneCxx form and checks the +// derived S against the independent closed form + the derived tangent vs FD + +// major symmetry (the Hessian property an energy-derived tangent must have). + +#include + +#include +#include +#include +#include +#include + +namespace numsim::codegen::examples { + +inline ConstitutiveModel make_svk_from_energy() { + using namespace numsim::cas; + + ConstitutiveModel model("SvkFromEnergy"); + + auto lambda = + model.add_parameter("lambda", /*default=*/1.0, "Lame first parameter"); + auto mu = model.add_parameter("mu", /*default=*/0.5, "Shear modulus"); + auto E = model.add_tensor_input("E", /*dim=*/3, /*rank=*/2, roles::Strain); + + // The ONLY physics the human writes: the strain-energy density ψ(E). + auto const trE = trace(E); + auto const psi = 0.5 * lambda * trE * trE + mu * dot(E); // dot(E) = E : E + + // Derive S = ∂ψ/∂E and the consistent tangent dS/dE = ∂²ψ/∂E². + model.add_hyperelastic_potential("S", psi, E, "dS_dE"); + + return model; +} + +// A second, deliberately NON-QUADRATIC potential so the derived tangent is a +// genuine function of E (not the constant SVK tangent) — this exercises the +// second-differentiation machinery on a curvature-bearing Hessian, which a +// linear-stress material cannot. +// +// ψ = ½λ (tr E)² + μ (E:E) + c (E:E)² +// S = ∂ψ/∂E = λ tr(E) I + 2μ E + 4c (E:E) E (cubic in E) +// dS/dE = ∂²ψ/∂E² (non-constant) +inline ConstitutiveModel make_nonlinear_from_energy() { + using namespace numsim::cas; + + ConstitutiveModel model("NonlinearFromEnergy"); + + auto lambda = model.add_parameter("lambda", 1.0, "Lame first parameter"); + auto mu = model.add_parameter("mu", 0.5, "Shear modulus"); + auto c = model.add_parameter("c", 0.4, "Quartic stiffening coefficient"); + auto E = model.add_tensor_input("E", 3, 2, roles::Strain); + + auto const trE = trace(E); + auto const EE = dot(E); // E : E + auto const psi = + 0.5 * lambda * trE * trE + mu * EE + c * EE * EE; + + model.add_hyperelastic_potential("S", psi, E, "dS_dE"); + + return model; +} + +// A potential in a NON-SYMMETRIC leaf F (roles::DeformationGradient), to exercise +// add_hyperelastic_potential's non-symmetric-leaf path and to give the +// minor-symmetry check teeth: ∂²ψ/∂F² is major-symmetric (a Hessian) but NOT +// minor-symmetric, because F carries no i↔j symmetry. +// +// ψ = ½μ (F:F) + c (F:F)² → P = ∂ψ/∂F, dP/dF = ∂²ψ/∂F² (minor-asymmetric) +inline ConstitutiveModel make_nonsymmetric_from_energy() { + using namespace numsim::cas; + + ConstitutiveModel model("NonsymmetricFromEnergy"); + + auto mu = model.add_parameter("mu", 0.5, "Shear modulus"); + auto c = model.add_parameter("c", 0.4, "Quartic stiffening coefficient"); + auto F = model.add_tensor_input("F", 3, 2, roles::DeformationGradient); + + auto const FF = dot(F); // F : F + auto const psi = 0.5 * mu * FF + c * FF * FF; + + model.add_hyperelastic_potential("P", psi, F, "dP_dF"); + + return model; +} + +} // namespace numsim::codegen::examples + +#endif // NUMSIM_CODEGEN_EXAMPLES_SVK_FROM_ENERGY_RECIPE_H diff --git a/include/numsim_codegen/recipe.h b/include/numsim_codegen/recipe.h index 025b5a8..921d606 100644 --- a/include/numsim_codegen/recipe.h +++ b/include/numsim_codegen/recipe.h @@ -25,6 +25,7 @@ #include #include #include // Phase 3b-1: diff(tensor, tensor) +#include // #108: diff(energy ψ, strain) → stress #include #include @@ -881,6 +882,95 @@ class ConstitutiveModel { m_tangents.push_back( {std::move(name), std::move(of_output), std::move(wrt_input)}); } + + // #108 material-compiler front-end: derive a hyperelastic material from its + // energy potential ψ instead of writing the stress by hand. Registers the + // stress output (∂ψ/∂ε) and the consistent tangent (∂²ψ/∂ε² = ∂stress/∂ε), + // both via cas::diff — the AceGen "differentiate the potential" paradigm. + // + // `stress_name` — name of the emitted stress output (roles::Stress). + // `energy` — the scalar strain-energy density ψ (a tensor-to-scalar + // expression built from `strain` + parameters). + // `strain` — the strain input to differentiate against. MUST be a + // registered input leaf (the handle from add_tensor_input), + // not a derived measure, and `energy` must depend on it — + // both are checked, loudly. Its name is recovered from the + // handle so the stress and tangent differentiate w.r.t. the + // SAME leaf (no separate, desyncable name argument). + // `tangent_name` — name of the emitted rank-4 consistent tangent. + // + // The derived tangent inherits its symmetry from `strain`: a symmetric + // (roles::Strain) leaf gives the minor-symmetric tangent; a non-symmetric leaf + // (roles::DeformationGradient) does not. The human writes only ψ; stress and + // tangent are derived. This is the front-end slice of #108 (hyperelastic); the + // plastic return-map front-end and the compile-time reduction pass are separate, + // LATER slices — and note this helper eagerly derives-and-discards ψ (it stores + // only the results), so those slices need a STORED-potential representation, not + // an extension of this eager pattern. + void add_hyperelastic_potential( + std::string stress_name, + cas::expression_holder const &energy, + cas::expression_holder const &strain, + std::string tangent_name) { + // Recover the strain input's registered name by node identity — deriving it + // from the handle (rather than taking a second string) makes the stress diff + // and the tangent diff use the same leaf by construction. + std::string strain_name; + bool is_input = false; + for (auto const &sym : m_tensor_symbols) { + if (sym.second.data().get() == strain.data().get()) { + strain_name = sym.first; + // m_tensor_symbols also holds state-variable current/_old handles, which + // are NOT valid strains; require Category::Input specifically. + for (auto const &in : m_inputs_cache) + if (in.name == strain_name) { + is_input = true; + break; + } + break; + } + } + if (!is_input) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': add_hyperelastic_potential's `strain` must be " + "a registered tensor INPUT (the handle returned by add_tensor_input) — " + "not a state variable or a derived expression.", + m_name)); + } + // The energy must at least mention the strain leaf, else ∂ψ/∂ε ≡ 0 emits an + // inert zero stress and tangent. This is a leaf-PRESENCE check (it does not + // catch a pathological energy where the strain appears but cancels to a zero + // gradient); it rejects the common "built ψ from the wrong input" mistake + // loudly rather than shipping a silently-null material. + LeafCollector lc; + lc.collect_t2s(energy); + bool depends = false; + for (auto const &leaf : lc.tensor_names()) + if (leaf == strain_name) { + depends = true; + break; + } + if (!depends) { + throw std::runtime_error(std::format( + "ConstitutiveModel '{}': the strain-energy passed to " + "add_hyperelastic_potential does not depend on strain '{}'; ∂ψ/∂{} " + "would be identically zero (an inert stress and tangent).", + m_name, strain_name, strain_name)); + } + auto stress = cas::diff(energy, strain); // σ = ∂ψ/∂ε + add_output(stress_name, std::move(stress), roles::Stress); + // add_algorithmic_tangent validates AFTER add_output has committed; roll the + // stress output back on throw so a bad tangent_name leaves the model unchanged + // (strong exception guarantee) instead of half-registered. + try { + add_algorithmic_tangent(std::move(tangent_name), std::move(stress_name), + std::move(strain_name)); // ∂σ/∂ε = ∂²ψ/∂ε² + } catch (...) { + m_outputs.pop_back(); + throw; + } + } + [[nodiscard]] auto algorithmic_tangent_enabled() const noexcept -> bool { return !m_tangents.empty(); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f0aa040..49439de 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -265,6 +265,49 @@ target_link_libraries(svk_check_driver gtest_discover_tests(svk_check_driver) +# ─── #108 material-compiler front-end end-to-end gate ──────────────────────── +# Emit the ENERGY-DERIVED St. Venant–Kirchhoff material (S=∂ψ/∂E, tangent +# ∂²ψ/∂E², both derived by cas::diff of the potential), compile it, and verify +# the derived stress against the independently hand-written closed form plus the +# derived tangent against FD. Proves the material compiler produces the material +# a human would have written by hand. +add_executable(generate_energy_check + generated/generate_energy_check.cpp) +target_include_directories(generate_energy_check + PRIVATE ${CMAKE_SOURCE_DIR}/examples) +target_link_libraries(generate_energy_check + PRIVATE numsim::codegen numsim_codegen_warnings) + +set(GENERATED_SVK_ENERGY_HEADER + ${CMAKE_CURRENT_BINARY_DIR}/generated/SvkEnergyCheck.h) +set(GENERATED_NL_ENERGY_HEADER + ${CMAKE_CURRENT_BINARY_DIR}/generated/NlEnergyCheck.h) +set(GENERATED_NS_ENERGY_HEADER + ${CMAKE_CURRENT_BINARY_DIR}/generated/NsEnergyCheck.h) +add_custom_command( + OUTPUT ${GENERATED_SVK_ENERGY_HEADER} ${GENERATED_NL_ENERGY_HEADER} + ${GENERATED_NS_ENERGY_HEADER} + COMMAND ${CMAKE_COMMAND} -E make_directory + ${CMAKE_CURRENT_BINARY_DIR}/generated + COMMAND $ + ${GENERATED_SVK_ENERGY_HEADER} ${GENERATED_NL_ENERGY_HEADER} + ${GENERATED_NS_ENERGY_HEADER} + DEPENDS generate_energy_check + COMMENT "Generating energy-derived material headers via numsim-codegen" + VERBATIM) + +add_executable(energy_check_driver + generated/energy_check_driver.cpp + ${GENERATED_SVK_ENERGY_HEADER} ${GENERATED_NL_ENERGY_HEADER} + ${GENERATED_NS_ENERGY_HEADER}) +target_include_directories(energy_check_driver + PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/generated) +target_link_libraries(energy_check_driver + PRIVATE numsim::codegen numsim_codegen_warnings_light + GTest::gtest GTest::gtest_main) + +gtest_discover_tests(energy_check_driver) + # ─── numsim-materials end-to-end gate (Phase B) ────────────────────────────── # # The string-asserting NumSimMaterialTargetTest proves the emitted TEXT matches diff --git a/tests/RecipeTest.cpp b/tests/RecipeTest.cpp index 503460d..eb4a9db 100644 --- a/tests/RecipeTest.cpp +++ b/tests/RecipeTest.cpp @@ -260,4 +260,66 @@ TEST(Recipe, ResidualRejectsForeignHandle) { } } +// ── #108 add_hyperelastic_potential guards ────────────────────────────────── + +TEST(Recipe, HyperelasticPotentialDerivesStressAndTangent) { + using namespace numsim::cas; + ConstitutiveModel m("Hyper"); + auto mu = m.add_parameter("mu", 0.5); + auto E = m.add_tensor_input("E", 3, 2, roles::Strain); + m.add_hyperelastic_potential("S", mu * dot(E), E, "dS_dE"); + // Registers exactly one stress output and one tangent. + EXPECT_EQ(m.outputs().size(), 1u); + EXPECT_EQ(m.tangents().size(), 1u); +} + +TEST(Recipe, HyperelasticPotentialRejectsNonInputStrain) { + using namespace numsim::cas; + ConstitutiveModel m("Hyper"); + auto mu = m.add_parameter("mu", 0.5); + auto E = m.add_tensor_input("E", 3, 2, roles::Strain); + // A derived expression, not the registered input leaf, is rejected loudly. + EXPECT_THROW(m.add_hyperelastic_potential("S", mu * dot(E), 2.0 * E, "dS_dE"), + std::runtime_error); +} + +TEST(Recipe, HyperelasticPotentialRejectsStrainIndependentEnergy) { + using namespace numsim::cas; + ConstitutiveModel m("Hyper"); + auto mu = m.add_parameter("mu", 0.5); + auto E = m.add_tensor_input("E", 3, 2, roles::Strain); + auto F = m.add_tensor_input("F", 3, 2, roles::Strain); + // ψ depends on F, not on the strain E we differentiate against → ∂ψ/∂E ≡ 0. + EXPECT_THROW(m.add_hyperelastic_potential("S", mu * dot(F), E, "dS_dE"), + std::runtime_error); +} + +TEST(Recipe, HyperelasticPotentialRejectsStateVariableStrain) { + using namespace numsim::cas; + ConstitutiveModel m("Hyper"); + auto mu = m.add_parameter("mu", 0.5); + auto ep = m.add_tensor_state_variable( + "eps_p", 3, 2, make_expression(std::size_t{3}, std::size_t{2})); + // A state-variable handle is registered in m_tensor_symbols but is NOT a + // kinematic input — it must be rejected, not silently differentiated against. + EXPECT_THROW( + m.add_hyperelastic_potential("S", mu * dot(ep.current), ep.current, "dS_dE"), + std::runtime_error); +} + +TEST(Recipe, HyperelasticPotentialRollsBackStressOutputOnBadTangentName) { + using namespace numsim::cas; + ConstitutiveModel m("Hyper"); + auto mu = m.add_parameter("mu", 0.5); + auto E = m.add_tensor_input("E", 3, 2, roles::Strain); + // An invalid tangent name makes add_algorithmic_tangent throw AFTER the stress + // output was added — the stress output must be rolled back, leaving the model + // clean so a corrected retry succeeds. + EXPECT_THROW(m.add_hyperelastic_potential("S", mu * dot(E), E, "bad name"), + std::runtime_error); + EXPECT_EQ(m.outputs().size(), 0u) << "stress output must be rolled back"; + EXPECT_NO_THROW(m.add_hyperelastic_potential("S", mu * dot(E), E, "dS_dE")); + EXPECT_EQ(m.outputs().size(), 1u); +} + } // namespace numsim::codegen diff --git a/tests/generated/energy_check_driver.cpp b/tests/generated/energy_check_driver.cpp new file mode 100644 index 0000000..bf9cc6c --- /dev/null +++ b/tests/generated/energy_check_driver.cpp @@ -0,0 +1,222 @@ +// #108 e2e driver: compile the ENERGY-DERIVED materials and verify the material +// compiler. Three materials: +// * SvkFromEnergy (quadratic ψ) — the DERIVED stress ∂ψ/∂E matches the +// INDEPENDENTLY hand-written closed form (per #108's "self-FD insufficient" +// caveat); its tangent is constant. +// * NonlinearFromEnergy (quartic ψ) — a genuinely NON-CONSTANT derived tangent, +// which exercises the second-differentiation machinery a linear-stress +// material cannot; checked stress-vs-closed-form AND tangent-vs-FD. +// * NonsymmetricFromEnergy (ψ of a non-symmetric leaf F) — exercises the +// non-symmetric-leaf path; its tangent is major- but NOT minor-symmetric, +// which also proves minor_symmetric() discriminates. +// Symmetry is checked two ways: MAJOR (C_ijkl=C_klij) is automatic for a Hessian +// so it only guards against emit/diff corruption; MINOR (C_ijkl=C_jikl=C_ijlk) is +// the load-bearing property the symmetric roles::Strain leaf confers. + +#include "NlEnergyCheck.h" +#include "NsEnergyCheck.h" +#include "SvkEnergyCheck.h" + +#include + +#include + +#include +#include + +namespace { + +using T2 = tmech::tensor; +using T4 = tmech::tensor; + +constexpr double kLambda = 1.3; +constexpr double kMu = 0.7; +constexpr double kC = 0.4; // quartic coefficient of NonlinearFromEnergy + +// Independent, hand-written St. Venant–Kirchhoff 2nd-PK stress (ψ = ½λ(trE)²+μE:E). +T2 hand_written_S(T2 const &E) { + auto I = tmech::eye(); + return tmech::eval(kLambda * tmech::trace(E) * I + 2.0 * kMu * E); +} + +// Independent, hand-written stress for ψ = ½λ(trE)² + μ(E:E) + c(E:E)²: +// S = λ tr(E) I + 2μ E + 4c (E:E) E. +T2 hand_written_S_nonlinear(T2 const &E) { + auto I = tmech::eye(); + const double EE = tmech::dcontract(E, E); // E : E + return tmech::eval(kLambda * tmech::trace(E) * I + 2.0 * kMu * E + + 4.0 * kC * EE * E); +} + +// True if the rank-4 tangent is MAJOR-symmetric: C_ijkl = C_klij. The pair-swap +// is a tmech::basis_change with the (ij)↔(kl) permutation <3,4,1,2>; it's an +// involution, so equality with the original is convention-independent. +bool major_symmetric(T4 const &C, double tol = 1e-12) { + return tmech::almost_equal( + C, tmech::eval(tmech::basis_change>(C)), tol); +} + +// True if the tangent is MINOR-symmetric: C_ijkl = C_jikl = C_ijlk — the two +// intra-pair swaps <2,1,3,4> and <1,2,4,3>. This is the property the symmetric +// (roles::Strain) leaf confers; unlike major symmetry (automatic for any Hessian) +// it can genuinely fail if the leaf's symmetry space is dropped, so it's the +// load-bearing symmetry check here. +bool minor_symmetric(T4 const &C, double tol = 1e-12) { + return tmech::almost_equal( + C, tmech::eval(tmech::basis_change>(C)), + tol) && + tmech::almost_equal( + C, tmech::eval(tmech::basis_change>(C)), + tol); +} + +// 12 symmetric directions (a redundant, over-determined set — the symmetric +// space is only 6-D). Contracting the tangent against all of them exercises its +// action on symmetric strains; note this is blind to any minor-ASYMMETRIC +// component (checked separately by minor_symmetric()). +std::array directions() { + std::array d; + double b[6][9] = {{1,0,0,0,0,0,0,0,0}, {0,0,0,0,1,0,0,0,0}, {0,0,0,0,0,0,0,0,1}, + {0,1,0,1,0,0,0,0,0}, {0,0,1,0,0,0,1,0,0}, {0,0,0,0,0,1,0,1,0}}; + double m[6][9] = {{1,0.3,0.2,0.3,1,0.1,0.2,0.1,1}, {2,-0.5,0,-0.5,1,0.4,0,0.4,-1}, + {0.1,0.7,-0.2,0.7,0.3,0.5,-0.2,0.5,0.9}, {-1,0.2,0.6,0.2,2,-0.3,0.6,-0.3,0.4}, + {0.5,-0.8,0.1,-0.8,-0.5,0.9,0.1,0.9,1.2}, {1.5,0.4,-0.6,0.4,-1,0.2,-0.6,0.2,0.7}}; + int k = 0; + for (; k < 6; ++k) d[k] = T2{b[k][0],b[k][1],b[k][2],b[k][3],b[k][4],b[k][5],b[k][6],b[k][7],b[k][8]}; + for (int j = 0; j < 6; ++j, ++k) d[k] = T2{m[j][0],m[j][1],m[j][2],m[j][3],m[j][4],m[j][5],m[j][6],m[j][7],m[j][8]}; + return d; +} + +} // namespace + +// The DERIVED stress equals the independently hand-written one. +TEST(MaterialCompilerE2E, DerivedStressMatchesHandWritten) { + for (T2 const &E : {T2{0.03,0.01,0.0, 0.01,-0.02,0.005, 0.0,0.005,0.04}, + T2{0.2,0.1,-0.05, 0.1,0.15,0.08, -0.05,0.08,-0.12}}) { + T2 S; T4 dS; + SvkFromEnergy_compute(kLambda, kMu, E, S, dS); + EXPECT_TRUE(tmech::almost_equal(S, hand_written_S(E), 1e-12)) + << "derived ∂ψ/∂E must equal the hand-written closed form"; + } +} + +// The DERIVED tangent dS/dE = ∂²ψ/∂E² matches central FD of the derived stress, +// across strain states and every symmetric direction. +TEST(MaterialCompilerE2E, DerivedTangentMatchesFD) { + for (T2 const &E : {T2{0.03,0.01,0.0, 0.01,-0.02,0.005, 0.0,0.005,0.04}, + T2{0.2,0.1,-0.05, 0.1,0.15,0.08, -0.05,0.08,-0.12}}) { + T2 S; T4 dS; + SvkFromEnergy_compute(kLambda, kMu, E, S, dS); + const double t = 1e-6; + for (auto const &dE : directions()) { + T2 sp, sm; T4 scratch; + SvkFromEnergy_compute(kLambda, kMu, T2(tmech::eval(E + t * dE)), sp, scratch); + SvkFromEnergy_compute(kLambda, kMu, T2(tmech::eval(E - t * dE)), sm, scratch); + auto fd = tmech::eval((sp - sm) / (2.0 * t)); + auto an = tmech::eval(tmech::dcontract(dS, dE)); + EXPECT_TRUE(tmech::almost_equal(an, fd, 1e-6)) << "derived tangent vs FD"; + } + } +} + +// ── NonlinearFromEnergy: quartic ψ → NON-CONSTANT tangent ──────────────────── + +// The derived stress of the quartic potential matches its hand-written form. +TEST(MaterialCompilerE2E, NonlinearStressMatchesHandWritten) { + for (T2 const &E : {T2{0.03,0.01,0.0, 0.01,-0.02,0.005, 0.0,0.005,0.04}, + T2{0.2,0.1,-0.05, 0.1,0.15,0.08, -0.05,0.08,-0.12}}) { + T2 S; T4 dS; + NonlinearFromEnergy_compute(kLambda, kMu, kC, E, S, dS); + EXPECT_TRUE(tmech::almost_equal(S, hand_written_S_nonlinear(E), 1e-12)); + } +} + +// The NON-CONSTANT derived tangent matches FD — the real test of the +// second-differentiation machinery (SVK's constant tangent cannot exercise it, +// since a linear stress is FD-exact regardless). +TEST(MaterialCompilerE2E, NonlinearTangentMatchesFD) { + for (T2 const &E : {T2{0.05,0.02,0.0, 0.02,-0.03,0.01, 0.0,0.01,0.06}, + T2{0.25,0.12,-0.06, 0.12,0.18,0.09, -0.06,0.09,-0.14}}) { + T2 S; T4 dS; + NonlinearFromEnergy_compute(kLambda, kMu, kC, E, S, dS); + const double t = 1e-6; + for (auto const &dE : directions()) { + T2 sp, sm; T4 scratch; + NonlinearFromEnergy_compute(kLambda, kMu, kC, T2(tmech::eval(E + t * dE)), sp, scratch); + NonlinearFromEnergy_compute(kLambda, kMu, kC, T2(tmech::eval(E - t * dE)), sm, scratch); + auto fd = tmech::eval((sp - sm) / (2.0 * t)); + auto an = tmech::eval(tmech::dcontract(dS, dE)); + EXPECT_TRUE(tmech::almost_equal(an, fd, 1e-6)) << "non-constant tangent vs FD"; + } + } +} + +// Both derived tangents are symmetric. Two distinct checks: +// * MAJOR symmetry (C_ijkl = C_klij) is automatic for any correct Hessian +// ∂²ψ/∂E² (Clairaut), so this mainly guards against a diff-non-commutation +// or an emit/CSE corruption of the rank-4 — NOT a value error (FD covers that). +// * MINOR symmetry (C_ijkl = C_jikl = C_ijlk) is the property the symmetric +// roles::Strain leaf actually confers; it is NOT automatic and would fail if +// the leaf's symmetry space were dropped — this is the load-bearing check. +TEST(MaterialCompilerE2E, DerivedTangentIsSymmetric) { + for (T2 const &E : {T2{0.03,0.01,0.0, 0.01,-0.02,0.005, 0.0,0.005,0.04}, + T2{0.25,0.12,-0.06, 0.12,0.18,0.09, -0.06,0.09,-0.14}}) { + T2 S; T4 dS; + SvkFromEnergy_compute(kLambda, kMu, E, S, dS); + EXPECT_TRUE(major_symmetric(dS)) << "SVK tangent not major-symmetric"; + EXPECT_TRUE(minor_symmetric(dS)) << "SVK tangent not minor-symmetric"; + NonlinearFromEnergy_compute(kLambda, kMu, kC, E, S, dS); + EXPECT_TRUE(major_symmetric(dS)) << "nonlinear tangent not major-symmetric"; + EXPECT_TRUE(minor_symmetric(dS)) << "nonlinear tangent not minor-symmetric"; + } +} + +// ── NonsymmetricFromEnergy: ψ(F) with a NON-SYMMETRIC leaf ─────────────────── + +namespace { +// 12 GENERAL (non-symmetric) directions — F is a deformation gradient (9 DOF). +std::array general_directions() { + std::array d; + double v[12][9] = { + {1,0,0,0,0,0,0,0,0}, {0,1,0,0,0,0,0,0,0}, {0,0,1,0,0,0,0,0,0}, + {0,0,0,1,0,0,0,0,0}, {0,0,0,0,1,0,0,0,0}, {0,0,0,0,0,1,0,0,0}, + {0,0,0,0,0,0,1,0,0}, {0,0,0,0,0,0,0,1,0}, {0,0,0,0,0,0,0,0,1}, + {1,0.4,-0.2,0.3,1,0.1,-0.2,0.1,1}, {0.5,-0.7,0.2,0.9,0.3,-0.4,0.1,0.6,0.8}, + {-0.3,0.8,0.5,-0.6,0.2,0.7,0.4,-0.9,0.1}}; + for (int k = 0; k < 12; ++k) + d[k] = T2{v[k][0],v[k][1],v[k][2],v[k][3],v[k][4],v[k][5],v[k][6],v[k][7],v[k][8]}; + return d; +} +} // namespace + +// The derived P = ∂ψ/∂F for a NON-symmetric leaf matches FD in general (9-DOF) +// directions — proving the non-symmetric-leaf path works, not just the symmetric. +TEST(MaterialCompilerE2E, NonsymmetricTangentMatchesFD) { + constexpr double mu = 0.5, c = 0.4; + for (T2 const &F : {T2{1.1,0.2,0.0, 0.05,0.95,0.1, 0.0,0.08,1.04}, + T2{1.3,0.2,0.1, -0.2,1.1,0.3, 0.05,-0.1,0.9}}) { + T2 P; T4 dP; + NonsymmetricFromEnergy_compute(mu, c, F, P, dP); + const double t = 1e-6; + for (auto const &dF : general_directions()) { + T2 pp, pm; T4 scratch; + NonsymmetricFromEnergy_compute(mu, c, T2(tmech::eval(F + t * dF)), pp, scratch); + NonsymmetricFromEnergy_compute(mu, c, T2(tmech::eval(F - t * dF)), pm, scratch); + auto fd = tmech::eval((pp - pm) / (2.0 * t)); + auto an = tmech::eval(tmech::dcontract(dP, dF)); + EXPECT_TRUE(tmech::almost_equal(an, fd, 1e-6)) << "non-symmetric-leaf tangent vs FD"; + } + } +} + +// ∂²ψ/∂F² is MAJOR-symmetric (Hessian) but NOT minor-symmetric (F is not a +// symmetric leaf). This confirms both that the non-symmetric path is correct AND +// that minor_symmetric() genuinely discriminates (does not pass vacuously). +TEST(MaterialCompilerE2E, NonsymmetricTangentIsMajorButNotMinorSymmetric) { + constexpr double mu = 0.5, c = 0.4; + T2 F{1.3,0.2,0.1, -0.2,1.1,0.3, 0.05,-0.1,0.9}; + T2 P; T4 dP; + NonsymmetricFromEnergy_compute(mu, c, F, P, dP); + EXPECT_TRUE(major_symmetric(dP)) << "Hessian must be major-symmetric"; + EXPECT_FALSE(minor_symmetric(dP)) << "a non-symmetric leaf must NOT give a minor-symmetric tangent"; +} diff --git a/tests/generated/generate_energy_check.cpp b/tests/generated/generate_energy_check.cpp new file mode 100644 index 0000000..7f0a50b --- /dev/null +++ b/tests/generated/generate_energy_check.cpp @@ -0,0 +1,30 @@ +// #108 e2e generator: emit the energy-derived materials to headers the driver +// compiles and verifies. Two materials: SvkFromEnergy (quadratic ψ → constant +// tangent, checked against the hand-written closed form) and NonlinearFromEnergy +// (quartic ψ → non-constant tangent, exercising the second-diff machinery). +#include "svk_from_energy_recipe.h" +#include +#include +#include + +namespace { +bool write_header(numsim::codegen::ConstitutiveModel const &m, char const *path) { + std::ofstream f(path); + for (auto const &e : numsim::codegen::StandaloneCxxTarget{}.emit(m)) + f << e.contents; + return static_cast(f); +} +} // namespace + +int main(int argc, char **argv) { + if (argc < 4) { + std::cerr << "usage: generate_energy_check " + "\n"; + return 2; + } + using namespace numsim::codegen::examples; + if (!write_header(make_svk_from_energy(), argv[1])) return 1; + if (!write_header(make_nonlinear_from_energy(), argv[2])) return 1; + if (!write_header(make_nonsymmetric_from_energy(), argv[3])) return 1; + return 0; +}