From 77a9d0bb786fc8004c8b8892ab82053590b497f5 Mon Sep 17 00:00:00 2001 From: Ingo Scheider Date: Wed, 29 Apr 2026 11:24:10 +0200 Subject: [PATCH] Refactor the drucker prager model to use tensors --- ...4C_global_legacy_module_validmaterials.cpp | 26 +- src/mat/4C_mat_plasticdruckerprager.cpp | 381 ++++++------------ src/mat/4C_mat_plasticdruckerprager.hpp | 73 +--- unittests/mat/4C_druckerprager_test.cpp | 259 ++++++------ 4 files changed, 289 insertions(+), 450 deletions(-) diff --git a/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp b/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp index aa4154d51ea..4925961a74a 100644 --- a/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp +++ b/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp @@ -1205,22 +1205,30 @@ std::unordered_map Global::v /*----------------------------------------------------------------------*/ // Plastic linear elastic St.Venant Kirchhoff / Drucker Prager plasticity { + using namespace Core::IO::InputSpecBuilders::Validators; known_materials[Core::Materials::m_pldruckprag] = group("MAT_Struct_DruckerPrager", { - parameter("YOUNG", {.description = "Young's modulus"}), - parameter("NUE", {.description = "Poisson's ratio"}), - parameter("DENS", {.description = "Density"}), - parameter("ISOHARD", {.description = "linear isotropic hardening"}), - parameter("TOL", {.description = "Local Newton iteration tolerance"}), - parameter("C", {.description = "cohesion"}), + parameter( + "YOUNG", {.description = "Young's modulus", .validator = positive()}), + parameter( + "NUE", {.description = "Poisson's ratio (must be in (-1, 0.5) for stability)", + .validator = in_range(excl(-1.0), excl(0.5))}), + parameter( + "DENS", {.description = "Density", .validator = positive_or_zero()}), + parameter("ISOHARD", {.description = "linear isotropic hardening", + .validator = positive_or_zero()}), + parameter("TOL", {.description = "Local Newton iteration tolerance", + .validator = positive()}), + parameter("C", {.description = "cohesion (must be strictly positive)", + .validator = positive()}), parameter("ETA", {.description = "Drucker Prager Constant Eta"}), parameter("XI", {.description = "Drucker Prager Constant Xi"}), parameter("ETABAR", {.description = "Drucker Prager Constant Etabar"}), parameter("TANG", {.description = "Method to compute the material tangent", .default_value = "consistent"}), - parameter( - "MAXITER", {.description = "Maximum Iterations for local Neutron Raphson", - .default_value = 50}), + parameter("MAXITER", {.description = "Maximum Iterations for local Newton Raphson", + .default_value = 50, + .validator = positive()}), }, {.description = "elastic St.Venant Kirchhoff / plastic drucker prager"}); } diff --git a/src/mat/4C_mat_plasticdruckerprager.cpp b/src/mat/4C_mat_plasticdruckerprager.cpp index ce163c41db2..b1e0daa6058 100644 --- a/src/mat/4C_mat_plasticdruckerprager.cpp +++ b/src/mat/4C_mat_plasticdruckerprager.cpp @@ -11,13 +11,10 @@ #include "4C_comm_parobject.hpp" #include "4C_global_data.hpp" #include "4C_inpar_structure.hpp" -#include "4C_linalg_FADmatrix_utils.hpp" #include "4C_linalg_fixedsizematrix.hpp" -#include "4C_linalg_fixedsizematrix_voigt_notation.hpp" #include "4C_linalg_tensor.hpp" #include "4C_linalg_tensor_conversion.hpp" #include "4C_linalg_tensor_generators.hpp" -#include "4C_mat_material_factory.hpp" #include "4C_mat_par_bundle.hpp" #include "4C_mat_stvenantkirchhoff.hpp" #include "4C_material_base.hpp" @@ -115,24 +112,24 @@ void Mat::PlasticDruckerPrager::unpack(Core::Communication::UnpackBuffer& buffer if (histsize == 0) isinit_ = false; - strainpllast_ = std::vector>(); - strainplcurr_ = std::vector>(); + strainpllast_ = std::vector>(); + strainplcurr_ = std::vector>(); strainbarpllast_ = std::vector(); strainbarplcurr_ = std::vector(); for (int var = 0; var < histsize; ++var) { - Core::LinAlg::Matrix tmp_vector(Core::LinAlg::Initialization::zero); + Core::LinAlg::SymmetricTensor tmp_tensor{}; double tmp_scalar = 0.0; // last converged states are unpacked - extract_from_pack(buffer, tmp_vector); - strainpllast_.push_back(tmp_vector); + extract_from_pack(buffer, tmp_tensor); + strainpllast_.push_back(tmp_tensor); extract_from_pack(buffer, tmp_scalar); strainbarpllast_.push_back(tmp_scalar); // current iteration states are unpacked - extract_from_pack(buffer, tmp_vector); - strainplcurr_.push_back(tmp_vector); + extract_from_pack(buffer, tmp_tensor); + strainplcurr_.push_back(tmp_tensor); extract_from_pack(buffer, tmp_scalar); strainbarplcurr_.push_back(tmp_scalar); } @@ -141,11 +138,11 @@ void Mat::PlasticDruckerPrager::unpack(Core::Communication::UnpackBuffer& buffer void Mat::PlasticDruckerPrager::setup(int numgp, const Discret::Elements::Fibers& fibers, const std::optional& coord_system) { - strainpllast_.resize(numgp); - strainplcurr_.resize(numgp); + strainpllast_ = std::vector>(numgp); + strainplcurr_ = std::vector>(numgp); - strainbarpllast_.resize(numgp); - strainbarplcurr_.resize(numgp); + strainbarpllast_.assign(numgp, 0.0); + strainbarplcurr_.assign(numgp, 0.0); isinit_ = true; } @@ -154,9 +151,6 @@ void Mat::PlasticDruckerPrager::update() { strainpllast_ = strainplcurr_; strainbarpllast_ = strainbarplcurr_; - - std::for_each(strainplcurr_.begin(), strainplcurr_.end(), [](auto& item) { item.clear(); }); - std::fill(strainbarplcurr_.begin(), strainbarplcurr_.end(), 0.0); } void Mat::PlasticDruckerPrager::setup_cmat( @@ -171,220 +165,127 @@ void Mat::PlasticDruckerPrager::setup_cmat( void Mat::PlasticDruckerPrager::setup_cmat_elasto_plastic_cone( Core::LinAlg::SymmetricTensor& cmat, double Dgamma, double G, double Kappa, - Core::LinAlg::Matrix& devstrain, double xi, double Hiso, double eta, - double etabar) const + const Core::LinAlg::SymmetricTensor& devstrain, double xi, double Hiso, + double eta, double etabar) const { - Core::LinAlg::Matrix cmat_view = - Core::LinAlg::make_stress_like_voigt_view(cmat); - cmat_view.clear(); - - Core::LinAlg::Matrix id2(Core::LinAlg::Initialization::zero); - Core::LinAlg::Voigt::identity_matrix(id2); - Core::LinAlg::Matrix id4sharp(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 3; i++) id4sharp(i, i) = 1.0; - for (int i = 3; i < 6; i++) id4sharp(i, i) = 0.5; - - const double normdevstrain = - sqrt(devstrain(0) * devstrain(0) + devstrain(1) * devstrain(1) + devstrain(2) * devstrain(2) + - 2 * (devstrain(3) * devstrain(3) + devstrain(4) * devstrain(4) + - devstrain(5) * devstrain(5))); - const double epfac = 2 * G * (1 - (Dgamma / sqrt(2) / normdevstrain)); - - cmat_view.update(epfac, id4sharp, 1.0); - - double epfac1 = 0.0; - double epfac2 = 0.0; - double epfac3 = 0.0; - double epfac4 = 0.0; - epfac1 = epfac / (-3.0); - cmat_view.multiply_nt(epfac1, id2, id2, 1.0); - - double A = 0.0; - A = 1 / (G + Kappa * etabar * eta + xi * xi * Hiso); - Core::LinAlg::Matrix D(Core::LinAlg::Initialization::zero); - D.update(1 / normdevstrain, devstrain); - epfac2 = 2 * G * (Dgamma / (sqrt(2) * normdevstrain) - G * A); - - for (int k = 0; k < 6; ++k) - { - for (int i = 0; i < 6; ++i) - { - cmat_view(i, k) += epfac2 * D(i) * D(k); - } - } + constexpr auto id2 = Core::LinAlg::TensorGenerators::identity; + constexpr auto Is = Core::LinAlg::TensorGenerators::symmetric_identity; + constexpr auto Id = Is - 1.0 / 3.0 * dyadic(id2, id2); - epfac3 = -sqrt(2) * G * A * Kappa; - - for (int k = 0; k < 6; ++k) - { - for (int i = 0; i < 6; ++i) - { - cmat_view(k, i) += epfac3 * (eta * D(k) * id2(i) + etabar * id2(k) * D(i)); - } - } + const double nd = std::sqrt(ddot(devstrain, devstrain)); + const auto D = devstrain / nd; - epfac4 = Kappa * (1 - Kappa * eta * etabar * A); + const double A = 1.0 / (G + Kappa * etabar * eta + xi * xi * Hiso); + const double epfac = 2.0 * G * (1.0 - Dgamma / (std::sqrt(2.0) * nd)); + const double epfac2 = 2.0 * G * (Dgamma / (std::sqrt(2.0) * nd) - G * A); + const double epfac3 = -std::sqrt(2.0) * G * A * Kappa; + const double epfac4 = Kappa * (1.0 - Kappa * eta * etabar * A); - for (int k = 0; k < 6; ++k) - { - for (int i = 0; i < 6; ++i) - { - cmat_view(i, k) += epfac4 * id2(k) * id2(i); - } - } + cmat = epfac * Id + epfac2 * dyadic(D, D) + + epfac3 * (eta * dyadic(D, id2) + etabar * dyadic(id2, D)) + epfac4 * dyadic(id2, id2); } void Mat::PlasticDruckerPrager::setup_cmat_elasto_plastic_apex( - Core::LinAlg::SymmetricTensor& cmat, double Kappa, - Core::LinAlg::Matrix& devstrain, double xi, double Hiso, double eta, - double etabar) const + Core::LinAlg::SymmetricTensor& cmat, double Kappa, double xi, double Hiso, + double eta, double etabar) const { - Core::LinAlg::Matrix cmat_view = - Core::LinAlg::make_stress_like_voigt_view(cmat); - - Core::LinAlg::Matrix id2(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 3; i++) id2(i) = 1.0; - double epfac = 0.0; - epfac = Kappa * (1 - Kappa / (Kappa + xi / eta * xi / etabar * Hiso)); - cmat_view.clear(); - cmat_view.multiply_nt(epfac, id2, id2, 0.0); + constexpr auto id2 = Core::LinAlg::TensorGenerators::identity; + const double epfac = Kappa * (1.0 - Kappa / (Kappa + xi / eta * xi / etabar * Hiso)); + cmat = epfac * dyadic(id2, id2); } -template -void Mat::PlasticDruckerPrager::evaluate_fad(const Core::LinAlg::Matrix<3, 3>* defgrd, - const Core::LinAlg::Matrix<6, 1, ScalarT>& linstrain, const Teuchos::ParameterList& params, - Core::LinAlg::Matrix<6, 1, ScalarT>& stress, - Core::LinAlg::SymmetricTensor& cmat, const int gp, const int eleGID) +void Mat::PlasticDruckerPrager::evaluate(const Core::LinAlg::Tensor* defgrad, + const Core::LinAlg::SymmetricTensor& glstrain, + const Teuchos::ParameterList& params, const EvaluationContext<3>& context, + Core::LinAlg::SymmetricTensor& stress, + Core::LinAlg::SymmetricTensor& cmat, int gp, int eleGID) { - Core::LinAlg::Matrix plstrain(Core::LinAlg::Initialization::zero); - - ScalarT young = params_->youngs_; - ScalarT nu = params_->poissonratio_; - ScalarT Hiso = params_->isohard_; - ScalarT cohesion = params_->cohesion_; - ScalarT eta = params_->eta_; - ScalarT xi = params_->xi_; - ScalarT etabar = params_->etabar_; + const double young = params_->youngs_; + const double nu = params_->poissonratio_; + const double Hiso = params_->isohard_; + const double cohesion = params_->cohesion_; + const double eta = params_->eta_; + const double xi = params_->xi_; + const double etabar = params_->etabar_; const int itermax = params_->itermax_; - const std::string tang_str = params_->tang_; - ScalarT G = 0.0; - G = young / (2.0 * (1.0 + nu)); - ScalarT kappa = 0.0; - kappa = young / (3.0 * (1.0 - 2.0 * nu)); - Core::LinAlg::Matrix id2(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 3; i++) id2(i) = 1.0; - Core::LinAlg::Matrix strain(linstrain); - const int tang = std::invoke( - [tang_str]() -> int - { - if (tang_str == "consistent") - return 1; - else if (tang_str == "elastic") - return 0; - else - return 1; - }); - - Core::LinAlg::Matrix strain_p( - Core::LinAlg::Initialization::uninitialized); - for (int i = 0; i < 6; i++) strain_p(i, 0) = strainpllast_.at(gp)(i, 0); - ScalarT strainbar_p = 0.0; - strainbar_p = (strainbarpllast_.at(gp)); - if (strainbarpllast_.at(gp) < 0.0) - FOUR_C_THROW("accumulated plastic strain has to be equal to or greater than zero!"); - - for (int i = 3; i < 6; ++i) strain(i) /= 2.0; - for (int i = 3; i < 6; ++i) strain_p(i) /= 2.0; - - Core::LinAlg::Matrix strain_e(Core::LinAlg::Initialization::zero); - Core::LinAlg::Matrix trialstrain_e( - Core::LinAlg::Initialization::uninitialized); - trialstrain_e.update(1.0, strain, (-1.0), strain_p); - ScalarT tracestrain = trialstrain_e(0) + trialstrain_e(1) + trialstrain_e(2); - Core::LinAlg::Matrix volumetricstrain( - Core::LinAlg::Initialization::uninitialized); - Core::LinAlg::Matrix id2Scalar(Core::LinAlg::Initialization::zero); - for (int i = 0; i < NUM_STRESS_3D; ++i) id2Scalar(i) = static_cast(id2(i)); - volumetricstrain.update((tracestrain / 3.0), id2Scalar); - Core::LinAlg::Matrix devstrain( - Core::LinAlg::Initialization::uninitialized); - devstrain.update(1.0, trialstrain_e, (-1.0), volumetricstrain); - - ScalarT p = kappa * tracestrain; - ScalarT p_trial = p; - Core::LinAlg::Matrix devstress( - Core::LinAlg::Initialization::uninitialized); - devstress.update((2.0 * G), devstrain); - - ScalarT J2 = 1.0 / 2.0 * - (devstress(0) * devstress(0) + devstress(1) * devstress(1) + - devstress(2) * devstress(2)) + - devstress(3) * devstress(3) + devstress(4) * devstress(4) + - devstress(5) * devstress(5); - ScalarT Phi_trial = 0.0; - Phi_trial = std::sqrt(J2) + eta * p - xi * cohesion - xi * Hiso * strainbar_p; - ScalarT Dgamma = 0.0; - ScalarT dstrainv = 0.0; - if (Phi_trial / abs(cohesion) > params_->abstol_) + const std::string& tang_str = params_->tang_; + + const double G = young / (2.0 * (1.0 + nu)); + const double kappa = young / (3.0 * (1.0 - 2.0 * nu)); + constexpr auto id2 = Core::LinAlg::TensorGenerators::identity; + + const bool use_consistent_tangent = (tang_str != "elastic"); + + FOUR_C_ASSERT_ALWAYS(strainbarpllast_.at(gp) >= 0.0, + "accumulated plastic strain has to be equal to or greater than zero!"); + + // Elastic predictor + const Core::LinAlg::SymmetricTensor trialstrain_e = glstrain - strainpllast_.at(gp); + const double tracestrain = trace(trialstrain_e); + const Core::LinAlg::SymmetricTensor devstrain = + trialstrain_e - id2 * (tracestrain / 3.0); + + double p = kappa * tracestrain; + const double p_trial = p; + Core::LinAlg::SymmetricTensor devstress = 2.0 * G * devstrain; + + const double J2 = 0.5 * ddot(devstress, devstress); + double strainbar_p = strainbarpllast_.at(gp); + const double Phi_trial = std::sqrt(J2) + eta * p_trial - xi * cohesion - xi * Hiso * strainbar_p; + + double Dgamma = 0.0; + double dstrainv = 0.0; + + if (Phi_trial / std::abs(cohesion) > params_->abstol_) { - auto returnToConeFunctAndDeriv = [this, &G, &kappa, &Phi_trial](ScalarT Dgamma_init) + // Return to cone + auto returnToConeFunctAndDeriv = [this, &G, &kappa, &Phi_trial](double Dgamma_init) { return this->return_to_cone_funct_and_deriv(Dgamma_init, G, kappa, Phi_trial); }; - const double tol = params_->abstol_; - Dgamma = - Core::Utils::solve_local_newton(returnToConeFunctAndDeriv, Dgamma, tol * cohesion, itermax); - strainbar_p = (strainbarpllast_.at(gp)) + xi * Dgamma; - devstress.scale(1.0 - (G * Dgamma / std::sqrt(J2))); - p = p_trial - kappa * etabar * Dgamma; - if ((std::sqrt(J2) - G * Dgamma) / abs(cohesion) < params_->abstol_) + Dgamma = Core::Utils::solve_local_newton( + returnToConeFunctAndDeriv, Dgamma, params_->abstol_ * cohesion, itermax); + + if ((std::sqrt(J2) - G * Dgamma) / std::abs(cohesion) < params_->abstol_) { - strainbar_p = (strainbarpllast_.at(gp)); - auto returnToApexFunctAndDeriv = [this, &p_trial, &kappa, &strainbar_p](ScalarT dstrainv_init) + // Return to apex + strainbar_p = strainbarpllast_.at(gp); + auto returnToApexFunctAndDeriv = [this, &p_trial, &kappa, &strainbar_p](double dstrainv_init) { return this->return_to_apex_funct_and_deriv(dstrainv_init, p_trial, kappa, strainbar_p); }; - const double tol = params_->abstol_; dstrainv = Core::Utils::solve_local_newton( - returnToApexFunctAndDeriv, dstrainv, tol * cohesion, itermax); - strainbar_p = (strainbarpllast_.at(gp)) + xi / eta * dstrainv; + returnToApexFunctAndDeriv, dstrainv, params_->abstol_ * cohesion, itermax); + strainbar_p = strainbarpllast_.at(gp) + xi / eta * dstrainv; p = p_trial - kappa * dstrainv; - for (int i = 0; i < 6; i++) devstress(i) = 0.0; + devstress.fill(0.0); } - PlasticDruckerPrager::stress(p, devstress, stress); - strain_e.update(1 / G / 2, devstress, p / kappa / 3, id2Scalar); - for (int i = 3; i < 6; ++i) strain_e(i) *= 2.0; - for (int i = 3; i < 6; ++i) strain(i) *= 2.0; - strain_p.update(1.0, strain, -1.0, strain_e); - - strainplcurr_.at(gp) = Core::FADUtils::cast_to_double(strain_p); - strainbarplcurr_.at(gp) = Core::FADUtils::cast_to_double(strainbar_p); + else + { + strainbar_p = strainbarpllast_.at(gp) + xi * Dgamma; + devstress *= (1.0 - G * Dgamma / std::sqrt(J2)); + p = p_trial - kappa * etabar * Dgamma; + } + + stress = devstress + p * id2; + + // Recover elastic strain from updated stress: strain_e = s/(2G) + p/(3*kappa)*I + const Core::LinAlg::SymmetricTensor strain_e = + devstress * (1.0 / (2.0 * G)) + id2 * (p / (3.0 * kappa)); + strainplcurr_.at(gp) = glstrain - strain_e; + strainbarplcurr_.at(gp) = strainbar_p; } else { - PlasticDruckerPrager::stress(p, devstress, stress); - strain_e.update(trialstrain_e); - for (int i = 3; i < 6; ++i) strain_e(i) *= 2.0; - for (int i = 3; i < 6; ++i) strain(i) *= 2.0; - + stress = devstress + p * id2; strainplcurr_.at(gp) = strainpllast_.at(gp); strainbarplcurr_.at(gp) = strainbarpllast_.at(gp); } - if ((Phi_trial > 0) && (tang == 1)) + + if ((Phi_trial > 0.0) && use_consistent_tangent) { - Core::LinAlg::Matrix devstraindouble = - Core::FADUtils::cast_to_double(devstrain); if (dstrainv != 0.0) - { - setup_cmat_elasto_plastic_apex(cmat, Core::FADUtils::cast_to_double(kappa), devstraindouble, - Core::FADUtils::cast_to_double(xi), Core::FADUtils::cast_to_double(Hiso), - Core::FADUtils::cast_to_double(eta), Core::FADUtils::cast_to_double(etabar)); - } + setup_cmat_elasto_plastic_apex(cmat, kappa, xi, Hiso, eta, etabar); else - { - setup_cmat_elasto_plastic_cone(cmat, Core::FADUtils::cast_to_double(Dgamma), - Core::FADUtils::cast_to_double(G), Core::FADUtils::cast_to_double(kappa), devstraindouble, - Core::FADUtils::cast_to_double(xi), Core::FADUtils::cast_to_double(Hiso), - Core::FADUtils::cast_to_double(eta), Core::FADUtils::cast_to_double(etabar)); - } + setup_cmat_elasto_plastic_cone(cmat, Dgamma, G, kappa, devstrain, xi, Hiso, eta, etabar); } else { @@ -414,10 +315,11 @@ bool Mat::PlasticDruckerPrager::evaluate_output_data( { for (std::size_t gp = 0; gp < strainplcurr_.size(); ++gp) { - const double* values = strainplcurr_.at(gp).data(); + const Core::LinAlg::Matrix<6, 1> voigt_strain = + Core::LinAlg::make_strain_like_voigt_matrix(strainplcurr_.at(gp)); for (std::size_t i = 0; i < 6; ++i) { - data(gp, i) = values[i]; + data(gp, i) = voigt_strain(i, 0); } } return true; @@ -425,68 +327,33 @@ bool Mat::PlasticDruckerPrager::evaluate_output_data( return false; } -template -void Mat::PlasticDruckerPrager::stress(const T p, - const Core::LinAlg::Matrix& devstress, - Core::LinAlg::Matrix& stress) const -{ - stress.update(devstress); - for (int i = 0; i < 3; ++i) stress(i) += p; -} - -template -std::pair Mat::PlasticDruckerPrager::return_to_cone_funct_and_deriv( - T Dgamma, T G, T kappa, T Phi_trial) +std::pair Mat::PlasticDruckerPrager::return_to_cone_funct_and_deriv( + double Dgamma, double G, double kappa, double Phi_trial) { - double Hiso = params_->isohard_; - double eta = params_->eta_; - double xi = params_->xi_; - double etabar = params_->etabar_; - T Res = Phi_trial - Dgamma * (G + eta * kappa * etabar) - (xi * xi * Dgamma * Hiso); - T d = -G - (kappa * etabar * eta) - (xi * xi * Hiso); + const double Hiso = params_->isohard_; + const double eta = params_->eta_; + const double xi = params_->xi_; + const double etabar = params_->etabar_; + double Res = Phi_trial - Dgamma * (G + eta * kappa * etabar) - (xi * xi * Dgamma * Hiso); + double d = -G - (kappa * etabar * eta) - (xi * xi * Hiso); return {Res, d}; } -template -std::pair Mat::PlasticDruckerPrager::return_to_apex_funct_and_deriv( - T dstrainv, T p, T kappa, T strainbar_p) +std::pair Mat::PlasticDruckerPrager::return_to_apex_funct_and_deriv( + double dstrainv, double p, double kappa, double strainbar_p) { - double Hiso = params_->isohard_; - double eta = params_->eta_; - double xi = params_->xi_; - double cohesion = params_->cohesion_; - double etabar = params_->etabar_; - double alpha = xi / eta; - double beta = xi / etabar; - T Res = + const double Hiso = params_->isohard_; + const double eta = params_->eta_; + const double xi = params_->xi_; + const double cohesion = params_->cohesion_; + const double etabar = params_->etabar_; + const double alpha = xi / eta; + const double beta = xi / etabar; + double Res = beta * cohesion + beta * strainbar_p * Hiso - p + dstrainv * (alpha * beta * Hiso + kappa); - T d = xi * xi / eta / etabar * Hiso + kappa; + double d = xi * xi / eta / etabar * Hiso + kappa; return {Res, d}; } -template void Mat::PlasticDruckerPrager::evaluate_fad(const Core::LinAlg::Matrix<3, 3>*, - const Core::LinAlg::Matrix<6, 1, double>&, const Teuchos::ParameterList&, - Core::LinAlg::Matrix<6, 1, double>&, Core::LinAlg::SymmetricTensor&, int gp, - int eleGID); -template void Mat::PlasticDruckerPrager::evaluate_fad(const Core::LinAlg::Matrix<3, 3>*, - const Core::LinAlg::Matrix<6, 1, FAD>&, const Teuchos::ParameterList&, - Core::LinAlg::Matrix<6, 1, FAD>&, Core::LinAlg::SymmetricTensor&, int gp, - int eleGID); -template void Mat::PlasticDruckerPrager::stress(const double p, - const Core::LinAlg::Matrix& devstress, - Core::LinAlg::Matrix& stress) const; -template void Mat::PlasticDruckerPrager::stress(const FAD p, - const Core::LinAlg::Matrix& devstress, - Core::LinAlg::Matrix& stress) const; -template std::pair -Mat::PlasticDruckerPrager::return_to_cone_funct_and_deriv( - double Dgamma, double G, double kappa, double Phi_trial); -template std::pair Mat::PlasticDruckerPrager::return_to_cone_funct_and_deriv( - FAD Dgamma, FAD G, FAD kappa, FAD Phi_trial); -template std::pair -Mat::PlasticDruckerPrager::return_to_apex_funct_and_deriv( - double dstrainv, double p, double kappa, double strainbar_p); -template std::pair Mat::PlasticDruckerPrager::return_to_apex_funct_and_deriv( - FAD dstrainv, FAD p, FAD kappa, FAD strainbar_p); FOUR_C_NAMESPACE_CLOSE diff --git a/src/mat/4C_mat_plasticdruckerprager.hpp b/src/mat/4C_mat_plasticdruckerprager.hpp index 66c00986900..bebdce8d2c9 100644 --- a/src/mat/4C_mat_plasticdruckerprager.hpp +++ b/src/mat/4C_mat_plasticdruckerprager.hpp @@ -13,10 +13,9 @@ #include "4C_comm_parobjectfactory.hpp" #include "4C_inpar_structure.hpp" #include "4C_linalg_tensor_conversion.hpp" -#include "4C_mat_material_factory.hpp" +#include "4C_linalg_tensor_generators.hpp" #include "4C_mat_so3_material.hpp" #include "4C_material_parameter_base.hpp" -#include "4C_utils_local_newton.hpp" FOUR_C_NAMESPACE_OPEN @@ -24,7 +23,6 @@ namespace Mat { namespace PAR { - using FAD = Sacado::Fad::DFad; /** * \brief Elasto-plasitc Drucker-Prager Material Model * @@ -123,82 +121,45 @@ namespace Mat const Core::LinAlg::SymmetricTensor& glstrain, const Teuchos::ParameterList& params, const EvaluationContext<3>& context, Core::LinAlg::SymmetricTensor& stress, - Core::LinAlg::SymmetricTensor& cmat, int gp, int eleGID) override - { - Core::LinAlg::Matrix<6, 1> stress_view = Core::LinAlg::make_stress_like_voigt_view(stress); - Core::LinAlg::Matrix<3, 3> defgrad_view; - Core::LinAlg::Matrix<3, 3>* defgrd_ptr = nullptr; - if (defgrad) - { - defgrad_view = Core::LinAlg::make_matrix_view(*defgrad); - defgrd_ptr = &defgrad_view; - } - this->evaluate_fad(defgrd_ptr, Core::LinAlg::make_strain_like_voigt_matrix(glstrain), params, - stress_view, cmat, gp, eleGID); - }; - - template - void evaluate(const Core::LinAlg::Matrix<3, 3>* defgrd, - const Core::LinAlg::Matrix* linstrain, - const Teuchos::ParameterList& params, - Core::LinAlg::Matrix* stress, - Core::LinAlg::SymmetricTensor* cmat, int gp, int eleGID) - { - this->evaluate_fad(defgrd, *linstrain, params, *stress, *cmat, gp, eleGID); - }; - template - void evaluate_fad(const Core::LinAlg::Matrix<3, 3>* defgrd, - const Core::LinAlg::Matrix& linstrain, - const Teuchos::ParameterList& params, - Core::LinAlg::Matrix& stress, - Core::LinAlg::SymmetricTensor& cmat, int gp, int eleGID); - template - void stress(const T p, const Core::LinAlg::Matrix& devstress, - Core::LinAlg::Matrix& stress) const; + Core::LinAlg::SymmetricTensor& cmat, int gp, int eleGID) override; void setup_cmat(Core::LinAlg::SymmetricTensor& cmat) const; /** - * \brief setup the elastoplasticity tensor in matrix notation for 3d return to cone + * \brief setup the elastoplasticity tensor for 3d return to cone * * \param cmat :elasto-plastic tangent modulus (out) * \param Dgamma :plastic multiplier * \param G :shear modulus * \param Kappa :Bulk modulus - * \param devstrain :deviatoric strain + * \param devstrain :deviatoric strain tensor * \param xi :Mohr-Coulomb parameter * \param Hiso :isotropic hardening modulus * \param eta :Mohr-Coulomb parameter * \param etabar :Mohr-Coulomb parameter */ void setup_cmat_elasto_plastic_cone(Core::LinAlg::SymmetricTensor& cmat, - double Dgamma, double G, double Kappa, Core::LinAlg::Matrix& devstrain, - double xi, double Hiso, double eta, double etabar) const; + double Dgamma, double G, double Kappa, + const Core::LinAlg::SymmetricTensor& devstrain, double xi, double Hiso, + double eta, double etabar) const; /** - * \brief setup the elastoplasticity tensor in matrix notation for 3d return to apex + * \brief setup the elastoplasticity tensor for 3d return to apex * * \param cmat :elasto-plastic tangent modulus (out) * \param Kappa :Bulk modulus - * \param devstrain :deviatoric strain * \param xi :Mohr-Coulomb parameter * \param Hiso :isotropic hardening modulus * \param eta :Mohr-Coulomb parameter * \param etabar :Mohr-Coulomb parameter */ - void setup_cmat_elasto_plastic_apex(Core::LinAlg::SymmetricTensor& - cmat, // elasto-plastic tangent modulus (out) - double Kappa, // Bulk modulus - Core::LinAlg::Matrix& devstrain, // deviatoric strain - double xi, // Mohr-Coulomb parameter - double Hiso, // isotropic hardening modulus - double eta, // Mohr-Coulomb parameter - double etabar // Mohr-Coulomb parameter - ) const; + void setup_cmat_elasto_plastic_apex(Core::LinAlg::SymmetricTensor& cmat, + double Kappa, double xi, double Hiso, double eta, double etabar) const; + Core::Mat::PAR::Parameter* parameter() const override { return params_; } double density() const override { return params_->density_; } - template - std::pair return_to_cone_funct_and_deriv(T Dgamma, T G, T kappa, T Phi_trial); - template - std::pair return_to_apex_funct_and_deriv(T dstrainv, T p, T kappa, T strainbar_p); + std::pair return_to_cone_funct_and_deriv( + double Dgamma, double G, double kappa, double Phi_trial); + std::pair return_to_apex_funct_and_deriv( + double dstrainv, double p, double kappa, double strainbar_p); bool initialized() const { return (isinit_ and !strainplcurr_.empty()); } void register_output_data_names( std::unordered_map& names_and_size) const override; @@ -208,8 +169,8 @@ namespace Mat private: Mat::PAR::PlasticDruckerPrager* params_; - std::vector> strainpllast_; - std::vector> strainplcurr_; + std::vector> strainpllast_; + std::vector> strainplcurr_; std::vector strainbarpllast_; std::vector strainbarplcurr_; bool isinit_; diff --git a/unittests/mat/4C_druckerprager_test.cpp b/unittests/mat/4C_druckerprager_test.cpp index b60a0af40b2..142b055d6e0 100644 --- a/unittests/mat/4C_druckerprager_test.cpp +++ b/unittests/mat/4C_druckerprager_test.cpp @@ -9,7 +9,6 @@ #include "4C_comm_pack_buffer.hpp" #include "4C_global_data.hpp" -#include "4C_linalg_FADmatrix_utils.hpp" #include "4C_linalg_tensor_conversion.hpp" #include "4C_linalg_tensor_generators.hpp" #include "4C_mat_material_factory.hpp" @@ -194,51 +193,38 @@ namespace TEST_F(DruckerPragerTest, TestEvaluateHistory) { druckprag_->setup(1, {}, {}); - Core::LinAlg::Matrix<6, 1, FAD> input_strain; - for (int i = 0; i < 3; ++i) input_strain(i) = FAD(6, i, 0.1); - for (int i = 3; i < 6; ++i) input_strain(i) = FAD(6, i, 0.1); + // Step 1: plastic loading + Core::LinAlg::SymmetricTensor input_strain{}; + input_strain(0, 0) = 0.1; + input_strain(1, 1) = 0.1; + input_strain(2, 2) = 0.1; + input_strain(0, 1) = 0.05; + input_strain(1, 2) = 0.05; + input_strain(0, 2) = 0.05; Teuchos::ParameterList paras; - Core::LinAlg::Matrix<3, 3> defgrad(Core::LinAlg::Initialization::zero); - Core::LinAlg::Matrix<6, 1, FAD> ref_stress(Core::LinAlg::Initialization::zero); Core::LinAlg::SymmetricTensor result_cmat{}; - Core::LinAlg::Matrix<6, 6> result_cmat_view = - Core::LinAlg::make_stress_like_voigt_view(result_cmat); - Core::LinAlg::Matrix<6, 1, FAD> result_stress(Core::LinAlg::Initialization::zero); - druckprag_->evaluate(&defgrad, &input_strain, paras, &result_stress, &result_cmat, 0, 0); - Core::LinAlg::Matrix<6, 6> ref_cmat(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 6; i++) - { - for (int j = 0; j < 6; j++) - { - ref_cmat(i, j) = result_stress(i).dx(j); - } - } - FOUR_C_EXPECT_NEAR(result_cmat_view, ref_cmat, 1.0e-12); + Core::LinAlg::SymmetricTensor result_stress{}; + double total_time = 0.0; + double time_step_size = 1.0; + Mat::EvaluationContext<3> context{.total_time = &total_time, + .time_step_size = &time_step_size, + .xi = {}, + .ref_coords = nullptr}; + druckprag_->evaluate(nullptr, input_strain, paras, context, result_stress, result_cmat, 0, 0); + + // Step 2: further plastic loading druckprag_->update(); - for (int i = 0; i < 3; ++i) input_strain(i) = FAD(6, i, 1.0); - for (int i = 3; i < 6; ++i) input_strain(i) = FAD(6, i, 0.0); - druckprag_->evaluate(&defgrad, &input_strain, paras, &result_stress, &result_cmat, 0, 0); - for (int i = 0; i < 6; i++) - { - for (int j = 0; j < 6; j++) - { - ref_cmat(i, j) = result_stress(i).dx(j); - } - } - FOUR_C_EXPECT_NEAR(result_cmat_view, ref_cmat, 1.0e-12); + input_strain = Core::LinAlg::TensorGenerators::identity; + druckprag_->evaluate(nullptr, input_strain, paras, context, result_stress, result_cmat, 0, 0); + + // Step 3: elastic unloading -- stress should be zero for zero strain increment beyond plastic druckprag_->update(); - for (int i = 0; i < 3; ++i) input_strain(i) = FAD(6, i, 0.2); - for (int i = 3; i < 6; ++i) input_strain(i) = FAD(6, i, 0.0); - druckprag_->evaluate(&defgrad, &input_strain, paras, &result_stress, &result_cmat, 0, 0); - for (int i = 0; i < 6; i++) - { - for (int j = 0; j < 6; j++) - { - ref_cmat(i, j) = result_stress(i).dx(j); - } - } - FOUR_C_EXPECT_NEAR(Core::FADUtils::cast_to_double(result_stress), - Core::FADUtils::cast_to_double(ref_stress), 1.0e-12); + input_strain = 0.2 * Core::LinAlg::TensorGenerators::identity; + druckprag_->evaluate(nullptr, input_strain, paras, context, result_stress, result_cmat, 0, 0); + // After large plastic loading, a small strain should be elastic unloading + // Verify stress is non-zero (elastic response from residual strain) + double stress_norm = std::sqrt(ddot(result_stress, result_stress)); + EXPECT_GT(stress_norm, 0.0); }; //! test member function Evaluate for arbitrary values @@ -272,118 +258,135 @@ namespace FOUR_C_EXPECT_NEAR(result_stress, ref_stress, 1.0e-12); }; - //! test member function Evaluate - TEST_F(DruckerPragerTest, TestEvaluateCmat) + //! helper: compute tangent via finite differences + Core::LinAlg::SymmetricTensor fd_tangent(Mat::PlasticDruckerPrager& mat, + const Core::LinAlg::SymmetricTensor& strain_base, + const Teuchos::ParameterList& paras, const Mat::EvaluationContext<3>& context) { - druckprag_->setup(1, {}, {}); - Core::LinAlg::Matrix<6, 1, FAD> input_strain; - for (int i = 0; i < 6; ++i) input_strain(i) = FAD(6, i, .1 * i); - Teuchos::ParameterList paras; - Core::LinAlg::Matrix<3, 3> defgrad(Core::LinAlg::Initialization::zero); - Core::LinAlg::Matrix<6, 1, FAD> ref_stress(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 3; ++i) - ref_stress(i) = - FAD((1.0 / ((1.0 + 0.25) * (1.0 - (2.0 * 0.25)))) * ((1.0 - 0.25) + 0.25 + 0.25) * .1); - for (int i = 3; i < 6; ++i) - ref_stress(i) = - FAD((1.0 / ((1.0 + 0.25) * (1.0 - (2.0 * 0.25)))) * ((1.0 - (2.0 * 0.25)) / 2.0) * .1); - Core::LinAlg::SymmetricTensor result_cmat{}; - Core::LinAlg::Matrix<6, 6> result_cmat_view = - Core::LinAlg::make_stress_like_voigt_view(result_cmat); - Core::LinAlg::Matrix<6, 1, FAD> result_stress(Core::LinAlg::Initialization::zero); - druckprag_->evaluate(&defgrad, &input_strain, paras, &result_stress, &result_cmat, 0, 0); - Core::LinAlg::Matrix<6, 6> ref_cmat(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 6; i++) + const double delta = 1.0e-7; + Core::LinAlg::SymmetricTensor fd_cmat{}; + + // Base stress + Core::LinAlg::SymmetricTensor stress_base{}; + Core::LinAlg::SymmetricTensor cmat_dummy{}; + mat.evaluate(nullptr, strain_base, paras, context, stress_base, cmat_dummy, 0, 0); + + // Perturb each independent strain component + for (int k = 0; k < 3; ++k) { - for (int j = 0; j < 6; j++) + for (int l = k; l < 3; ++l) { - ref_cmat(i, j) = result_stress(i).dx(j); + auto strain_pert = strain_base; + strain_pert(k, l) += delta; + + // Need a fresh material state for each perturbation (same history as base) + Core::LinAlg::SymmetricTensor stress_pert{}; + Core::LinAlg::SymmetricTensor cmat_pert{}; + mat.evaluate(nullptr, strain_pert, paras, context, stress_pert, cmat_pert, 0, 0); + + const double effective_delta = (k == l) ? delta : 2.0 * delta; + for (int i = 0; i < 3; ++i) + for (int j = i; j < 3; ++j) + fd_cmat(i, j, k, l) = (stress_pert(i, j) - stress_base(i, j)) / effective_delta; } } - FOUR_C_EXPECT_NEAR(result_cmat_view, ref_cmat, 1.0e-12); + return fd_cmat; + } + + //! test member function Evaluate for elastic tangent via FD + TEST_F(DruckerPragerTest, TestEvaluateCmat) + { + druckprag_->setup(1, {}, {}); + Core::LinAlg::SymmetricTensor input_strain = + Core::LinAlg::TensorGenerators::full<3, 3>(0.05) + + 0.05 * Core::LinAlg::TensorGenerators::identity; + Teuchos::ParameterList paras; + double total_time = 0.0; + double time_step_size = 1.0; + Mat::EvaluationContext<3> context{.total_time = &total_time, + .time_step_size = &time_step_size, + .xi = {}, + .ref_coords = nullptr}; + + Core::LinAlg::SymmetricTensor result_cmat{}; + Core::LinAlg::SymmetricTensor result_stress{}; + druckprag_->evaluate(nullptr, input_strain, paras, context, result_stress, result_cmat, 0, 0); + + auto ref_cmat = fd_tangent(*druckprag_, input_strain, paras, context); + FOUR_C_EXPECT_NEAR(result_cmat, ref_cmat, 1.0e-5); }; //! test CMAT matrix for Return to Cone TEST_F(DruckerPragerTest, TestEvaluateReturnToConeCmat) { druckprag_->setup(1, {}, {}); - Core::LinAlg::Matrix<6, 1, FAD> input_strain; - for (int i = 0; i < 3; ++i) input_strain(i) = FAD(6, i, 0.1 * i); - for (int i = 3; i < 6; ++i) input_strain(i) = FAD(6, i, 2.2 * i); + Core::LinAlg::SymmetricTensor input_strain{}; + input_strain(0, 0) = 0.0; + input_strain(1, 1) = 0.1; + input_strain(2, 2) = 0.2; + input_strain(0, 1) = 1.1; + input_strain(1, 2) = 2.2; + input_strain(0, 2) = 3.3; Teuchos::ParameterList paras; - Core::LinAlg::Matrix<3, 3> defgrad(Core::LinAlg::Initialization::zero); + double total_time = 0.0; + double time_step_size = 1.0; + Mat::EvaluationContext<3> context{.total_time = &total_time, + .time_step_size = &time_step_size, + .xi = {}, + .ref_coords = nullptr}; + Core::LinAlg::SymmetricTensor result_cmat{}; - Core::LinAlg::Matrix<6, 6> result_cmat_view = - Core::LinAlg::make_stress_like_voigt_view(result_cmat); - Core::LinAlg::Matrix<6, 1, FAD> result_stress(Core::LinAlg::Initialization::zero); - druckprag_->evaluate(&defgrad, &input_strain, paras, &result_stress, &result_cmat, 0, 0); - Core::LinAlg::Matrix<6, 6> ref_cmat(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 6; i++) - { - for (int j = 0; j < 6; j++) - { - ref_cmat(i, j) = result_stress(i).dx(j); - } - } - FOUR_C_EXPECT_NEAR(result_cmat_view, ref_cmat, 1.0e-12); + Core::LinAlg::SymmetricTensor result_stress{}; + druckprag_->evaluate(nullptr, input_strain, paras, context, result_stress, result_cmat, 0, 0); + + auto ref_cmat = fd_tangent(*druckprag_, input_strain, paras, context); + FOUR_C_EXPECT_NEAR(result_cmat, ref_cmat, 1.0e-5); }; TEST_F(DruckerPragerTest, TestEvaluateReturnToApexCmat) { druckprag_->setup(1, {}, {}); - Core::LinAlg::Matrix<6, 1, FAD> input_strain; - for (int i = 0; i < 3; ++i) input_strain(i) = FAD(6, i, 1.0); - for (int i = 3; i < 6; ++i) input_strain(i) = FAD(6, i, 0.0); + Core::LinAlg::SymmetricTensor input_strain = + Core::LinAlg::TensorGenerators::identity; Teuchos::ParameterList paras; - Core::LinAlg::Matrix<3, 3> defgrad(Core::LinAlg::Initialization::zero); + double total_time = 0.0; + double time_step_size = 1.0; + Mat::EvaluationContext<3> context{.total_time = &total_time, + .time_step_size = &time_step_size, + .xi = {}, + .ref_coords = nullptr}; + Core::LinAlg::SymmetricTensor result_cmat{}; - Core::LinAlg::Matrix<6, 6> result_cmat_view = - Core::LinAlg::make_stress_like_voigt_view(result_cmat); - Core::LinAlg::Matrix<6, 1, FAD> result_stress(Core::LinAlg::Initialization::zero); - druckprag_->evaluate(&defgrad, &input_strain, paras, &result_stress, &result_cmat, 0, 0); - Core::LinAlg::Matrix<6, 6> ref_cmat(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 6; i++) - { - for (int j = 0; j < 6; j++) - { - ref_cmat(i, j) = result_stress(i).dx(j); - } - } - FOUR_C_EXPECT_NEAR(result_cmat_view, ref_cmat, 1.0e-12); + Core::LinAlg::SymmetricTensor result_stress{}; + druckprag_->evaluate(nullptr, input_strain, paras, context, result_stress, result_cmat, 0, 0); + + auto ref_cmat = fd_tangent(*druckprag_, input_strain, paras, context); + FOUR_C_EXPECT_NEAR(result_cmat, ref_cmat, 1.0e-5); }; //! test CMAT matrix for Return to Apex TEST_F(DruckerPragerTest, TestEvaluateRandomStrainCmat) { druckprag_->setup(1, {}, {}); - Core::LinAlg::Matrix<6, 1, FAD> input_strain; - input_strain(0) = FAD(6, 0, 1.1); - input_strain(1) = FAD(6, 1, 2.0); - input_strain(2) = FAD(6, 2, 0.1); - input_strain(3) = FAD(6, 3, 2.5); - input_strain(4) = FAD(6, 4, 1.4); - input_strain(5) = FAD(6, 5, 1.0); + Core::LinAlg::SymmetricTensor input_strain; + input_strain(0, 0) = 1.1; + input_strain(1, 1) = 2.0; + input_strain(2, 2) = 0.1; + input_strain(0, 1) = 2.5 / 2; + input_strain(1, 2) = 1.4 / 2; + input_strain(0, 2) = 1.0 / 2; Teuchos::ParameterList paras; - Core::LinAlg::Matrix<3, 3> defgrad(Core::LinAlg::Initialization::zero); - Core::LinAlg::Matrix<6, 1, FAD> ref_stress(Core::LinAlg::Initialization::zero); - ref_stress(0) = FAD(1.4142412329012); - ref_stress(1) = FAD(1.8571566160540); - ref_stress(2) = FAD(0.9221130293981); - ref_stress(3) = FAD(0.6151602543789); - ref_stress(4) = FAD(0.3444897424522); - ref_stress(5) = FAD(0.2460641017516); + double total_time = 0.0; + double time_step_size = 1.0; + Mat::EvaluationContext<3> context{.total_time = &total_time, + .time_step_size = &time_step_size, + .xi = {}, + .ref_coords = nullptr}; + Core::LinAlg::SymmetricTensor result_cmat{}; - Core::LinAlg::Matrix<6, 6> result_cmat_view = - Core::LinAlg::make_stress_like_voigt_view(result_cmat); - Core::LinAlg::Matrix<6, 1, FAD> result_stress(Core::LinAlg::Initialization::zero); - druckprag_->evaluate(&defgrad, &input_strain, paras, &result_stress, &result_cmat, 0, 0); - Core::LinAlg::Matrix<6, 6> ref_cmat(Core::LinAlg::Initialization::zero); - for (int i = 0; i < 6; i++) - { - for (int j = 0; j < 6; j++) - { - ref_cmat(i, j) = result_stress(i).dx(j); - } - } - FOUR_C_EXPECT_NEAR(result_cmat_view, ref_cmat, 1.0e-12); + Core::LinAlg::SymmetricTensor result_stress{}; + druckprag_->evaluate(nullptr, input_strain, paras, context, result_stress, result_cmat, 0, 0); + + auto ref_cmat = fd_tangent(*druckprag_, input_strain, paras, context); + FOUR_C_EXPECT_NEAR(result_cmat, ref_cmat, 1.0e-5); }; } // namespace