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.

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.

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).

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).

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.

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.

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.

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 from 2b993ef to 10b82d1 Compare April 29, 2026 10:12

Refactor the drucker prager model to use tensors

47c364d

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

amgebauer approved these changes May 5, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Refactor the drucker prager model to use tensors instead of matrices#2007

Refactor the drucker prager model to use tensors instead of matrices#2007
ischeider wants to merge 1 commit into4C-multiphysics:mainfrom
ischeider:refactor_druckerprager

ischeider commented Apr 29, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI Apr 29, 2026

Uh oh!

Copilot AI Apr 29, 2026

Uh oh!

Copilot AI Apr 29, 2026

Uh oh!

Copilot AI Apr 29, 2026

Uh oh!

amgebauer left a comment

Uh oh!

amgebauer May 5, 2026

Uh oh!

amgebauer May 5, 2026

Uh oh!

amgebauer May 5, 2026

Uh oh!

BishrMaradni May 5, 2026

Uh oh!

amgebauer May 5, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

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

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!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

amgebauer left a comment

Choose a reason for hiding this comment

Uh oh!

amgebauer May 5, 2026

Choose a reason for hiding this comment

Uh oh!

amgebauer May 5, 2026

Choose a reason for hiding this comment

Uh oh!

amgebauer May 5, 2026

Choose a reason for hiding this comment

Uh oh!

BishrMaradni May 5, 2026

Choose a reason for hiding this comment

Uh oh!

amgebauer May 5, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants