materials: plasticity builds its own elastic tangent

include/numsim-materials/materials/drucker_prager_plasticity.h

@@ -41,11 +41,9 @@ class drucker_prager_plasticity final
         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)),
@@ -58,12 +56,15 @@ class drucker_prager_plasticity final
   {
     m_yf = yield_fn(base::template get_parameter<value_type>("eta"),
                     base::template get_parameter<value_type>("beta"),
-                    base::template get_parameter<value_type>("K_bulk"));
+                    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>();
@@ -80,7 +81,7 @@ class drucker_prager_plasticity final
   }
 
   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;
@@ -151,7 +152,16 @@ class drucker_prager_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.
@@ -187,13 +197,20 @@ class drucker_prager_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{};
+
+  /// 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{};
 };
 
 

include/numsim-materials/materials/j2_plasticity.h

@@ -45,9 +45,10 @@ namespace numsim::materials {
 /// may be nonlinear (see exponential_isotropic_hardening).
 ///
 /// Parameters:
-///   "name", "elastic_source", "hardening_source", "strain_source",
-///   "solver_source", "G", "sigma_0"
-/// -- the same set small_strain_plasticity takes, so this is a drop-in.
+///   "name", "hardening_source", "strain_source", "solver_source",
+///   "K", "G", "sigma_0"
+///
+/// No elastic_source: the stiffness is built from K and G here.
 template <typename Traits>
 class j2_plasticity final
     : public material_base<j2_plasticity<Traits>, Traits> {
@@ -71,11 +72,9 @@ class j2_plasticity final
             "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(base::template get_parameter<value_type>("K")),
         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)),
@@ -85,26 +84,41 @@ class j2_plasticity final
         m_dH(base::template add_input<value_type>(
             base::template get_parameter<std::string>("hardening_source"),
             "hardening_modulus", EdgeKind::Local)),
-        m_IIdev(plasticity_detail::make_IIdev<value_type, Dim>())
+        m_IIdev(plasticity_detail::make_IIdev<value_type, Dim>()),
+        m_C_e(build_elastic_tangent(m_K, m_G, m_IIdev))
   {}
 
+  /// The elastic stiffness, built here rather than read from another material.
+  ///
+  /// The closed forms below REQUIRE an isotropic C_e -- that is what makes
+  /// C_e : N = 2G N and N : C_e : N = 3G true. Accepting an arbitrary rank-4
+  /// tangent from an elastic_source advertised a generality this material
+  /// cannot honour, and dragged a linear_elasticity into every plasticity graph
+  /// whose own "stress" output (C : eps, ignoring eps_p) is meaningless once
+  /// yielding starts.
+  static tensor4 build_elastic_tangent(value_type K, value_type G,
+                                       const tensor4& IIdev) {
+    const auto I = tmech::eye<value_type, Dim, 2>();
+    return value_type{3} * K * (tmech::otimes(I, I) / value_type{Dim}) +
+           value_type{2} * G * IIdev;
+  }
+
   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>("K").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>();
     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;
@@ -143,8 +157,13 @@ class j2_plasticity final
 
     m_eps_p.new_value() = m_eps_p.old_value() + dlambda * N;
     m_kappa.new_value() = kappa_n + dlambda;
-    m_stress = tmech::dcontract(
-        C_e, tensor2(m_strain.get() - m_eps_p.new_value()));
+
+    // sigma = C_e : (eps - eps_p_new)
+    //       = C_e : (eps - eps_p_old) - dl (C_e : N)
+    //       = sig_trial - 2G dl N,           since C_e : N = 2G N.
+    // The same identity the tangent uses. Avoids a second rank-4 : rank-2
+    // contraction, which measures 33 ns against ~245 for the whole step.
+    m_stress = sig_trial - value_type{2} * m_G * dlambda * N;
 
     m_H.update_source();
     const auto GG = m_G * m_G;
@@ -162,14 +181,15 @@ class j2_plasticity final
 
   const value_type& m_G;
   const value_type& m_sigma_0;
+  const value_type& m_K;
   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;
 
   const tensor4 m_IIdev;
+  const tensor4 m_C_e;
 };
 
 }  // namespace numsim::materials

include/numsim-materials/materials/rk_plasticity.h

@@ -37,23 +37,26 @@ class rk_plasticity final
         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_K(base::template get_parameter<value_type>("K")),
         m_G(base::template get_parameter<value_type>("G")),
         m_sigma_0(base::template get_parameter<value_type>("sigma_0")),
         m_tol(base::template get_parameter<value_type>("tolerance")),
         m_max_iter(base::template get_parameter<int>("max_iter")),
         m_tableau(base::template get_parameter<const butcher_tableau*>("tableau")),
-        m_elastic_source(base::template get_parameter<std::string>("elastic_source")),
         m_hardening_source(base::template get_parameter<std::string>("hardening_source")),
         m_strain_source(base::template get_parameter<std::string>("strain_source")),
-        m_C_e(base::template add_input<tensor4>(
-            m_elastic_source, "tangent", EdgeKind::Global)),
         m_strain(base::template add_input<tensor2>(
             m_strain_source, "strain", EdgeKind::Global)),
         m_H(base::template add_input<value_type>(
             m_hardening_source, "hardening_stress", EdgeKind::Local)),
         m_dH(base::template add_input<value_type>(
             m_hardening_source, "hardening_modulus", EdgeKind::Local))
   {
+    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 * IIvol +
+            value_type{2} * m_G * plasticity_detail::make_IIdev<value_type, Dim>();
+
     const int s = m_tableau->stages();
     m_dlambda.resize(s, value_type{0});
     m_N_stage.resize(s);
@@ -67,9 +70,9 @@ class rk_plasticity final
 
   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<value_type>("K").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>();
     para.template insert<value_type>("tolerance")
@@ -80,7 +83,7 @@ class rk_plasticity final
   }
 
   void compute() {
-    const auto& C_e = m_C_e.get();
+    const auto& C_e = m_C_e;
     const auto& eps = m_strain.get();
     const auto kappa_n = m_kappa.old_value();
     const auto eps_p_n = m_eps_p.old_value();
@@ -181,20 +184,26 @@ class rk_plasticity final
   history_property<tensor2>& m_eps_p;
   history_property<value_type>& m_kappa;
 
+  const value_type& m_K;
   const value_type& m_G;
   const value_type& m_sigma_0;
   const value_type& m_tol;
   const int& m_max_iter;
   const butcher_tableau* m_tableau;
-  const std::string& m_elastic_source;
   const std::string& m_hardening_source;
   const std::string& m_strain_source;
 
-  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;
 
+  /// Built here, not read from an elastic_source. compute_tangent's collapse
+  /// (C_e : X = 2G X) and effective_modulus(G) = 3G both require an isotropic
+  /// C_e, so taking an arbitrary rank-4 tangent promised more than this class
+  /// can deliver -- and pulled a linear_elasticity into the graph whose own
+  /// "stress" output is meaningless once eps_p is nonzero.
+  tensor4 m_C_e{};
+
   yield_fn m_yf{};
   std::vector<value_type> m_dlambda;
   std::vector<tensor2> m_N_stage;

tests/plot_data.cpp

@@ -78,10 +78,10 @@ run_result run_j2(T increment, int steps,
 
   p.clear();
   p.insert<std::string>("name", "j2");
-  p.insert<std::string>("elastic_source", "elastic");
   p.insert<std::string>("hardening_source", "hardening");
   p.insert<std::string>("strain_source", "stepper");
   p.insert<std::string>("solver_source", "solver");
+  p.insert<T>("K", K_val);
   p.insert<T>("G", G_val);
   p.insert<T>("sigma_0", sigma_0);
   ctx.create<j2_plasticity>(p);
@@ -144,7 +144,6 @@ run_result run_dp(T increment, int steps,
 
   p.clear();
   p.insert<std::string>("name", "dp");
-  p.insert<std::string>("elastic_source", "elastic");
   p.insert<std::string>("hardening_source", "hardening");
   p.insert<std::string>("strain_source", "stepper");
   p.insert<std::string>("solver_source", "solver");

tests/test_drucker_prager.cpp

@@ -116,7 +116,6 @@ class DruckerPragerTest : public ::testing::Test {
 
     p.clear();
     p.insert<std::string>("name", "dp");
-    p.insert<std::string>("elastic_source", "elastic");
     p.insert<std::string>("hardening_source", "hardening");
     p.insert<std::string>("strain_source", "stepper");
     p.insert<std::string>("solver_source", "solver");
@@ -221,7 +220,6 @@ class DPTangentTest : public ::testing::Test {
 
     p.clear();
     p.insert<std::string>("name", "dp");
-    p.insert<std::string>("elastic_source", "elastic");
     p.insert<std::string>("hardening_source", "hardening");
     p.insert<std::string>("strain_source", "stepper");
     p.insert<std::string>("solver_source", "solver");
@@ -295,7 +293,6 @@ T run_dp_max_tangent_error(T increment, int steps) {
 
   p.clear();
   p.insert<std::string>("name", "dp");
-  p.insert<std::string>("elastic_source", "elastic");
   p.insert<std::string>("hardening_source", "hardening");
   p.insert<std::string>("strain_source", "stepper");
   p.insert<std::string>("solver_source", "solver");
@@ -383,7 +380,6 @@ T max_tangent_error(std::vector<double> direction, T increment, int steps) {
 
   p.clear();
   p.insert<std::string>("name", "dp");
-  p.insert<std::string>("elastic_source", "elastic");
   p.insert<std::string>("hardening_source", "hardening");
   p.insert<std::string>("strain_source", "stepper");
   p.insert<std::string>("solver_source", "solver");
@@ -481,7 +477,6 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) {
 
   p.clear();
   p.insert<std::string>("name", "dp");
-  p.insert<std::string>("elastic_source", "elastic");
   p.insert<std::string>("hardening_source", "hardening");
   p.insert<std::string>("strain_source", "stepper");
   p.insert<std::string>("solver_source", "solver");
@@ -540,7 +535,6 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) {
   ctx.create<numsim::materials::linear_isotropic_hardening<policy>>(p);
   p.clear();
   p.insert<std::string>("name", "dp");
-  p.insert<std::string>("elastic_source", "elastic");
   p.insert<std::string>("hardening_source", "hardening");
   p.insert<std::string>("strain_source", "stepper");
   p.insert<std::string>("solver_source", "solver");

tests/test_j2_plasticity.cpp

@@ -56,10 +56,10 @@ class J2PlasticityTest : public ::testing::Test {
     // J2 plasticity — solver passed as pointer
     p.clear();
     p.insert<std::string>("name", "j2");
-    p.insert<std::string>("elastic_source", "elastic");
-    p.insert<std::string>("hardening_source", "hardening");
+      p.insert<std::string>("hardening_source", "hardening");
     p.insert<std::string>("strain_source", "stepper");
     p.insert<std::string>("solver_source", "solver");
+    p.insert<T>("K", K);
     p.insert<T>("G", G);
     p.insert<T>("sigma_0", sigma_0);
     ctx.create<numsim::materials::j2_plasticity<policy>>(p);
@@ -168,10 +168,10 @@ class J2TangentTest : public ::testing::Test {
 
     p.clear();
     p.insert<std::string>("name", "j2");
-    p.insert<std::string>("elastic_source", "elastic");
-    p.insert<std::string>("hardening_source", "hardening");
+      p.insert<std::string>("hardening_source", "hardening");
     p.insert<std::string>("strain_source", "stepper");
     p.insert<std::string>("solver_source", "solver");
+    p.insert<T>("K", T{166.67});
     p.insert<T>("G", T{76.92});
     p.insert<T>("sigma_0", T{50.0});
     ctx.create<numsim::materials::j2_plasticity<policy>>(p);
@@ -258,9 +258,9 @@ class RKPlasticityTest : public ::testing::Test {
 
     p.clear();
     p.insert<std::string>("name", "j2");
-    p.insert<std::string>("elastic_source", "elastic");
     p.insert<std::string>("hardening_source", "hardening");
     p.insert<std::string>("strain_source", "stepper");
+    p.insert<T>("K", T{166.67});
     p.insert<T>("G", T{76.92});
     p.insert<T>("sigma_0", T{50.0});
     p.insert<const numsim::materials::butcher_tableau*>("tableau", &m_tab);