From 24643bb7a6ddc3f31218693bccc604cc26122106 Mon Sep 17 00:00:00 2001 From: petlenz Date: Fri, 4 Sep 2026 22:55:15 +0200 Subject: [PATCH 1/2] 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 2c94a8567a6600ed8eaeca9c01fe54846fb3d428 Mon Sep 17 00:00:00 2001 From: petlenz Date: Sun, 6 Sep 2026 20:43:19 +0200 Subject: [PATCH 2/2] headers: name each include guard after its own file The 1 header(s) this branch introduces follow the convention set on feature/drucker-prager: the guard is the file's own name, no NUMSIM_MATERIALS_ prefix. Checked against every dependency header and /usr/include for a prior #define of each new name -- none. --- .../numsim-materials/materials/drucker_prager_plasticity.h | 6 +++--- 1 file changed, 3 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..1e73725 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/include/numsim-materials/materials/drucker_prager_plasticity.h @@ -1,5 +1,5 @@ -#ifndef NUMSIM_MATERIALS_DRUCKER_PRAGER_PLASTICITY_H -#define NUMSIM_MATERIALS_DRUCKER_PRAGER_PLASTICITY_H +#ifndef DRUCKER_PRAGER_PLASTICITY_H +#define DRUCKER_PRAGER_PLASTICITY_H #include #include @@ -199,4 +199,4 @@ class drucker_prager_plasticity final } // namespace numsim::materials -#endif // NUMSIM_MATERIALS_DRUCKER_PRAGER_PLASTICITY_H +#endif // DRUCKER_PRAGER_PLASTICITY_H