Add vector_newton: coupled Newton for mixed scalar/tensor local systems

CMakeLists.txt

@@ -25,14 +25,29 @@ FetchContent_Declare(
 FetchContent_MakeAvailable(numsim-core)
 
 # --- Dependencies (tmech tensor library) ---
-find_package(tmech QUIET)
-if(NOT tmech_FOUND)
+# Fetched by default rather than found. find_package(tmech) succeeds against
+# ANY installed version, including one predating features used here (the Mandel
+# adaptor tag in solvers/unknown_layout.h), and tmech's config exposes no
+# version constraint that could express the difference — so a stale system
+# install would silently win and fail at compile time.
+option(NUMSIM_USE_SYSTEM_TMECH "Use an installed tmech instead of fetching it" OFF)
+
+if(NUMSIM_USE_SYSTEM_TMECH)
+  find_package(tmech REQUIRED)
+else()
+  # tmech defaults TMECH_BUILD_TESTS and TMECH_BUILD_EXAMPLES to ON; as a
+  # subproject we want neither. (TMECH_INSTALL already defaults off when tmech
+  # is not the top-level project.)
+  set(TMECH_BUILD_TESTS OFF CACHE BOOL "" FORCE)
+  set(TMECH_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE)
+  set(TMECH_BUILD_BENCHMARK OFF CACHE BOOL "" FORCE)
   FetchContent_Declare(
     tmech
     GIT_REPOSITORY https://github.com/petlenz/tmech
     GIT_TAG master
   )
-  set(BUILD_TESTS OFF CACHE BOOL "" FORCE)
+  # To use a local tmech checkout, configure with:
+  #   -DFETCHCONTENT_SOURCE_DIR_TMECH=/path/to/tmech
   FetchContent_MakeAvailable(tmech)
 endif()
 

README.md

@@ -70,6 +70,7 @@ cmake --build build -j$(nproc)
 | `NUMSIM_BUILD_TESTS` | ON | Build GTest unit tests |
 | `NUMSIM_BUILD_EXAMPLES` | ON | Build examples |
 | `ENABLE_PLOTTING` | OFF | Qt6 + QCustomPlot live plotting |
+| `NUMSIM_USE_SYSTEM_TMECH` | OFF | Use an installed tmech instead of fetching it |
 
 ### Run tests
 
@@ -83,7 +84,9 @@ cd build && ctest --output-on-failure
 
 - C++23 (GCC 14+ or Clang 18+)
 - [numsim-core](https://github.com/NumSim-Stack/numsim-core) — fetched automatically via CMake
-- [tmech](https://github.com/petlenz/tmech) — tensor library (`find_package(tmech REQUIRED)`)
+- [tmech](https://github.com/petlenz/tmech) — tensor library; fetched automatically via CMake
+  (set `NUMSIM_USE_SYSTEM_TMECH=ON` to use an installed copy instead, or
+  `-DFETCHCONTENT_SOURCE_DIR_TMECH=/path/to/tmech` for a local checkout)
 - [nlohmann/json](https://github.com/nlohmann/json) — for JSON configuration (optional)
 - Qt6 + QCustomPlot — for live plotting (optional)
 

include/numsim-materials/core/unknown_spec.h

@@ -0,0 +1,34 @@
+#ifndef NUMSIM_MATERIALS_UNKNOWN_SPEC_H
+#define NUMSIM_MATERIALS_UNKNOWN_SPEC_H
+
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace numsim::materials {
+
+/// Runtime description of the unknown set of a coupled local system.
+///
+/// Lives in core/ rather than solvers/ so that the JSON layer can read it
+/// without depending on the solver or on tmech. A bounded switch in
+/// vector_newton turns each spec into the corresponding compile-time kind in
+/// solvers/unknown_layout.h; the set is deliberately small and closed so that
+/// switch stays exhaustive.
+///
+/// Symmetry has to be declared here because it cannot be recovered from the
+/// C++ type: tmech's tensor<T,Dim,2> is symmetry-agnostic storage, and
+/// property_traits carries no shape metadata.
+enum class unknown_kind { scalar, sym_tensor };
+
+struct unknown_spec {
+  std::string name;
+  unknown_kind kind{unknown_kind::scalar};
+};
+
+/// Identifies a Jacobian block by the names of its row and column unknowns.
+/// Used to declare structurally-zero blocks, which are then never wired.
+using block_ref = std::pair<std::string, std::string>;
+
+} // namespace numsim::materials
+
+#endif // NUMSIM_MATERIALS_UNKNOWN_SPEC_H

include/numsim-materials/default_materials.h

@@ -4,6 +4,7 @@
 #include <numsim-core/object_registry.h>
 #include "numsim-materials/core/material_base.h"
 #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/linear_elasticity.h"
 #include "numsim-materials/materials/autocatalytic_reaction.h"
@@ -80,6 +81,9 @@ void register_default_materials() {
   factory.template register_type<strain_threshold_yield<Traits>>("strain_threshold_yield");
   factory.template register_type<exponential_damage_law<Traits>>("exponential_damage_law");
   factory.template register_type<isotropic_damage<Traits>>("isotropic_damage");
+  // A single registered type covers every unknown combination — the layout is
+  // dispatched from the "unknowns" parameter, not baked into the template.
+  factory.template register_type<vector_newton<Traits>>("vector_newton");
 }
 
 } // namespace numsim::materials

include/numsim-materials/io/json_parameter_converter.h

@@ -13,6 +13,7 @@
 #include <variant>
 #include <vector>
 #include <numsim-core/input_parameter_controller.h>
+#include "numsim-materials/core/unknown_spec.h"
 
 namespace numsim::materials {
 
@@ -92,6 +93,7 @@ json_reader_registry<JsonType> make_default_json_registry() {
   reg.template add<double>();
   reg.template add<float>();
   reg.template add<int>();
+  reg.template add<bool>();
   reg.template add<std::string>();
 
   // Vectors
@@ -132,6 +134,48 @@ json_reader_registry<JsonType> make_default_json_registry() {
         return result;
       });
 
+  // unknown_spec list for coupled local systems (vector_newton):
+  //   "unknowns": [{"name": "dgamma", "kind": "scalar"},
+  //                {"name": "backstress", "kind": "sym_tensor"}]
+  // No "dim" — the dimension is fixed by the Traits policy, which keeps the
+  // kind set small enough for an exhaustive switch on the solver side.
+  reg.template add<std::vector<unknown_spec>>(
+      [](const JsonType& j) -> std::any {
+        std::vector<unknown_spec> result;
+        for (const auto& elem : j) {
+          unknown_spec s;
+          s.name = adapter::template get<std::string>(adapter::at(elem, "name"));
+          const auto kind =
+              adapter::template get<std::string>(adapter::at(elem, "kind"));
+          if (kind == "scalar")           s.kind = unknown_kind::scalar;
+          else if (kind == "sym_tensor")  s.kind = unknown_kind::sym_tensor;
+          else
+            throw std::runtime_error(
+                "unknown_spec: unrecognised kind '" + kind +
+                "' (expected \"scalar\" or \"sym_tensor\")");
+          result.push_back(std::move(s));
+        }
+        return result;
+      });
+
+  // Structurally-zero Jacobian blocks:
+  //   "zero_blocks": [["dgamma", "backstress"]]
+  reg.template add<std::vector<block_ref>>(
+      [](const JsonType& j) -> std::any {
+        std::vector<block_ref> result;
+        for (const auto& elem : j) {
+          std::vector<std::string> names;
+          for (const auto& part : elem)
+            names.push_back(adapter::template get<std::string>(part));
+          if (names.size() != 2)
+            throw std::runtime_error(
+                "zero_blocks: each entry must be a [row, column] pair of "
+                "unknown names");
+          result.emplace_back(std::move(names[0]), std::move(names[1]));
+        }
+        return result;
+      });
+
   return reg;
 }
 

include/numsim-materials/materials/small_strain_plasticity.h

@@ -150,7 +150,8 @@ class small_strain_plasticity final
       return {m_yf.residual(phi_trial, dl, G_eff, m_sigma_0, m_H.get()),
               m_yf.jacobian(G_eff, m_dH.get())};
     };
-    return m_solver.get().solve(eval);
+    // Non-negative: Δλ < 0 would be backward plastic flow (KKT).
+    return m_solver.get().solve_nonnegative(eval);
   }
 
   /// Smooth-cone return Newton: phi = modified_sig_eq, G_eff from yield function.

include/numsim-materials/solvers/backward_euler.h

@@ -83,19 +83,34 @@ class backward_euler final
   /// Direct call: another material provides eval(x) → {residual, jacobian}.
   /// Used when the caller drives the iteration (e.g., plasticity return mapping).
   /// Sets m_converged to indicate whether the iteration converged.
-  /// The returned value is clamped to be non-negative — for plasticity, a
-  /// negative plastic-multiplier increment is unphysical (backward plastic flow).
+  ///
+  /// This is a general scalar Newton solver: the root is returned as found,
+  /// including negative roots. Callers needing the plasticity KKT constraint
+  /// Δγ ≥ 0 must ask for it explicitly via solve_nonnegative().
   template<typename Eval>
   value_type solve(Eval&& eval, value_type x0 = value_type{0}) {
     auto x = x0;
     for (int i = 0; i < m_max_iter; ++i) {
       auto [r, dr] = eval(x);
-      if (std::abs(r) < m_tol) { m_converged = true; return std::max(x, value_type{0}); }
-      if (std::abs(dr) < value_type{1e-30}) { m_converged = false; return std::max(x, value_type{0}); }
+      if (std::abs(r) < m_tol) { m_converged = true; return x; }
+      if (std::abs(dr) < value_type{1e-30}) { m_converged = false; return x; }
       x -= r / dr;
     }
     m_converged = false;
-    return std::max(x, value_type{0});
+    return x;
+  }
+
+  /// solve() with the plasticity non-negativity projection applied.
+  ///
+  /// For a return mapping, a negative plastic-multiplier increment is
+  /// unphysical (backward plastic flow), so the converged root is clamped.
+  /// The clamp is deliberately NOT applied on the failure paths: a
+  /// non-converged iterate is returned raw so that a failed solve cannot be
+  /// mistaken for a plausible non-negative value. Check converged().
+  template<typename Eval>
+  value_type solve_nonnegative(Eval&& eval, value_type x0 = value_type{0}) {
+    const auto x = solve(std::forward<Eval>(eval), x0);
+    return m_converged ? std::max(x, value_type{0}) : x;
   }
 
   /// Whether the last solve() call converged.