materials: split J2 from Drucker-Prager, and give plasticity its own elastic tangent

include/numsim-materials/materials/small_strain_plasticity.h → include/numsim-materials/materials/drucker_prager_plasticity.h

@@ -1,5 +1,5 @@
-#ifndef NUMSIM_MATERIALS_SMALL_STRAIN_PLASTICITY_H
-#define NUMSIM_MATERIALS_SMALL_STRAIN_PLASTICITY_H
+#ifndef NUMSIM_MATERIALS_DRUCKER_PRAGER_PLASTICITY_H
+#define NUMSIM_MATERIALS_DRUCKER_PRAGER_PLASTICITY_H
 
 #include <cmath>
 #include <concepts>
@@ -8,57 +8,42 @@
 #include <tmech/tmech.h>
 #include "numsim-materials/core/material_base.h"
 #include "numsim-materials/core/material_ref.h"
-#include "numsim-materials/materials/yield_functions.h"
 #include "numsim-materials/materials/drucker_prager_yield_function.h"
 #include "numsim-materials/materials/plasticity_utils.h"
 #include "numsim-materials/solvers/backward_euler.h"
 
 namespace numsim::materials {
 
-/// Concept for yield functions that support an apex return branch.
-/// All five methods must be present; checking a single sentinel is insufficient.
-template<typename YF, typename T, std::size_t Dim>
-concept has_apex_return = requires(const YF& yf,
-    const tmech::tensor<T, Dim, 2>& t2, T v) {
-  { yf.needs_apex_return(v, v, v) }      -> std::convertible_to<bool>;
-  { yf.apex_modified_sig_eq(t2) }        -> std::convertible_to<T>;
-  { yf.apex_effective_modulus() }         -> std::convertible_to<T>;
-  { yf.apex_plastic_strain(t2, t2, v) };
-  { yf.apex_tangent(v) };
-};
-
 /// Single-stage implicit Euler plasticity (classical return mapping).
 ///
 /// Uses solver.solve() for the Newton iteration. No tableau, no stage
 /// vectors, no overhead. This is the standard radial return for J2.
-template<typename Traits, typename YieldFunction>
-class small_strain_plasticity final
-    : public material_base<small_strain_plasticity<Traits, YieldFunction>, Traits> {
+template<typename Traits>
+class drucker_prager_plasticity final
+    : public material_base<drucker_prager_plasticity<Traits>, Traits> {
 public:
-  using base = material_base<small_strain_plasticity<Traits, YieldFunction>, Traits>;
+  using base = material_base<drucker_prager_plasticity<Traits>, Traits>;
   using value_type = typename base::value_type;
   using input_parameter_controller = typename base::input_parameter_controller;
   static constexpr auto Dim = base::Dim;
   using tensor2 = tmech::tensor<value_type, Dim, 2>;
   using tensor4 = tmech::tensor<value_type, Dim, 4>;
-  using yield_fn = YieldFunction;
+  using yield_fn = drucker_prager_yield_function<value_type, base::Dim>;
   using solver_type = backward_euler<Traits>;
 
   template <typename... Args>
-  explicit small_strain_plasticity(Args&&... args)
+  explicit drucker_prager_plasticity(Args&&... args)
       : base(std::forward<Args>(args)...),
         m_stress(base::template add_output<tensor2>(
-            "stress", &small_strain_plasticity::compute)),
+            "stress", &drucker_prager_plasticity::compute)),
         m_tangent(base::template add_output<tensor4>("tangent")),
         m_eps_p(base::template add_history_output<tensor2>("plastic_strain")),
         m_kappa(base::template add_history_output<value_type>("equivalent_plastic_strain")),
         m_G(base::template get_parameter<value_type>("G")),
         m_sigma_0(base::template get_parameter<value_type>("sigma_0")),
+        m_K_bulk(base::template get_parameter<value_type>("K_bulk")),
         m_solver(base::template add_material_ref<solver_type>(
             base::template get_parameter<std::string>("solver_source"))),
-        m_C_e(base::template add_input<tensor4>(
-            base::template get_parameter<std::string>("elastic_source"),
-            "tangent", EdgeKind::Global)),
         m_strain(base::template add_input<tensor2>(
             base::template get_parameter<std::string>("strain_source"),
             "strain", EdgeKind::Global)),
@@ -69,23 +54,34 @@ class small_strain_plasticity final
             base::template get_parameter<std::string>("hardening_source"),
             "hardening_modulus", EdgeKind::Local))
   {
-    if (base::m_parameter_handler.contains("yield_function"))
-      m_yf = base::template get_parameter<yield_fn>("yield_function");
+    m_yf = yield_fn(base::template get_parameter<value_type>("eta"),
+                    base::template get_parameter<value_type>("beta"),
+                    m_K_bulk);
+    const auto I = tmech::eye<value_type, Dim, 2>();
+    const tensor4 IIvol{tmech::otimes(I, I) / value_type{Dim}};
+    m_C_e = value_type{3} * m_K_bulk * IIvol +
+            value_type{2} * m_G * plasticity_detail::make_IIdev<value_type, Dim>();
   }
 
   static input_parameter_controller parameters() {
     input_parameter_controller para{base::parameters()};
-    para.template insert<std::string>("elastic_source").template add<is_required>();
     para.template insert<std::string>("hardening_source").template add<is_required>();
     para.template insert<std::string>("strain_source").template add<is_required>();
     para.template insert<std::string>("solver_source").template add<is_required>();
     para.template insert<value_type>("G").template add<is_required>();
     para.template insert<value_type>("sigma_0").template add<is_required>();
+    // The cone's friction, dilatancy and bulk modulus, as plain scalars.
+    // They used to arrive inside a "yield_function" C++ object, which the JSON
+    // reader cannot convert -- so the material could not be configured from a
+    // document at all (see #33).
+    para.template insert<value_type>("eta").template add<is_required>();
+    para.template insert<value_type>("beta").template add<is_required>();
+    para.template insert<value_type>("K_bulk").template add<is_required>();
     return para;
   }
 
   void compute() {
-    const auto& C_e = m_C_e.get();
+    const auto& C_e = m_C_e;
     const auto kappa_n = m_kappa.old_value();
 
     m_kappa.new_value() = kappa_n;
@@ -104,28 +100,21 @@ class small_strain_plasticity final
     // Conservative pre-check using the zero-hardening dlambda bound:
     //   dlambda_max = F_trial / G_eff  ≥  true dlambda  (for H' ≥ 0)
     // If even this upper bound triggers apex, smooth will also.
-    if constexpr (has_apex_return<yield_fn, value_type, Dim>) {
-      const auto G_eff = m_yf.effective_modulus(m_G);
-      const auto dlambda_max = ts.eval.F / G_eff;
-      if (m_yf.needs_apex_return(m_G, dlambda_max, ts.eval.sig_eq)) {
-        do_apex_return(ts.eval.sig, C_e, kappa_n);
-        return;
-      }
+    const auto G_eff_pre = m_yf.effective_modulus(m_G);
+    if (m_yf.needs_apex_return(m_G, ts.eval.F / G_eff_pre, ts.eval.sig_eq)) {
+      do_apex_return(ts.eval.sig, C_e, kappa_n);
+      return;
     }
 
     const auto dlambda = solve_smooth_newton(ts.eval.modified_sig_eq, kappa_n);
 
     // If smooth Newton fails and apex is available, try apex as fallback.
     if (!m_solver.get().converged()) {
-      if constexpr (has_apex_return<yield_fn, value_type, Dim>) {
-        do_apex_return(ts.eval.sig, C_e, kappa_n);
-        if (!m_solver.get().converged())
-          throw std::runtime_error(
-              "small_strain_plasticity: both smooth and apex Newton failed");
-        return;
-      }
-      throw std::runtime_error(
-          "small_strain_plasticity: smooth return-mapping Newton failed");
+      do_apex_return(ts.eval.sig, C_e, kappa_n);
+      if (!m_solver.get().converged())
+        throw std::runtime_error(
+            "drucker_prager_plasticity: both smooth and apex Newton failed");
+      return;
     }
 
     do_smooth_return(ts.eval, C_e, kappa_n, dlambda);
@@ -163,22 +152,28 @@ class small_strain_plasticity final
                         const tensor4& C_e, value_type kappa_n, value_type dlambda) {
     m_eps_p.new_value() = m_eps_p.old_value() + dlambda * ts.N;
     m_kappa.new_value() = kappa_n + dlambda;
-    m_stress = tmech::dcontract(C_e, m_strain.get() - m_eps_p.new_value());
+
+    // sigma = C_e : (eps - eps_p_new) = sig_trial - dl (C_e : N).
+    // For isotropic C_e that is 2G dev(N) + K tr(N) I -- no rank-4 : rank-2
+    // contraction, which measures 33 ns against ~720 for the whole step. The
+    // flow is non-associative, so tr(N) = beta is generally nonzero and the
+    // volumetric term does not drop out as it does for J2.
+    const auto I = tmech::eye<value_type, Dim, 2>();
+    const tensor2 Ce_N{value_type{2} * m_G * tmech::dev(ts.N) +
+                       m_K_bulk * tmech::trace(ts.N) * I};
+    m_stress = ts.sig - dlambda * Ce_N;
 
     // Tangent at trial state (return mapping uses N_trial).
     // For J2, trial = converged. For DP, they differ.
     m_H.update_source();
     m_tangent = plasticity_detail::compute_tangent<value_type, Dim>(
-        m_yf, ts.sig_dev, ts.N, ts.sig_eq, dlambda, m_dH.get(), C_e);
+        m_yf, ts.sig_dev, ts.N, ts.sig_eq, dlambda, m_dH.get(), C_e, m_G);
   }
 
   /// Apex return: deviatoric stress vanishes, only volumetric Newton.
   /// dev(ε_p) = dev(ε), tr(ε_p) += β·Δκ. Tangent is rank-1 volumetric.
-  /// Compiled only when the yield function provides apex support.
   void do_apex_return(const tensor2& sig_trial, const tensor4& C_e,
-                      value_type kappa_n)
-    requires has_apex_return<yield_fn, value_type, Dim>
-  {
+                      value_type kappa_n) {
     const auto phi_apex = m_yf.apex_modified_sig_eq(sig_trial);
     const auto G_eff_apex = m_yf.apex_effective_modulus();
     const auto dkappa = solve_scalar_return(phi_apex, G_eff_apex, kappa_n);
@@ -202,25 +197,23 @@ class small_strain_plasticity final
 
   const value_type& m_G;
   const value_type& m_sigma_0;
+  const value_type& m_K_bulk;
   material_ref<solver_type, Traits>& m_solver;
 
-  const input_property<tensor4, property_traits>& m_C_e;
   const input_property<tensor2, property_traits>& m_strain;
   const input_property<value_type, property_traits>& m_H;
   const input_property<value_type, property_traits>& m_dH;
   yield_fn m_yf{};
-};
 
-template<typename Traits>
-using j2_plasticity = small_strain_plasticity<Traits,
-    j2_yield_function<typename Traits::value_type, Traits::Dim>>;
+  /// The elastic stiffness, built here rather than read from another material.
+  /// The tangent collapse requires an isotropic C_e -- C_e : X = 2G X for
+  /// deviatoric X -- so accepting an arbitrary rank-4 tangent advertised a
+  /// generality this material cannot honour. K_bulk and G are already
+  /// parameters, so nothing new is asked of the caller.
+  tensor4 m_C_e{};
+};
 
-/// Drucker-Prager plasticity. The yield function (with η, β, K_bulk) must be
-/// supplied via the "yield_function" parameter at construction.
-template<typename Traits>
-using drucker_prager_plasticity = small_strain_plasticity<Traits,
-    drucker_prager_yield_function<typename Traits::value_type, Traits::Dim>>;
 
 } // namespace numsim::materials
 
-#endif // NUMSIM_MATERIALS_SMALL_STRAIN_PLASTICITY_H
+#endif // NUMSIM_MATERIALS_DRUCKER_PRAGER_PLASTICITY_H