#108 MVP: derive a hyperelastic material from its energy potential

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)

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 <numsim_codegen/targets/standalone_cxx.h>
+#include <iostream>
+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;
+}

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 <numsim_codegen/recipe.h>
+
+#include <numsim_cas/scalar/scalar_operators.h>
+#include <numsim_cas/tensor/tensor_definitions.h>
+#include <numsim_cas/tensor/tensor_operators.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_functions.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_operators.h>
+
+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

include/numsim_codegen/recipe.h

@@ -25,6 +25,7 @@
 #include <numsim_cas/scalar/scalar_operators.h>
 #include <numsim_cas/tensor/tensor_definitions.h>
 #include <numsim_cas/tensor/tensor_diff.h> // Phase 3b-1: diff(tensor, tensor)
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_diff.h> // #108: diff(energy ψ, strain) → stress
 
 #include <cstddef>
 #include <format>
@@ -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<cas::tensor_to_scalar_expression> const &energy,
+      cas::expression_holder<cas::tensor_expression> 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();
   }

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 $<TARGET_FILE:generate_energy_check>
+            ${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

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<tensor_zero>(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