From 9229131fd9e437e6b3df81266caef27ecc7df10b Mon Sep 17 00:00:00 2001 From: petlenz Date: Fri, 4 Sep 2026 22:38:14 +0200 Subject: [PATCH 1/6] materials: split J2 out of the generic plasticity class; 5x faster small_strain_plasticity is parameterised over a yield function so it can serve Drucker-Prager. J2 used none of that generality and paid for all of it: - a yield normal distinct from the flow normal, identical for associative J2 - a "modified" equivalent stress carrying pressure coupling J2 does not have - an apex branch a cylinder cannot reach The cost was not only readability. The general consistent tangent forms C_e : (dN/dsigma) : C_e -- two rank-4 x rank-4 contractions per plastic step. For isotropic elasticity N is deviatoric, so C_e : N = 2G N and N : C_e : N = 3G, and the expression collapses to the standard closed form: C = C_e - (6G^2 dl/sig_eq) IIdev + (4G^2 dl/sig_eq - 4G^2/(3G+H')) N (x) N algebraically identical to what the general path computes, not an approximation of it. old small_strain_plasticity : 1262.2 ns/step new j2_plasticity : 245.2 ns/step (5.15x) Equivalence over 60 steps on the same path: max |dstress| = 0.000e+00 (bit-identical) max |dalpha| = 0.000e+00 (bit-identical) max |dtangent| = 8.760e-14 (roundoff; different order of operations) All nine pre-existing J2 tests pass unchanged, including the tangent checker. ISOTROPY is not newly assumed. small_strain_plasticity already required it through effective_modulus(G) = 3G in its residual; the closed form states it instead of implying it. The scalar Newton is kept, because hardening is a separate graph node and may be nonlinear. Also fixes a blind spot the refactor exposed. J2TangentTest bounded the whole run by 0.1, because the elastic->plastic transition step is genuinely inexact (a central difference straddling the yield surface averages two tangents). That made it blind to real errors: a 0.5% error in the tangent lands at 2.4e-4 and passed. Fully plastic steps are now bounded at 1e-8 and the crossing step separately at 0.1. Verified: the mutation that passed before now fails. --- .../materials/j2_plasticity.h | 177 ++++++++++++++++++ .../materials/small_strain_plasticity.h | 7 +- tests/plot_data.cpp | 1 + tests/test_j2_plasticity.cpp | 40 +++- 4 files changed, 213 insertions(+), 12 deletions(-) create mode 100644 include/numsim-materials/materials/j2_plasticity.h diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h new file mode 100644 index 0000000..7322fd4 --- /dev/null +++ b/include/numsim-materials/materials/j2_plasticity.h @@ -0,0 +1,177 @@ +#ifndef NUMSIM_MATERIALS_J2_PLASTICITY_H +#define NUMSIM_MATERIALS_J2_PLASTICITY_H + +#include +#include +#include +#include +#include "numsim-materials/core/material_base.h" +#include "numsim-materials/core/material_ref.h" +#include "numsim-materials/materials/plasticity_utils.h" +#include "numsim-materials/solvers/backward_euler.h" + +namespace numsim::materials { + +/// J2 (von Mises) plasticity with isotropic hardening, radial return. +/// +/// Split out of small_strain_plasticity. That class +/// is parameterised over a yield function so it can also serve Drucker-Prager, +/// and J2 paid for the generality without using it: +/// +/// - a yield normal distinct from the flow normal, which for associative J2 +/// is the same tensor; +/// - a "modified" equivalent stress carrying pressure coupling J2 does not +/// have; +/// - an apex branch a cylinder cannot reach. +/// +/// The cost was not only readability. The general consistent tangent forms +/// C_e : (dN/dsigma) : C_e -- two rank-4 x rank-4 contractions per plastic +/// step. For isotropic elasticity N is deviatoric, so C_e : N = 2G N and +/// N : C_e : N = 3G, and the whole expression collapses: +/// +/// C = C_e - (6G^2 dl / sig_eq) IIdev +/// + (4G^2 dl / sig_eq - 4G^2 / (3G + H')) N (x) N +/// +/// which is the standard closed form. It is algebraically identical to what +/// the general path computes, not an approximation of it. +/// +/// ISOTROPY. The closed form uses C_e : N = 2G N. That is not a new +/// restriction: small_strain_plasticity already assumes isotropic elasticity +/// through effective_modulus(G) = 3G in its residual, so both paths are +/// equally limited. Here it is stated rather than implied. +/// +/// The scalar Newton is kept. Linear hardening converges in one iteration and +/// has a closed form, but the hardening material is a separate graph node and +/// 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. +template +class j2_plasticity final + : public material_base, Traits> { +public: + using base = material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + using base::Dim; + using tensor2 = tmech::tensor; + using tensor4 = tmech::tensor; + using solver_type = backward_euler; + + template + explicit j2_plasticity(Args&&... args) + : base(std::forward(args)...), + m_stress(base::template add_output( + "stress", &j2_plasticity::compute)), + m_tangent(base::template add_output("tangent")), + m_eps_p(base::template add_history_output("plastic_strain")), + m_kappa(base::template add_history_output( + "equivalent_plastic_strain")), + m_G(base::template get_parameter("G")), + m_sigma_0(base::template get_parameter("sigma_0")), + m_solver(base::template add_material_ref( + base::template get_parameter("solver_source"))), + m_C_e(base::template add_input( + base::template get_parameter("elastic_source"), + "tangent", EdgeKind::Global)), + m_strain(base::template add_input( + base::template get_parameter("strain_source"), + "strain", EdgeKind::Global)), + m_H(base::template add_input( + base::template get_parameter("hardening_source"), + "hardening_stress", EdgeKind::Local)), + m_dH(base::template add_input( + base::template get_parameter("hardening_source"), + "hardening_modulus", EdgeKind::Local)), + m_IIdev(plasticity_detail::make_IIdev()) + {} + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + para.template insert("elastic_source") + .template add(); + para.template insert("hardening_source") + .template add(); + para.template insert("strain_source") + .template add(); + para.template insert("solver_source") + .template add(); + para.template insert("G").template add(); + para.template insert("sigma_0").template add(); + return para; + } + + void compute() { + const auto& C_e = m_C_e.get(); + const auto kappa_n = m_kappa.old_value(); + + m_kappa.new_value() = kappa_n; + m_H.update_source(); + + // Trial state: freeze the plastic strain and load elastically. + const tensor2 sig_trial{ + tmech::dcontract(C_e, tensor2(m_strain.get() - m_eps_p.old_value()))}; + const tensor2 s{tmech::dev(sig_trial)}; + const auto sig_eq = + std::sqrt(value_type{1.5} * tmech::dcontract(s, s)); + + if (sig_eq - m_sigma_0 - m_H.get() <= value_type{0}) { + m_stress = sig_trial; + m_tangent = C_e; + m_eps_p.new_value() = m_eps_p.old_value(); + m_kappa.new_value() = kappa_n; + return; + } + + // Radial return: the flow direction is fixed by the trial state, so only + // the magnitude is solved for. + const tensor2 N{value_type{1.5} * s / sig_eq}; + const auto G_eff = value_type{3} * m_G; + + auto eval = [&](value_type dl) -> std::pair { + m_kappa.new_value() = kappa_n + dl; + m_H.update_source(); + return {sig_eq - G_eff * dl - m_sigma_0 - m_H.get(), + -G_eff - m_dH.get()}; + }; + const auto dlambda = m_solver.get().solve(eval); + if (!m_solver.get().converged()) + throw std::runtime_error( + "j2_plasticity: return-mapping Newton failed to converge"); + + 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())); + + m_H.update_source(); + const auto GG = m_G * m_G; + const auto a = value_type{6} * GG * dlambda / sig_eq; + const auto b = value_type{4} * GG * dlambda / sig_eq - + value_type{4} * GG / (G_eff + m_dH.get()); + m_tangent = C_e - a * m_IIdev + b * tmech::otimes(N, N); + } + +private: + tensor2& m_stress; + tensor4& m_tangent; + history_property& m_eps_p; + history_property& m_kappa; + + const value_type& m_G; + const value_type& m_sigma_0; + material_ref& m_solver; + + const input_property& m_C_e; + const input_property& m_strain; + const input_property& m_H; + const input_property& m_dH; + + const tensor4 m_IIdev; +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_J2_PLASTICITY_H diff --git a/include/numsim-materials/materials/small_strain_plasticity.h b/include/numsim-materials/materials/small_strain_plasticity.h index 67e1f91..67a24a9 100644 --- a/include/numsim-materials/materials/small_strain_plasticity.h +++ b/include/numsim-materials/materials/small_strain_plasticity.h @@ -211,9 +211,10 @@ class small_strain_plasticity final yield_fn m_yf{}; }; -template -using j2_plasticity = small_strain_plasticity>; +// j2_plasticity moved to materials/j2_plasticity.h as a dedicated class: J2 is +// associative and has no apex, so it used none of this class's generality and +// paid two rank-4 contractions per step for a tangent that has a closed form. +// This class now serves the pressure-dependent models it was written for. /// Drucker-Prager plasticity. The yield function (with η, β, K_bulk) must be /// supplied via the "yield_function" parameter at construction. diff --git a/tests/plot_data.cpp b/tests/plot_data.cpp index 73dc427..35fd728 100644 --- a/tests/plot_data.cpp +++ b/tests/plot_data.cpp @@ -7,6 +7,7 @@ #include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/drucker_prager_yield_function.h" #include "numsim-materials/materials/small_strain_plasticity.h" +#include "numsim-materials/materials/j2_plasticity.h" #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/postprocessing/numerical_diff_checker.h" diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index 9c1c6fc..3fabd3f 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include "numsim-materials/core/material_context.h" @@ -6,6 +7,7 @@ #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/small_strain_plasticity.h" +#include "numsim-materials/materials/j2_plasticity.h" #include "numsim-materials/materials/rk_plasticity.h" #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/solvers/butcher_tableau.h" @@ -192,19 +194,39 @@ class J2TangentTest : public ::testing::Test { }; TEST_F(J2TangentTest, ConsistentTangentAllSteps) { - T max_rel_error = 0; + // The elastic->plastic transition step is genuinely inexact: the tangent is + // discontinuous across the yield surface, so a central difference straddling + // it averages two different tangents. That step gets a loose bound. + // + // Every OTHER step is fully plastic and matches at ~1e-10. Bounding the whole + // run by the transition step's 10% made the check blind: a 0.5% error in the + // closed-form tangent still landed at 2.4e-4 and passed. The two regimes are + // now bounded separately. + T worst_plastic = 0, worst_transition = 0; + T alpha_before = 0; + int plastic_steps = 0; for (int i = 0; i < 20; ++i) { ctx.update(); - auto rel = ctx.get("checker", "rel_error"); - auto alpha = ctx.get("j2", "equivalent_plastic_strain"); - std::println(" step {:2d}: rel={:.2e} alpha={:.4e}", i, rel, alpha); - if (rel > max_rel_error) max_rel_error = rel; + const auto rel = ctx.get("checker", "rel_error"); + const auto alpha = ctx.get("j2", "equivalent_plastic_strain"); + const bool was_plastic = alpha_before > T{1e-10}; + const bool is_plastic = alpha > T{1e-10}; + if (was_plastic && is_plastic) { // fully inside the plastic regime + worst_plastic = std::max(worst_plastic, rel); + ++plastic_steps; + } else if (is_plastic) { // the crossing step + worst_transition = std::max(worst_transition, rel); + } + alpha_before = alpha; ctx.commit(); } - // Transition steps (elastic→plastic) show ~5% error due to yield surface crossing. - // Fully elastic and fully plastic steps match at machine precision. - EXPECT_LT(max_rel_error, 0.1) - << "Consistent tangent should match numerical derivative"; + + ASSERT_GT(plastic_steps, 5) << "the path must spend real time yielding"; + EXPECT_LT(worst_plastic, 1e-8) + << "the consistent tangent must match the numerical derivative on fully " + "plastic steps (worst " << worst_plastic << ")"; + EXPECT_LT(worst_transition, 0.1) + << "the yield-crossing step is inexact but should stay bounded"; } // --- RK plasticity: multi-stage return mapping via tableau parameter --- From 01a2403ea7affdf5335bbb8d83a2b7709105f074 Mon Sep 17 00:00:00 2001 From: petlenz Date: Fri, 4 Sep 2026 22:47:26 +0200 Subject: [PATCH 2/6] materials: collapse the two rank-4 contractions in the consistent tangent compute_tangent formed C_e : (dN/dsigma) : C_e explicitly -- two rank-4 x rank-4 contractions per plastic step, and the dominant cost in both Drucker-Prager and the RK integrators. Every flow normal here is a deviatoric term plus a constant volumetric one, so dN/dsigma is deviatoric in BOTH index pairs. For isotropic C_e that gives C_e : X = 2G X and X : C_e = 2G X, so the whole expression is C_e : (dN/dsigma) : C_e == 4 G^2 (dN/dsigma) Verified against the explicit form at 2.5e-16 for the J2 and Drucker-Prager derivatives independently, before changing any material. Isotropy is not newly assumed: effective_modulus() already required it, as 3G for J2 and G + K*eta*beta for Drucker-Prager. Drucker-Prager: 1178.9 -> 761.3 ns/step (1.55x) Bit-identical output, to all 17 digits, on a 20-step path: before: alpha=0.042067012632805115 s00=108.98966076999213 C0000=251.4917329876391 after : alpha=0.042067012632805115 s00=108.98966076999213 C0000=251.4917329876391 rk_plasticity shares compute_tangent and gets the same reduction; its three tests pass unchanged. The DP tangent is independently checked against a numerical derivative on six load paths including the apex branch, which is what would catch the collapse being wrong rather than merely faster. --- .../materials/plasticity_utils.h | 16 ++++++++++++++-- .../numsim-materials/materials/rk_plasticity.h | 2 +- .../materials/small_strain_plasticity.h | 2 +- 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/include/numsim-materials/materials/plasticity_utils.h b/include/numsim-materials/materials/plasticity_utils.h index 222e9e4..89c2f85 100644 --- a/include/numsim-materials/materials/plasticity_utils.h +++ b/include/numsim-materials/materials/plasticity_utils.h @@ -111,7 +111,8 @@ tmech::tensor compute_tangent( T sig_eq, T total_dlambda, T dH_val, - const tmech::tensor& C_e) + const tmech::tensor& C_e, + T G) { using tensor2 = tmech::tensor; using tensor4 = tmech::tensor; @@ -140,7 +141,18 @@ tmech::tensor compute_tangent( // flow_normal_stress_derivative takes (sig_dev, sig_eq) to avoid // cancellation error from reconstructing s from N. const tensor4 dN_dsig{yf.flow_normal_stress_derivative(sig_dev, sig_eq)}; - const tensor4 C_dN_C{tmech::dcontract(C_e, tmech::dcontract(dN_dsig, C_e))}; + + // C_e : (dN/dsig) : C_e == 4G^2 (dN/dsig). + // + // Every flow normal here is deviatoric plus a constant volumetric part, so + // dN/dsig is deviatoric in BOTH index pairs. For isotropic C_e that makes + // C_e : X = 2G X and X : C_e = 2G X, and the two rank-4 x rank-4 + // contractions collapse to a scalar multiply. Verified to 2.5e-16 against + // the explicit form for both the J2 and Drucker-Prager derivatives. + // + // Isotropy is not newly assumed: effective_modulus() above already requires + // it (3G for J2, G + K*eta*beta for DP). + const tensor4 C_dN_C{T{4} * G * G * dN_dsig}; const tensor4 A{C_e - total_dlambda * C_dN_C}; return A + tmech::otimes(dsig_ddlambda, dlambda_deps); diff --git a/include/numsim-materials/materials/rk_plasticity.h b/include/numsim-materials/materials/rk_plasticity.h index 909c1dc..cf5b0b1 100644 --- a/include/numsim-materials/materials/rk_plasticity.h +++ b/include/numsim-materials/materials/rk_plasticity.h @@ -172,7 +172,7 @@ class rk_plasticity final m_yf, eps, eps_p_new, C_e, m_sigma_0, m_H.get()); m_tangent = plasticity_detail::compute_tangent( m_yf, converged.sig_dev, converged.N, converged.sig_eq, - total_dlambda, m_dH.get(), C_e); + total_dlambda, m_dH.get(), C_e, m_G); } private: diff --git a/include/numsim-materials/materials/small_strain_plasticity.h b/include/numsim-materials/materials/small_strain_plasticity.h index 67a24a9..74af287 100644 --- a/include/numsim-materials/materials/small_strain_plasticity.h +++ b/include/numsim-materials/materials/small_strain_plasticity.h @@ -169,7 +169,7 @@ class small_strain_plasticity final // For J2, trial = converged. For DP, they differ. m_H.update_source(); m_tangent = plasticity_detail::compute_tangent( - 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. From ba8d1c37be258ca5147bb0f0ea07fd2f5a7796f3 Mon Sep 17 00:00:00 2001 From: petlenz Date: Fri, 4 Sep 2026 22:55:15 +0200 Subject: [PATCH 3/6] materials: collapse small_strain_plasticity into a dedicated Drucker-Prager After J2 moved out, small_strain_plasticity had exactly one instantiation. A template parameter with one argument is not generality, it is indirection, and it cost: - a has_apex_return concept plus three if constexpr / requires sites, guarding a branch the only remaining user always has; - a yield function passed as a C++ OBJECT through a "yield_function" parameter. That second one was not just noise. The JSON reader has no converter for the object, so Drucker-Prager could not be configured from a document at all -- the blocker behind #33, where it had to stay unregistered because a document naming it would silently get a default-constructed cone (eta = beta = k = 0), which builds, runs, never yields, and looks like elasticity. eta, beta and K_bulk are now ordinary scalar parameters. Verified: a Drucker-Prager model built entirely from a JSON document reproduces the C++ reference bit-for-bit, alpha = 0.042067012632805115 either way. The apex is now unconditional -- the concept and every if constexpr are gone -- because a cone always has one. Drucker-Prager: 1178.9 -> 721.9 ns/step over the two commits (1.63x) Bit-identical to the pre-refactor implementation on a 20-step path, all 17 digits of alpha, stress and tangent. The yield function survives as an internal member rather than a template parameter: it holds the verified apex algebra, and rewriting that to save a file would have traded a real risk for a cosmetic gain. --- ...asticity.h => drucker_prager_plasticity.h} | 85 +++++++------------ tests/debug_apex.cpp | 3 +- tests/plot_data.cpp | 7 +- tests/test_drucker_prager.cpp | 32 ++++--- tests/test_j2_plasticity.cpp | 2 +- 5 files changed, 55 insertions(+), 74 deletions(-) rename include/numsim-materials/materials/{small_strain_plasticity.h => drucker_prager_plasticity.h} (72%) diff --git a/include/numsim-materials/materials/small_strain_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h similarity index 72% rename from include/numsim-materials/materials/small_strain_plasticity.h rename to include/numsim-materials/materials/drucker_prager_plasticity.h index 74af287..7797fd8 100644 --- a/include/numsim-materials/materials/small_strain_plasticity.h +++ b/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 #include @@ -8,47 +8,34 @@ #include #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 -concept has_apex_return = requires(const YF& yf, - const tmech::tensor& t2, T v) { - { yf.needs_apex_return(v, v, v) } -> std::convertible_to; - { yf.apex_modified_sig_eq(t2) } -> std::convertible_to; - { yf.apex_effective_modulus() } -> std::convertible_to; - { 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 -class small_strain_plasticity final - : public material_base, Traits> { +template +class drucker_prager_plasticity final + : public material_base, Traits> { public: - using base = material_base, Traits>; + using base = material_base, 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; using tensor4 = tmech::tensor; - using yield_fn = YieldFunction; + using yield_fn = drucker_prager_yield_function; using solver_type = backward_euler; template - explicit small_strain_plasticity(Args&&... args) + explicit drucker_prager_plasticity(Args&&... args) : base(std::forward(args)...), m_stress(base::template add_output( - "stress", &small_strain_plasticity::compute)), + "stress", &drucker_prager_plasticity::compute)), m_tangent(base::template add_output("tangent")), m_eps_p(base::template add_history_output("plastic_strain")), m_kappa(base::template add_history_output("equivalent_plastic_strain")), @@ -69,8 +56,9 @@ class small_strain_plasticity final base::template get_parameter("hardening_source"), "hardening_modulus", EdgeKind::Local)) { - if (base::m_parameter_handler.contains("yield_function")) - m_yf = base::template get_parameter("yield_function"); + m_yf = yield_fn(base::template get_parameter("eta"), + base::template get_parameter("beta"), + base::template get_parameter("K_bulk")); } static input_parameter_controller parameters() { @@ -81,6 +69,13 @@ class small_strain_plasticity final para.template insert("solver_source").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); + // 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("eta").template add(); + para.template insert("beta").template add(); + para.template insert("K_bulk").template add(); return para; } @@ -104,28 +99,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) { - 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) { - 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); @@ -174,11 +162,8 @@ class small_strain_plasticity final /// 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 - { + 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); @@ -211,17 +196,7 @@ class small_strain_plasticity final yield_fn m_yf{}; }; -// j2_plasticity moved to materials/j2_plasticity.h as a dedicated class: J2 is -// associative and has no apex, so it used none of this class's generality and -// paid two rank-4 contractions per step for a tangent that has a closed form. -// This class now serves the pressure-dependent models it was written for. - -/// Drucker-Prager plasticity. The yield function (with η, β, K_bulk) must be -/// supplied via the "yield_function" parameter at construction. -template -using drucker_prager_plasticity = small_strain_plasticity>; } // namespace numsim::materials -#endif // NUMSIM_MATERIALS_SMALL_STRAIN_PLASTICITY_H +#endif // NUMSIM_MATERIALS_DRUCKER_PRAGER_PLASTICITY_H diff --git a/tests/debug_apex.cpp b/tests/debug_apex.cpp index 2dce047..52e84d2 100644 --- a/tests/debug_apex.cpp +++ b/tests/debug_apex.cpp @@ -5,7 +5,7 @@ #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/drucker_prager_yield_function.h" -#include "numsim-materials/materials/small_strain_plasticity.h" +#include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/solvers/backward_euler.h" using policy = numsim::materials::material_policy_default; @@ -22,7 +22,6 @@ int main() { const T lambda = K - T{2}*G/T{3}; // 115385 const T G_eff = G + K*eta*beta; // 84423 - dp_yield yf(eta, beta, K); std::println("Elastic constants: lambda={:.1f}, G={:.1f}, K={:.1f}", lambda, G, K); std::println("G_eff = {:.1f}", G_eff); diff --git a/tests/plot_data.cpp b/tests/plot_data.cpp index 35fd728..b32ed3f 100644 --- a/tests/plot_data.cpp +++ b/tests/plot_data.cpp @@ -6,7 +6,7 @@ #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/drucker_prager_yield_function.h" -#include "numsim-materials/materials/small_strain_plasticity.h" +#include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/materials/j2_plasticity.h" #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/postprocessing/numerical_diff_checker.h" @@ -141,7 +141,6 @@ run_result run_dp(T increment, int steps, p.insert("K", H_mod); ctx.create>(p); - dp_yield yf(T{0.3}, T{0.15}, K_val); p.clear(); p.insert("name", "dp"); @@ -151,7 +150,9 @@ run_result run_dp(T increment, int steps, p.insert("solver_source", "solver"); p.insert("G", G_val); p.insert("sigma_0", sigma_0); - p.insert("yield_function", yf); + p.insert("eta", T{0.3}); + p.insert("beta", T{0.15}); + p.insert("K_bulk", K_val); ctx.create(p); p.clear(); diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp index ebe7c7d..aa8f4fc 100644 --- a/tests/test_drucker_prager.cpp +++ b/tests/test_drucker_prager.cpp @@ -8,7 +8,7 @@ #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_isotropic_hardening.h" #include "numsim-materials/materials/drucker_prager_yield_function.h" -#include "numsim-materials/materials/small_strain_plasticity.h" +#include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/materials/rk_plasticity.h" #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/solvers/butcher_tableau.h" @@ -113,7 +113,6 @@ class DruckerPragerTest : public ::testing::Test { ctx.create>(p); // Drucker-Prager yield function with friction and dilatancy - dp_yield yf(dp_eta, dp_beta, K); p.clear(); p.insert("name", "dp"); @@ -123,7 +122,9 @@ class DruckerPragerTest : public ::testing::Test { p.insert("solver_source", "solver"); p.insert("G", G); p.insert("sigma_0", cohesion); - p.insert("yield_function", yf); + p.insert("eta", dp_eta); + p.insert("beta", dp_beta); + p.insert("K_bulk", K); ctx.create(p); ctx.finalize(); @@ -217,7 +218,6 @@ class DPTangentTest : public ::testing::Test { p.insert("K", T{500.0}); ctx.create>(p); - dp_yield yf(T{0.1}, T{0.05}, T{166.67}); p.clear(); p.insert("name", "dp"); @@ -227,7 +227,9 @@ class DPTangentTest : public ::testing::Test { p.insert("solver_source", "solver"); p.insert("G", T{76.92}); p.insert("sigma_0", T{20.0}); - p.insert("yield_function", yf); + p.insert("eta", T{0.1}); + p.insert("beta", T{0.05}); + p.insert("K_bulk", T{166.67}); ctx.create(p); p.clear(); @@ -290,7 +292,6 @@ T run_dp_max_tangent_error(T increment, int steps) { p.insert("K", T{500.0}); ctx.create>(p); - dp_yield yf(T{0.1}, T{0.05}, T{166.67}); p.clear(); p.insert("name", "dp"); @@ -300,7 +301,9 @@ T run_dp_max_tangent_error(T increment, int steps) { p.insert("solver_source", "solver"); p.insert("G", T{76.92}); p.insert("sigma_0", T{20.0}); - p.insert("yield_function", yf); + p.insert("eta", T{0.1}); + p.insert("beta", T{0.05}); + p.insert("K_bulk", T{166.67}); ctx.create(p); p.clear(); @@ -378,7 +381,6 @@ T max_tangent_error(std::vector direction, T increment, int steps) { p.insert("K", T{500.0}); ctx.create>(p); - dp_yield yf(T{0.1}, T{0.05}, T{166.67}); p.clear(); p.insert("name", "dp"); p.insert("elastic_source", "elastic"); @@ -387,7 +389,9 @@ T max_tangent_error(std::vector direction, T increment, int steps) { p.insert("solver_source", "solver"); p.insert("G", T{76.92}); p.insert("sigma_0", T{20.0}); - p.insert("yield_function", yf); + p.insert("eta", T{0.1}); + p.insert("beta", T{0.05}); + p.insert("K_bulk", T{166.67}); ctx.create(p); p.clear(); @@ -475,7 +479,6 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) { p.insert("K", H_mod); ctx.create>(p); - dp_yield yf(dp_eta, dp_beta, K); p.clear(); p.insert("name", "dp"); p.insert("elastic_source", "elastic"); @@ -484,7 +487,9 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) { p.insert("solver_source", "solver"); p.insert("G", G); p.insert("sigma_0", cohesion); - p.insert("yield_function", yf); + p.insert("eta", dp_eta); + p.insert("beta", dp_beta); + p.insert("K_bulk", K); ctx.create(p); ctx.finalize(); @@ -533,7 +538,6 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) { p.insert("source", "dp"); p.insert("K", H_mod); ctx.create>(p); - dp_yield yf(dp_eta, dp_beta, K); p.clear(); p.insert("name", "dp"); p.insert("elastic_source", "elastic"); @@ -541,7 +545,9 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) { p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); p.insert("G", G); p.insert("sigma_0", cohesion); - p.insert("yield_function", yf); + p.insert("eta", dp_eta); + p.insert("beta", dp_beta); + p.insert("K_bulk", K); ctx.create(p); ctx.finalize(); diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index 3fabd3f..03ff086 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -6,7 +6,7 @@ #include "numsim-materials/materials/tensor_component_stepper.h" #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_isotropic_hardening.h" -#include "numsim-materials/materials/small_strain_plasticity.h" +#include "numsim-materials/materials/drucker_prager_plasticity.h" #include "numsim-materials/materials/j2_plasticity.h" #include "numsim-materials/materials/rk_plasticity.h" #include "numsim-materials/solvers/backward_euler.h" From 7ff63457f3bdd231efdfb14e923b0a18a156b147 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 10:38:00 +0200 Subject: [PATCH 4/6] materials: drop the redundant stress contraction in the return mapping Profiling the step showed one rank-4 : rank-2 contraction at ~33 ns against ~245 ns for a whole J2 step, and the return map did two of them: once for the trial stress, once for the returned stress. The second is redundant: sigma = C_e : (eps - eps_p_old - dl N) = sig_trial - dl (C_e : N) and for isotropic C_e, C_e : N = 2G dev(N) + K tr(N) I -- the same identity the tangent collapse already uses. J2's flow is deviatoric so tr(N) = 0 and it reduces to 2G N; Drucker-Prager's is not, so the volumetric term stays. J2 (same-process A/B): 245.2 -> 195.4 ns/step DP (interleaved A/B): ~8%, medians 1037 -> 953 ns/step NOT bit-identical this time, unlike the earlier steps: the result agrees to about 1 ULP (s00 108.7500531151484 vs ...843) because the arithmetic is algebraically equal but ordered differently. All 50 tests pass, including the FD tangent checks on six load paths. The apex return keeps its contraction: there eps_p comes from apex_plastic_strain rather than dl*N, so the identity does not apply, and the branch is rare. Measurement note: absolute ns figures in this branch's earlier commits were single runs taken at different times, and this machine's load moves them by 30%+. The J2 comparison above is a same-process A/B; the DP one is interleaved across alternating runs. Ratios are meaningful, absolutes are not comparable across commits. --- .../materials/drucker_prager_plasticity.h | 13 ++++++++++++- include/numsim-materials/materials/j2_plasticity.h | 9 +++++++-- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h index 7797fd8..89ef2a8 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/include/numsim-materials/materials/drucker_prager_plasticity.h @@ -41,6 +41,7 @@ class drucker_prager_plasticity final m_kappa(base::template add_history_output("equivalent_plastic_strain")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), + m_K_bulk(base::template get_parameter("K_bulk")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), m_C_e(base::template add_input( @@ -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(); + 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,6 +197,7 @@ class drucker_prager_plasticity final const value_type& m_G; const value_type& m_sigma_0; + const value_type& m_K_bulk; material_ref& m_solver; const input_property& m_C_e; diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h index 7322fd4..07b352e 100644 --- a/include/numsim-materials/materials/j2_plasticity.h +++ b/include/numsim-materials/materials/j2_plasticity.h @@ -143,8 +143,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; From 32f573c9377d372008eab2d4c0e0107088bf8b32 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 11:07:57 +0200 Subject: [PATCH 5/6] materials: plasticity builds its own elastic tangent Both plasticity classes read the elastic stiffness from an elastic_source's "tangent" property. That looked like flexibility and was not: the closed forms REQUIRE an isotropic C_e -- it is what makes C_e : N = 2G N and N : C_e : N = 3G true -- so accepting an arbitrary rank-4 tangent advertised a generality neither class can honour. It also dragged a linear_elasticity into every plasticity graph, and that material's own "stress" output is not merely unused there, it is WRONG: it is C : eps with eps_p ignored, so once yielding starts it is not the stress of anything. Measured on a 200-step path it over-predicts by a growing margin: step 40 alpha 0.0094 j2 106.25 elastic 107.69 +1.4% step 120 alpha 0.1094 j2 306.25 elastic 323.08 +5.5% step 200 alpha 0.2094 j2 506.25 elastic 538.46 +6.4% A postprocessor logging elastic::stress from a plasticity graph gets that, under a name that reads as authoritative. Both classes now build C_e from moduli they already hold or now take: j2_plasticity gains a "K" parameter; drucker_prager_plasticity needs nothing new, since K_bulk and G were already required for the cone. elastic_source is gone from both, and a plasticity graph no longer contains an elastic material at all. J2: 195.4 -> 168.1 ns/step DP: ~720 -> ~608 ns/step Values agree with the pre-refactor implementation to about 1 ULP on alpha, stress and tangent. linear_elasticity is untouched and keeps its users: for a genuinely elastic model C : eps IS the answer, and isotropic_damage consumes both its stress and its tangent legitimately. rk_plasticity still takes an elastic_source; it is the remaining templated class and is left alone here. An earlier draft of this reasoning proposed dead-property elimination in the property engine, with a declared-outputs mechanism to tell an unread property from one a host reads through ctx.get(). That was solving the symptom. The property should not be in the graph. --- .../materials/drucker_prager_plasticity.h | 20 ++++++---- .../materials/j2_plasticity.h | 37 +++++++++++++------ tests/plot_data.cpp | 3 +- tests/test_drucker_prager.cpp | 6 --- tests/test_j2_plasticity.cpp | 8 ++-- 5 files changed, 44 insertions(+), 30 deletions(-) diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h index 89ef2a8..67aed65 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/include/numsim-materials/materials/drucker_prager_plasticity.h @@ -44,9 +44,6 @@ class drucker_prager_plasticity final m_K_bulk(base::template get_parameter("K_bulk")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), - m_C_e(base::template add_input( - base::template get_parameter("elastic_source"), - "tangent", EdgeKind::Global)), m_strain(base::template add_input( base::template get_parameter("strain_source"), "strain", EdgeKind::Global)), @@ -59,12 +56,15 @@ class drucker_prager_plasticity final { m_yf = yield_fn(base::template get_parameter("eta"), base::template get_parameter("beta"), - base::template get_parameter("K_bulk")); + m_K_bulk); + const auto I = tmech::eye(); + 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(); } static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("elastic_source").template add(); para.template insert("hardening_source").template add(); para.template insert("strain_source").template add(); para.template insert("solver_source").template add(); @@ -81,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; @@ -200,11 +200,17 @@ class drucker_prager_plasticity final const value_type& m_K_bulk; material_ref& m_solver; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& 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{}; }; diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h index 07b352e..dee151c 100644 --- a/include/numsim-materials/materials/j2_plasticity.h +++ b/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 class j2_plasticity final : public material_base, Traits> { @@ -71,11 +72,9 @@ class j2_plasticity final "equivalent_plastic_strain")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), + m_K(base::template get_parameter("K")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), - m_C_e(base::template add_input( - base::template get_parameter("elastic_source"), - "tangent", EdgeKind::Global)), m_strain(base::template add_input( base::template get_parameter("strain_source"), "strain", EdgeKind::Global)), @@ -85,26 +84,41 @@ class j2_plasticity final m_dH(base::template add_input( base::template get_parameter("hardening_source"), "hardening_modulus", EdgeKind::Local)), - m_IIdev(plasticity_detail::make_IIdev()) + m_IIdev(plasticity_detail::make_IIdev()), + 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(); + 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("elastic_source") - .template add(); para.template insert("hardening_source") .template add(); para.template insert("strain_source") .template add(); para.template insert("solver_source") .template add(); + para.template insert("K").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); 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; @@ -167,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& m_solver; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& m_dH; const tensor4 m_IIdev; + const tensor4 m_C_e; }; } // namespace numsim::materials diff --git a/tests/plot_data.cpp b/tests/plot_data.cpp index b32ed3f..b32db4b 100644 --- a/tests/plot_data.cpp +++ b/tests/plot_data.cpp @@ -78,10 +78,10 @@ run_result run_j2(T increment, int steps, p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", K_val); p.insert("G", G_val); p.insert("sigma_0", sigma_0); ctx.create(p); @@ -144,7 +144,6 @@ run_result run_dp(T increment, int steps, p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp index aa8f4fc..c495c32 100644 --- a/tests/test_drucker_prager.cpp +++ b/tests/test_drucker_prager.cpp @@ -116,7 +116,6 @@ class DruckerPragerTest : public ::testing::Test { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -221,7 +220,6 @@ class DPTangentTest : public ::testing::Test { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -295,7 +293,6 @@ T run_dp_max_tangent_error(T increment, int steps) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -383,7 +380,6 @@ T max_tangent_error(std::vector direction, T increment, int steps) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -481,7 +477,6 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -540,7 +535,6 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) { ctx.create>(p); p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index 03ff086..bc88a7a 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -56,10 +56,10 @@ class J2PlasticityTest : public ::testing::Test { // J2 plasticity — solver passed as pointer p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); - p.insert("hardening_source", "hardening"); + p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", K); p.insert("G", G); p.insert("sigma_0", sigma_0); ctx.create>(p); @@ -168,10 +168,10 @@ class J2TangentTest : public ::testing::Test { p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); - p.insert("hardening_source", "hardening"); + p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", T{166.67}); p.insert("G", T{76.92}); p.insert("sigma_0", T{50.0}); ctx.create>(p); From e428ccbeb7f9eb5633377d15a33bdfed9ed6e500 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sat, 5 Sep 2026 11:13:58 +0200 Subject: [PATCH 6/6] materials: rk_plasticity builds its own elastic tangent too Same change as the other two plasticity classes, for the same reason. rk_plasticity read the stiffness from an elastic_source while assuming isotropy twice over: effective_modulus(G) = 3G in its stage residuals, and compute_tangent's C_e : X = 2G X collapse. An arbitrary rank-4 tangent could not have been honoured by either. Takes "K" and builds C_e in the constructor; elastic_source is gone. Bit- identical over a 40-step SDIRK3 path, all 17 digits: before: 0.01062302967272258 108.7500531151484 250.00062497968736 76.554019397938148 after : 0.01062302967272258 108.7500531151484 250.00062497968736 76.554019397938148 That comparison matters here because the suite's rk-vs-j2 equivalence test now has both sides changed; this one holds rk against its own pre-change output. SDIRK3TangentCheck, which compares the tangent against a numerical derivative, is the independent check and passes unchanged. No plasticity material takes an elastic_source now. The one remaining consumer is isotropic_damage, which reads elastic::stress legitimately -- damage scales an elastic stress, and with no plastic strain in that model C : eps IS the stress. --- .../materials/rk_plasticity.h | 23 +++++++++++++------ tests/test_j2_plasticity.cpp | 2 +- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/include/numsim-materials/materials/rk_plasticity.h b/include/numsim-materials/materials/rk_plasticity.h index cf5b0b1..05d40bd 100644 --- a/include/numsim-materials/materials/rk_plasticity.h +++ b/include/numsim-materials/materials/rk_plasticity.h @@ -37,16 +37,14 @@ class rk_plasticity final m_tangent(base::template add_output("tangent")), m_eps_p(base::template add_history_output("plastic_strain")), m_kappa(base::template add_history_output("equivalent_plastic_strain")), + m_K(base::template get_parameter("K")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), m_tol(base::template get_parameter("tolerance")), m_max_iter(base::template get_parameter("max_iter")), m_tableau(base::template get_parameter("tableau")), - m_elastic_source(base::template get_parameter("elastic_source")), m_hardening_source(base::template get_parameter("hardening_source")), m_strain_source(base::template get_parameter("strain_source")), - m_C_e(base::template add_input( - m_elastic_source, "tangent", EdgeKind::Global)), m_strain(base::template add_input( m_strain_source, "strain", EdgeKind::Global)), m_H(base::template add_input( @@ -54,6 +52,11 @@ class rk_plasticity final m_dH(base::template add_input( m_hardening_source, "hardening_modulus", EdgeKind::Local)) { + const auto I = tmech::eye(); + 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(); + 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("elastic_source").template add(); para.template insert("hardening_source").template add(); para.template insert("strain_source").template add(); + para.template insert("K").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); para.template insert("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& m_eps_p; history_property& 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& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& 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 m_dlambda; std::vector m_N_stage; diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index bc88a7a..a7b3429 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -258,9 +258,9 @@ class RKPlasticityTest : public ::testing::Test { p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); + p.insert("K", T{166.67}); p.insert("G", T{76.92}); p.insert("sigma_0", T{50.0}); p.insert("tableau", &m_tab);