[GeoMechanicsApplication] Convert the reordering of Mohr-Coulomb principal stresses to averaging by mnabideltares · Pull Request #13595 · KratosMultiphysics/Kratos

mnabideltares · 2025-07-04T11:16:08Z

📝 Description
In order to make the currently implemented Mohr-Coulomb consistent to UDSM of Plaxis, we need to change the reordering of the principal stresse to averaging type.

🆕 Changelog

Averaging is implemented replacing ot reordering
The mapping directions are averaged

With returning to the original trial stresses, average the stresses, average the mapping direction, then map it again

markelov208

Hi Mohamed, finally your hard work can be used by others. Here are my comments. Hope some of them may be useful.

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.cpp

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

...icsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_yield_surface.cpp

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

...tions/GeoMechanicsApplication/custom_constitutive/interface_coulomb_with_tension_cut_off.cpp

…hen needed.

… stresses only once.

…sses now outside if else condition again

…ulomb yield surface and the classes using coulomb yield surface.

…ment.

markelov208

Hi Wijtze Pieter, thank you very much for the changes, especially, enum class. They made the code is cleaner and more understandable. I have a few minor comments.

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

markelov208 · 2025-07-09T09:35:13Z

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

+                mCoulombWithTensionCutOffImpl.DoReturnMapping(r_prop, trial_sigma_tau, averaging_type);
+            mapped_principal_stress_vector = StressStrainUtilities::TransformSigmaTauToPrincipalStresses(
+                mapped_sigma_tau, averaged_principal_trial_stress_vector);
+            mapped_principal_stress_vector[1] = mapped_principal_stress_vector[static_cast<std::size_t>(averaging_type)];


I am afraid using the static_cast is error prone. I would suggest to use switch with enum. Then there is no need to assign enum values explicitly,

changed to a switch with enum

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

...icsApplication/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp

markelov208 · 2025-07-09T10:04:48Z

applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_yield_surface.cpp

+    auto           p_coulomb_yield_surface =
+        std::make_unique<CoulombYieldSurface>(friction_angle, cohesion, dilatancy_angle);


I am just curious, what is the reason for this change? It requires change in line 66.

This was done for a previous different attempt of refactoring, I have reverted the change now.

…ction instead.

rfaasse

Thank you for all the work on this (both the coding as well as all investigations and debugging sessions). It looks clear to me, although I don't know how to check all formulations etc.

I have some minor suggestions, but I think the main question is if we need to change the documentation

applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.cpp

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

rfaasse · 2025-07-09T14:35:23Z

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

+                averaged_principal_trial_stress_vector);
+            mapped_sigma_tau =
+                mCoulombWithTensionCutOffImpl.DoReturnMapping(r_prop, trial_sigma_tau, averaging_type);
+            mapped_principal_stress_vector = StressStrainUtilities::TransformSigmaTauToPrincipalStresses(


When there is averaging, the TransformSigmaTauToPrincipalStresses function is called twice and the result of the first call is overwritten. We might be able to remove the redundant call in case of averaging.

Before the result of the first call is overwritten, the result of the first call is used to determine ( with FindAveragingType() ) whether the averaging route should be taken. There is no redundant call.

I am not sure what you mean. We can discuss it when you are available.

applications/GeoMechanicsApplication/tests/test_mohr_coulomb_with_tension_cutoff.py

…hard.

WPK4FEM

I think we are ready.

WPK4FEM · 2025-07-07T10:56:18Z

applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp

+bool IsStressAtCornerReturnZone(const Vector& rTrialSigmaTau,
+                                const Vector& rDerivativeOfFlowFunction,
+                                const Vector& rCornerPoint)
 {
-    return rTrialSigmaTau[0] - rCornerPoint[0] -
-               (rTrialSigmaTau[1] - rCornerPoint[1]) * std::sin(DilatancyAngleInRadians) >=
+    return (rTrialSigmaTau[0] - rCornerPoint[0]) * rDerivativeOfFlowFunction[1] -
+               (rTrialSigmaTau[1] - rCornerPoint[1]) * rDerivativeOfFlowFunction[0] >=
           0.0;


This is now general, such that it can be used for all return directions.

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.cpp

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

WPK4FEM · 2025-07-10T08:40:16Z

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

+                averaged_principal_trial_stress_vector);
+            mapped_sigma_tau =
+                mCoulombWithTensionCutOffImpl.DoReturnMapping(r_prop, trial_sigma_tau, averaging_type);
+            mapped_principal_stress_vector = StressStrainUtilities::TransformSigmaTauToPrincipalStresses(


Before the result of the first call is overwritten, the result of the first call is used to determine ( with FindAveragingType() ) whether the averaging route should be taken. There is no redundant call.

applications/GeoMechanicsApplication/tests/test_mohr_coulomb_with_tension_cutoff.py

WPK4FEM

With the updated documentation, for me this is ready to go.

mnabideltares added 7 commits May 16, 2025 00:07

Applied averaging

e20e5d3

updated the reference for test case

123353e

Merge branch 'master' into geo/13371-convert-reordering-to-averaging-new

60635e2

Merge branch 'master' into geo/13371-convert-reordering-to-averaging-new

f949e95

Changed the mapping to be same as UDSM

2e4eecd

With returning to the original trial stresses, average the stresses, average the mapping direction, then map it again

refactoring 1

7d822b3

refactoring 2

0cf33f3

mnabideltares requested a review from a team as a code owner July 4, 2025 11:16

mnabideltares linked an issue Jul 4, 2025 that may be closed by this pull request

[GeoMechanicsApplication] Convert the reordering of Mohr-Coulomb principal stresses to averaging #13371

Closed

mnabideltares self-assigned this Jul 4, 2025

mnabideltares added the Feature label Jul 4, 2025

mnabideltares added this to GeoMechanicsApplication (Deltares) Jul 4, 2025

github-project-automation bot moved this to In progress in GeoMechanicsApplication (Deltares) Jul 4, 2025

mnabideltares requested review from WPK4FEM and avdg81 July 4, 2025 11:19

markelov208 reviewed Jul 4, 2025

View reviewed changes

Fixes based on review comments

79cc3ef

mnabideltares requested a review from markelov208 July 7, 2025 11:24

WPK4FEM added 12 commits July 7, 2025 14:38

refactoring: less operations when no yield function violated

5ee1a94

refactoring: negation removed from branching, actual averating only w…

6df01ee

…hen needed.

refactoring: no double to int conversion and compare mapped principal…

2347626

… stresses only once.

refactoring: renaming variables to reflect purpose

973eb6c

refactoring: clang tidy remarks

15c201a

refactoring: renaming variables to their purpose, RotatePrincipalStre…

3b0321f

…sses now outside if else condition again

refactoring: natural order and Sonarcloud codesmell unreachable break

6ccab5b

refactoring: variable naming and scope

03fa072

refactoring: misspelled trial

d9b40d3

refactoring: use std::fill for averaged principal stress setting

9a8b9ba

refactoring: mapping type --> averaging type, only implemented for co…

08a58e8

…ulomb yield surface and the classes using coulomb yield surface.

refactoring: the nature of AveragingType is enum behaviour.

c888d43

refactoring: first return mapping with constant, removed outdated com…

92bdda8

…ment.

markelov208 reviewed Jul 9, 2025

View reviewed changes

WPK4FEM added 4 commits July 9, 2025 14:16

most review comments by Gennady

fdb6071

removed empty statement and formatted.

595fb9c

refactoring: do not use the enum value as index, but a switch constru…

1cd462f

…ction instead.

refactoring: static functions to anonymous namespace.

7144edb

rfaasse reviewed Jul 9, 2025

View reviewed changes

WPK4FEM and others added 6 commits July 10, 2025 10:43

refactoring: replaced if branching with switch. Review comment of Ric…

813e5a9

…hard.

Documentation is modified based on the applied averaging method

edd24b1

pdf file of documentation is modified

96f1c9a

fixes in readme

c464ba5

another fix in readme

62f225d

fix

3a6c667

mnabideltares requested review from markelov208 and rfaasse July 10, 2025 12:41

WPK4FEM reviewed Jul 10, 2025

View reviewed changes

WPK4FEM approved these changes Jul 10, 2025

View reviewed changes

github-project-automation bot moved this from In progress to Reviewer approved in GeoMechanicsApplication (Deltares) Jul 10, 2025

markelov208 approved these changes Jul 10, 2025

View reviewed changes

WPK4FEM merged commit e99fb96 into master Jul 10, 2025
10 checks passed

WPK4FEM deleted the geo/13371-convert-reordering-to-averaging-new branch July 10, 2025 14:58

github-project-automation bot moved this from Reviewer approved to Done in GeoMechanicsApplication (Deltares) Jul 10, 2025

		auto p_coulomb_yield_surface =
		std::make_unique<CoulombYieldSurface>(friction_angle, cohesion, dilatancy_angle);

Conversation

mnabideltares commented Jul 4, 2025

Uh oh!

markelov208 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

markelov208 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rfaasse left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

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!

Uh oh!

WPK4FEM left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

WPK4FEM left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants