materials: isotropic_tangent + linear_stress — a stiffness the graph can order

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)

include/numsim-materials/default_materials.h

@@ -6,7 +6,10 @@
 #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"
 #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 +72,9 @@ void register_default_materials() {
   auto& factory = material_factory<Traits>::instance();
   factory.template register_type<scalar_stepper<Traits>>("scalar_stepper");
   factory.template register_type<linear_elasticity<Traits>>("linear_elasticity");
+  factory.template register_type<constant_scalar<Traits>>("constant_scalar");
+  factory.template register_type<isotropic_tangent<Traits>>("isotropic_tangent");
+  factory.template register_type<linear_stress<Traits>>("linear_stress");
   factory.template register_type<autocatalytic_reaction<Traits>>("autocatalytic_reaction");
   factory.template register_type<backward_euler<Traits>>("backward_euler");
   factory.template register_type<tensor_component_stepper<1, Traits>>("tensor_component_stepper_rank1");

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 <typename Traits>
+class constant_scalar final
+    : public material_base<constant_scalar<Traits>, Traits> {
+public:
+  using base = material_base<constant_scalar<Traits>, Traits>;
+  using value_type = typename base::value_type;
+  using input_parameter_controller = typename base::input_parameter_controller;
+
+  template <typename... Args>
+  explicit constant_scalar(Args&&... args)
+      : base(std::forward<Args>(args)...),
+        m_value(base::template add_output<value_type>("value")) {
+    m_value = base::template get_parameter<value_type>("value");
+  }
+
+  static input_parameter_controller parameters() {
+    input_parameter_controller para{base::parameters()};
+    para.template insert<value_type>("value").template add<is_required>();
+    return para;
+  }
+
+private:
+  value_type& m_value;
+};
+
+}  // namespace numsim::materials
+
+#endif  // NUMSIM_MATERIALS_CONSTANT_SCALAR_H

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 <typename Traits>
 class isotropic_damage final
     : public material_base<isotropic_damage<Traits>, Traits> {
@@ -50,8 +52,13 @@ class isotropic_damage final
         // inputs
         m_stress(base::template add_input<tensor2>(
             m_elastic_source, "stress", EdgeKind::Global)),
+        // 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<tensor4>(
-            m_elastic_source, "tangent", EdgeKind::Global)),
+            base::m_parameter_handler.contains("tangent_source")
+                ? base::template get_parameter<std::string>("tangent_source")
+                : m_elastic_source,
+            "tangent", EdgeKind::Global)),
         m_damage(base::template add_input<value_type>(
             m_damage_source, "damage", EdgeKind::Global)),
         m_d_damage(base::template add_input<value_type>(
@@ -65,6 +72,8 @@ class isotropic_damage final
   static input_parameter_controller parameters() {
     input_parameter_controller para{base::parameters()};
     para.template insert<std::string>("elastic_source").template add<is_required>();
+    // No check == optional; declared so the JSON schema still knows the key.
+    para.template insert<std::string>("tangent_source");
     para.template insert<std::string>("damage_source").template add<is_required>();
     para.template insert<std::string>("state_source").template add<is_required>();
     para.template insert<std::string>("yield_source").template add<is_required>();

include/numsim-materials/materials/isotropic_tangent.h

@@ -0,0 +1,85 @@
+#ifndef NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H
+#define NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H
+
+#include <tmech/tmech.h>
+#include "numsim-materials/core/material_base.h"
+#include "numsim-materials/materials/plasticity_utils.h"
+
+namespace numsim::materials {
+
+/// 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.
+///
+/// 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 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.
+///
+/// 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 <typename Traits>
+class isotropic_tangent final
+    : public material_base<isotropic_tangent<Traits>, Traits> {
+public:
+  using base = material_base<isotropic_tangent<Traits>, Traits>;
+  using value_type = typename base::value_type;
+  using input_parameter_controller = typename base::input_parameter_controller;
+  using base::Dim;
+  using tensor4 = tmech::tensor<value_type, Dim, 4>;
+
+  template <typename... Args>
+  explicit isotropic_tangent(Args&&... args)
+      : base(std::forward<Args>(args)...),
+        m_C(base::template add_output<tensor4>(
+            "tangent", &isotropic_tangent::update_tangent)),
+        m_K(base::template add_input<value_type>(
+            base::template get_parameter<std::string>("K_source"),
+            base::template get_parameter<std::string>("K_property"),
+            EdgeKind::Global)),
+        m_G(base::template add_input<value_type>(
+            base::template get_parameter<std::string>("G_source"),
+            base::template get_parameter<std::string>("G_property"),
+            EdgeKind::Global)) {}
+
+  static input_parameter_controller parameters() {
+    input_parameter_controller para{base::parameters()};
+    para.template insert<std::string>("K_source").template add<is_required>();
+    para.template insert<std::string>("G_source").template add<is_required>();
+    // 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<std::string>("K_property")
+        .template add<set_default>(std::string{"value"});
+    para.template insert<std::string>("G_property")
+        .template add<set_default>(std::string{"value"});
+    return para;
+  }
+
+  void update_tangent() {
+    const auto I{tmech::eye<value_type, Dim, 2>()};
+    const auto IIvol{tmech::otimes(I, I) / Dim};
+    m_C = 3 * m_K.get() * IIvol +
+          2 * m_G.get() * plasticity_detail::make_IIdev<value_type, Dim>();
+  }
+
+private:
+  tensor4& m_C;
+  const input_property<value_type, property_traits>& m_K;
+  const input_property<value_type, property_traits>& m_G;
+};
+
+}  // namespace numsim::materials
+
+#endif  // NUMSIM_MATERIALS_ISOTROPIC_TANGENT_H

include/numsim-materials/materials/linear_stress.h

@@ -0,0 +1,59 @@
+#ifndef NUMSIM_MATERIALS_LINEAR_STRESS_H
+#define NUMSIM_MATERIALS_LINEAR_STRESS_H
+
+#include <tmech/tmech.h>
+#include "numsim-materials/core/material_base.h"
+
+namespace numsim::materials {
+
+/// sigma = C : eps, with C supplied by another material.
+///
+/// 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 with any tangent generator. linear_elasticity stays the better choice
+/// when the moduli are fixed — one fewer material, tangent computed once.
+template <typename Traits>
+class linear_stress final
+    : public material_base<linear_stress<Traits>, Traits> {
+public:
+  using base = material_base<linear_stress<Traits>, Traits>;
+  using value_type = typename base::value_type;
+  using input_parameter_controller = typename base::input_parameter_controller;
+  using base::Dim;
+  using tensor2 = tmech::tensor<value_type, Dim, 2>;
+  using tensor4 = tmech::tensor<value_type, Dim, 4>;
+
+  template <typename... Args>
+  explicit linear_stress(Args&&... args)
+      : base(std::forward<Args>(args)...),
+        m_sig(base::template add_output<tensor2>(
+            "stress", &linear_stress::update_stress)),
+        m_C(base::template add_input<tensor4>(
+            base::template get_parameter<std::string>("tangent_source"),
+            "tangent", EdgeKind::Global)),
+        m_eps(base::template add_input<tensor2>(
+            base::template get_parameter<std::string>("strain_source"),
+            "strain", EdgeKind::Global)) {}
+
+  static input_parameter_controller parameters() {
+    input_parameter_controller para{base::parameters()};
+    para.template insert<std::string>("tangent_source")
+        .template add<is_required>();
+    para.template insert<std::string>("strain_source")
+        .template add<is_required>();
+    return para;
+  }
+
+  void update_stress() { m_sig = tmech::dcontract(m_C.get(), m_eps.get()); }
+
+private:
+  tensor2& m_sig;
+  const input_property<tensor4, property_traits>& m_C;
+  const input_property<tensor2, property_traits>& m_eps;
+};
+
+}  // namespace numsim::materials
+
+#endif  // NUMSIM_MATERIALS_LINEAR_STRESS_H

include/numsim-materials/materials/weighted_sum.h

@@ -1,6 +1,7 @@
 #ifndef WEIGHTED_SUM_H
 #define WEIGHTED_SUM_H
 
+#include <stdexcept>
 #include <tmech/tmech.h>
 #include "numsim-materials/core/material_base.h"
 
@@ -18,6 +19,8 @@ namespace numsim::materials {
 /// Parameters:
 ///   "name":  material name
 ///   "terms": vector<pair<string,string>> — (weight_name, constituent_name) pairs
+///   "tangent_sources": optional vector<string> — one per term, naming a
+///       different material for that term's tangent; "" keeps the term's own.
 template <typename Traits>
 class weighted_sum final
     : public material_base<weighted_sum<Traits>, Traits> {
@@ -38,20 +41,44 @@ class weighted_sum final
         m_terms_param(base::template get_parameter<terms_type>("terms")),
         m_weight_property(base::template get_parameter<std::string>("weight_property")),
         m_stress_property(base::template get_parameter<std::string>("stress_property")),
-        m_tangent_property(base::template get_parameter<std::string>("tangent_property"))
+        m_tangent_property(base::template get_parameter<std::string>("tangent_property")),
+        m_tangent_sources(base::m_parameter_handler.contains("tangent_sources")
+                              ? base::template get_parameter<std::vector<std::string>>("tangent_sources")
+                              : std::vector<std::string>{})
   {
+    // 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) {
+      const bool overridden =
+          i < m_tangent_sources.size() && !m_tangent_sources[i].empty();
+      const std::string& tangent_owner =
+          overridden ? m_tangent_sources[i] : mat_name;
       auto& w = base::template add_input<value_type>(weight_name, m_weight_property, EdgeKind::Global);
       auto& s = base::template add_input<tensor2>(mat_name, m_stress_property, EdgeKind::Global);
-      auto& c = base::template add_input<tensor4>(mat_name, m_tangent_property, EdgeKind::Global);
+      auto& c = base::template add_input<tensor4>(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()};
+    // 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<std::vector<std::string>>("tangent_sources");
     para.template insert<terms_type>("terms").template add<is_required>();
     para.template insert<std::string>("weight_property")
         .template add<set_default>(std::string{"value"});
@@ -94,6 +121,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<std::string> m_tangent_sources;
   std::vector<term> m_terms;
 };