Refactor the drucker prager model to use tensors instead of matrices by ischeider · Pull Request #2007 · 4C-multiphysics/4C

ischeider · 2026-04-29T09:41:10Z

Description and Context

Now also the Drucker Prager Material model is supposed to use tensors instead of matrices.
The changes are significant though (all the peculiar scalar types are gone, and the FAD routines vanished then as well).
GitHub Copilot assisted a lot ;-)

Interested parties

@BishrMaradni

Copilot

Pull request overview

Refactors the structural Plastic Drucker–Prager material implementation to operate directly on Core::LinAlg tensor types (symmetric 2nd/4th order tensors) rather than Voigt matrices, aligning it with the codebase’s ongoing tensor-based constitutive model direction.

Changes:

Reworked PlasticDruckerPrager::evaluate() to use tensor operations end-to-end (stress, deviatoric split, J2, updates).
Updated history storage/serialization from Voigt vectors to SymmetricTensor<double,3,3>.
Reimplemented consistent tangent assembly for cone/apex returns using tensor generators (identity, symmetric_identity, dyadic, ddot).

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.

File	Description
`src/mat/4C_mat_plasticdruckerprager.hpp`	Updates public API signatures and internal history storage to tensor types; removes FAD-based interfaces.
`src/mat/4C_mat_plasticdruckerprager.cpp`	Implements tensor-based return mapping and tangent setup; updates pack/unpack and output extraction for tensor history.

Copilot · 2026-04-29T09:46:12Z

+  if (Phi_trial / std::abs(cohesion) > params_->abstol_)
  {


Phi_trial / std::abs(cohesion) will divide by zero when the input parameter C (cohesion) is 0.0, which is currently allowed by the material input spec. Consider using a cohesion-independent yield tolerance (e.g., compare Phi_trial against abstol_ * max(1, |cohesion|)), or explicitly validate/throw when cohesion == 0.0 if that case is unsupported.

input validators included.

Copilot · 2026-04-29T09:46:12Z

+    Dgamma = Core::Utils::solve_local_newton(
+        returnToConeFunctAndDeriv, Dgamma, params_->abstol_ * cohesion, itermax);
+    strainbar_p = strainbarpllast_.at(gp) + xi * Dgamma;


solve_local_newton(..., params_->abstol_ * cohesion, ...) can pass a zero or negative tolerance (if cohesion is 0 or negative), which breaks the Newton convergence criterion and can lead to non-termination until max-iterations/throw. Use a strictly positive tolerance (e.g., abstol_ * max(1, |cohesion|) or just abstol_) and/or validate cohesion > 0 up-front.

input validator included as well.

Copilot · 2026-04-29T09:46:12Z

+    strainbar_p = strainbarpllast_.at(gp) + xi * Dgamma;
+    devstress *= (1.0 - G * Dgamma / std::sqrt(J2));
    p = p_trial - kappa * etabar * Dgamma;


In the cone-return branch, devstress is rescaled by 1 - G*Dgamma/std::sqrt(J2). If the trial state is (nearly) hydrostatic then J2 == 0 and this will divide by zero before the code checks whether it should switch to the apex return. Consider detecting J2 near zero early and going directly to the apex return (or otherwise avoiding any 1/sqrt(J2) operations until J2 is known to be safely > 0).

Since sqrt(J2) - G*Dgamma/abs(cohesion) > abstol, J2 cannot be close to zero.

Copilot · 2026-04-29T09:46:12Z

+  const double nd = std::sqrt(ddot(devstrain, devstrain));
+  const auto D = devstrain / nd;



nd = std::sqrt(ddot(devstrain, devstrain)) can be zero for a purely volumetric (hydrostatic) strain state. The subsequent normalization D = devstrain / nd (and the 1/nd terms in epfac/epfac2) will then produce NaNs in the consistent tangent. Add a guard for nd near zero (e.g., treat as apex case, fall back to elastic tangent, or assert if cone return is not valid for nd==0).

This only happens in the cone region of plastic state, where the stress state cannot be hydrostatic.

amgebauer

Thank you @ischeider for migrating this material to tensors!

amgebauer · 2026-05-05T07:04:40Z

    for (int var = 0; var < histsize; ++var)
    {
-      Core::LinAlg::Matrix<NUM_STRESS_3D, 1> tmp_vector(Core::LinAlg::Initialization::zero);
+      Core::LinAlg::SymmetricTensor<double, 3, 3> 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);
    }


You could also directly extract (and pack) the full vectors instead of each component separately. This would simplify this pack logic a lot.

This is not clear to me. If you think I should pack/unpack the variables for all gauss points at once: I looked at other materials, but they all do it this way. Can you give any example for a leaner method?

You can just write

extract_from_pack(buffer, strainpllast_); extract_from_pack(buffer, strainbarpllast_); extract_from_pack(buffer, strainplcurr_); extract_from_pack(buffer, strainbarplcurr_);

and correspondingly for packing

add_to_pack(data, strainpllast_); add_to_pack(data, strainbarpllast_); add_to_pack(data, strainplcurr_); add_to_pack(data, strainbarplcurr_);

This avoids the manual loop over all gauss points.

Note, previously, this was often not possible. That's why other materials often also do this manual loop.

I agree, this would make packing/unpacking histsize unnecessary as well.

amgebauer · 2026-05-05T07:05:45Z

+  strainbarpllast_.assign(numgp, 0.0);
+  strainbarplcurr_.assign(numgp, 0.0);
+
+  for (int i = 0; i < numgp; ++i)
+  {
+    strainpllast_.at(i).fill(0.0);
+    strainplcurr_.at(i).fill(0.0);
+  }


I don't thunk that you need this, the vector should ale´ready be initialized with zero tensors.

amgebauer · 2026-05-05T07:08:31Z

  strainbarpllast_ = strainbarplcurr_;

-  std::for_each(strainplcurr_.begin(), strainplcurr_.end(), [](auto& item) { item.clear(); });
+  std::for_each(strainplcurr_.begin(), strainplcurr_.end(), [](auto& item) { item.fill(0.0); });


I don't know the details of the material, but do you really need to reset the plastic strain here? Often it makes sense to keep it as an initial value so that the local Newton algorithm is started close to the solution. But I haven't looked into the details of the material, so I might be wrong here.

Based on the reference they set it to 0 at setup, but I think this is not necessary as the plastic strain is calculated as the total strain - elastic strain. At the same time, the newton iterator is solving for the elastic component of the strain as far as I recall.

I don't think it's necessary either ;-)

amgebauer · 2026-05-05T07:11:26Z

+    input_strain(0, 0) = 1.0;
+    input_strain(1, 1) = 1.0;
+    input_strain(2, 2) = 1.0;
+    input_strain(0, 1) = 0.0;
+    input_strain(1, 2) = 0.0;
+    input_strain(0, 2) = 0.0;


Suggested change

input_strain(0, 0) = 1.0;

input_strain(1, 1) = 1.0;

input_strain(2, 2) = 1.0;

input_strain(0, 1) = 0.0;

input_strain(1, 2) = 0.0;

input_strain(0, 2) = 0.0;

input_strain = Core::LinAlg::TensorGenerators::identity<double, 3, 3>;

the same applies at some other places.

rjoussen · 2026-06-18T10:25:49Z

+                                          .default_value = 50,
+                                          .validator = positive<int>()}),
        },
        {.description = "elastic St.Venant Kirchhoff / plastic drucker prager"});


Suggested change

{.description = "elastic St.Venant Kirchhoff / plastic drucker prager"});

{.description = "Linear-elasto-plastic material with a Drucker-Prager yield surface: $\\phi = \\sqrt{J_2} + \\eta*\\mathrm{tr}(\\mathbf{\\sigma}) - \\xi\\cdot(c + H_\\mathrm{iso}\\varepsilon_\\mathrm{p}^\\mathrm{acc}$ and a potentially non-associated plastic potential: = \\sqrt{J_2} + \\overline{\\eta}*\\mathrm{tr}(\\mathbf{\\sigma}) - \\xi\\cdot(c + H_\\mathrm{iso}\\varepsilon_\\mathrm{p}^\\mathrm{acc}$"});

I believe this is what is happening in the material.

rjoussen · 2026-06-18T10:32:50Z

+                           .validator = in_range(excl(-1.0), excl(0.5))}),
+            parameter<double>(
+                "DENS", {.description = "Density", .validator = positive_or_zero<double>()}),
+            parameter<double>("ISOHARD", {.description = "linear isotropic hardening",


Suggested change

parameter<double>("ISOHARD", {.description = "linear isotropic hardening",

parameter<double>("ISOHARD", {.description = "Linear isotropic hardening modulus $H_\\mathrm{iso}$",

same for the other parameters.

rjoussen · 2026-06-18T11:59:05Z

            parameter<double>("ETA", {.description = "Drucker Prager Constant Eta"}),
            parameter<double>("XI", {.description = "Drucker Prager Constant Xi"}),
            parameter<double>("ETABAR", {.description = "Drucker Prager Constant Etabar"}),
            parameter<std::string>("TANG", {.description = "Method to compute the material tangent",


the "TANG" parameter should be an enum.

rjoussen · 2026-06-18T12:05:35Z

    for (int var = 0; var < histsize; ++var)
    {
-      Core::LinAlg::Matrix<NUM_STRESS_3D, 1> tmp_vector(Core::LinAlg::Initialization::zero);
+      Core::LinAlg::SymmetricTensor<double, 3, 3> 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);
    }


I agree, this would make packing/unpacking histsize unnecessary as well.

rjoussen · 2026-06-18T12:06:48Z

+  strainbarpllast_.assign(numgp, 0.0);
+  strainbarplcurr_.assign(numgp, 0.0);

  isinit_ = true;


Setup is called unconditionally, so isinit_ is not needed here I think. The variable can vanish.

rjoussen · 2026-06-18T12:07:57Z

+        double Dgamma, double G, double kappa, double Phi_trial);
+    std::pair<double, double> return_to_apex_funct_and_deriv(
+        double dstrainv, double p, double kappa, double strainbar_p);
    bool initialized() const { return (isinit_ and !strainplcurr_.empty()); }


I think the initialized() function is unnecessary

rjoussen · 2026-06-18T12:15:44Z

+    std::pair<double, double> return_to_cone_funct_and_deriv(
+        double Dgamma, double G, double kappa, double Phi_trial);
+    std::pair<double, double> return_to_apex_funct_and_deriv(
+        double dstrainv, double p, double kappa, double strainbar_p);


Suggested change

std::pair<double, double> return_to_cone_funct_and_deriv(

double Dgamma, double G, double kappa, double Phi_trial);

std::pair<double, double> return_to_apex_funct_and_deriv(

double dstrainv, double p, double kappa, double strainbar_p);

std::pair<double, double> return_to_cone_funct_and_deriv(

const double Dgamma, const double G, const double kappa, const double Phi_trial);

std::pair<double, double> return_to_apex_funct_and_deriv(

const double dstrainv, const double p, const double kappa, const double strainbar_p);

rjoussen · 2026-06-18T12:20:49Z

-        Core::LinAlg::Matrix<NUM_STRESS_3D, 1, T>& stress) const;
+        Core::LinAlg::SymmetricTensor<double, 3, 3, 3, 3>& cmat, int gp, int eleGID) override;

    void setup_cmat(Core::LinAlg::SymmetricTensor<double, 3, 3, 3, 3>& cmat) const;


I would remove this and use Mat::StVenantKirchhoff::evaluate_stress_linearization(young, nu) directly.

rjoussen · 2026-06-18T12:23:21Z

-    Core::LinAlg::Matrix<NUM_STRESS_3D, 1>& devstrain, double xi, double Hiso, double eta,
-    double etabar) const
+    const Core::LinAlg::SymmetricTensor<double, 3, 3>& devstrain, double xi, double Hiso,
+    double eta, double etabar) const


Suggested change

double eta, double etabar) const

Core::LinAlg::SymmetricTensor<double, 3, 3, 3, 3> Mat::PlasticDruckerPrager::setup_cmat_elasto_plastic_cone(

const double Dgamma, const double G, const double Kappa,

const Core::LinAlg::SymmetricTensor<double, 3, 3>& devstrain, const double xi, const double Hiso,

const double eta, const double etabar) const

I think this would be clearer, same for the other cmat function.

rjoussen · 2026-06-18T12:40:39Z

 #include "4C_global_data.hpp"
 #include "4C_inpar_structure.hpp"
-#include "4C_linalg_FADmatrix_utils.hpp"
 #include "4C_linalg_fixedsizematrix.hpp"


Suggested change

#include "4C_global_data.hpp"

#include "4C_inpar_structure.hpp"

#include "4C_linalg_FADmatrix_utils.hpp"

#include "4C_linalg_fixedsizematrix.hpp"

#include "4C_global_data.hpp"

#include "4C_linalg_fixedsizematrix.hpp"

#include "4C_inpar_structure.hpp" is not needed, same for #include "4C_utils_enum.hpp"

rjoussen · 2026-06-18T12:43:25Z

 #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"


Suggested change

#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_inpar_structure.hpp"

#include "4C_mat_so3_material.hpp"

#include "4C_linalg_tensor_conversion.hpp" and #include "4C_linalg_tensor_generators.hpp" are not needed.

ischeider requested a review from Copilot April 29, 2026 09:41

Copilot started reviewing on behalf of ischeider April 29, 2026 09:42 View session

Copilot AI reviewed Apr 29, 2026

View reviewed changes

ischeider force-pushed the refactor_druckerprager branch 2 times, most recently from 10b82d1 to 47c364d Compare April 29, 2026 10:54

amgebauer previously approved these changes May 5, 2026

View reviewed changes

rjoussen reviewed May 9, 2026

View reviewed changes

Comment thread src/mat/4C_mat_plasticdruckerprager.cpp Outdated

Refactor the drucker prager model to use tensors

77a9d0b

ischeider dismissed amgebauer’s stale review via 77a9d0b May 26, 2026 08:07

ischeider force-pushed the refactor_druckerprager branch from 47c364d to 77a9d0b Compare May 26, 2026 08:07

ischeider requested review from BishrMaradni and rjoussen May 26, 2026 14:34

ischeider self-assigned this May 26, 2026

rjoussen reviewed Jun 18, 2026

View reviewed changes

		const double nd = std::sqrt(ddot(devstrain, devstrain));
		const auto D = devstrain / nd;

	parameter<double>("ISOHARD", {.description = "linear isotropic hardening",
	parameter<double>("ISOHARD", {.description = "Linear isotropic hardening modulus $H_\\mathrm{iso}$",

-    double eta, double etabar) const
+Core::LinAlg::SymmetricTensor<double, 3, 3, 3, 3> Mat::PlasticDruckerPrager::setup_cmat_elasto_plastic_cone(
+    const double Dgamma, const double G, const double Kappa,
+    const Core::LinAlg::SymmetricTensor<double, 3, 3>& devstrain, const double xi, const double Hiso,
+    const double eta, const double etabar) const

Conversation

ischeider commented Apr 29, 2026

Description and Context

Interested parties

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

amgebauer left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

amgebauer Jun 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rjoussen Jun 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rjoussen Jun 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rjoussen Jun 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

amgebauer Jun 9, 2026 •

edited

Loading

rjoussen Jun 18, 2026 •

edited

Loading

rjoussen Jun 18, 2026 •

edited

Loading

rjoussen Jun 18, 2026 •

edited

Loading