Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
WPK4FEM merged 31 commits into from
Jul 10, 2025
Merged

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

Show file tree
Hide file tree
Changes from 25 commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
e20e5d3
Applied averaging
mnabideltares May 15, 2025
123353e
updated the reference for test case
mnabideltares May 16, 2025
60635e2
Merge branch 'master' into geo/13371-convert-reordering-to-averaging-new
mnabideltares May 19, 2025
f949e95
Merge branch 'master' into geo/13371-convert-reordering-to-averaging-new
mnabideltares Jun 27, 2025
2e4eecd
Changed the mapping to be same as UDSM
mnabideltares Jul 4, 2025
7d822b3
refactoring 1
mnabideltares Jul 4, 2025
0cf33f3
refactoring 2
mnabideltares Jul 4, 2025
79cc3ef
Fixes based on review comments
mnabideltares Jul 7, 2025
5ee1a94
refactoring: less operations when no yield function violated
WPK4FEM Jul 7, 2025
6df01ee
refactoring: negation removed from branching, actual averating only w…
WPK4FEM Jul 7, 2025
2347626
refactoring: no double to int conversion and compare mapped principal…
WPK4FEM Jul 7, 2025
973eb6c
refactoring: renaming variables to reflect purpose
WPK4FEM Jul 7, 2025
15c201a
refactoring: clang tidy remarks
WPK4FEM Jul 7, 2025
3b0321f
refactoring: renaming variables to their purpose, RotatePrincipalStre…
WPK4FEM Jul 7, 2025
6ccab5b
refactoring: natural order and Sonarcloud codesmell unreachable break
WPK4FEM Jul 8, 2025
03fa072
refactoring: variable naming and scope
WPK4FEM Jul 8, 2025
d9b40d3
refactoring: misspelled trial
WPK4FEM Jul 8, 2025
9a8b9ba
refactoring: use std::fill for averaged principal stress setting
WPK4FEM Jul 8, 2025
08a58e8
refactoring: mapping type --> averaging type, only implemented for co…
WPK4FEM Jul 8, 2025
c888d43
refactoring: the nature of AveragingType is enum behaviour.
WPK4FEM Jul 9, 2025
92bdda8
refactoring: first return mapping with constant, removed outdated com…
WPK4FEM Jul 9, 2025
fdb6071
most review comments by Gennady
WPK4FEM Jul 9, 2025
595fb9c
removed empty statement and formatted.
WPK4FEM Jul 9, 2025
1cd462f
refactoring: do not use the enum value as index, but a switch constru…
WPK4FEM Jul 9, 2025
7144edb
refactoring: static functions to anonymous namespace.
WPK4FEM Jul 9, 2025
813e5a9
refactoring: replaced if branching with switch. Review comment of Ric…
WPK4FEM Jul 10, 2025
edd24b1
Documentation is modified based on the applied averaging method
mnabideltares Jul 10, 2025
96f1c9a
pdf file of documentation is modified
mnabideltares Jul 10, 2025
c464ba5
fixes in readme
mnabideltares Jul 10, 2025
62f225d
another fix in readme
mnabideltares Jul 10, 2025
3a6c667
fix
mnabideltares Jul 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 11 additions & 6 deletions ...cations/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,10 +54,12 @@ bool IsStressAtTensionCutoffReturnZone(const Vector& rTrialSigmaTau, double Tens
rCornerPoint[1] - rTrialSigmaTau[1] - rCornerPoint[0] + rTrialSigmaTau[0] > 0.0;
}

bool IsStressAtCornerReturnZone(const Vector& rTrialSigmaTau, double DilatancyAngleInRadians, const Vector& rCornerPoint)
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;
Comment on lines +57 to 63
Copy link
Copy Markdown
Contributor

@WPK4FEM WPK4FEM Jul 7, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

}

Expand Down Expand Up @@ -112,7 +114,9 @@ bool CoulombWithTensionCutOffImpl::IsAdmissibleSigmaTau(const Vector& rTrialSigm
return coulomb_yield_function_value < coulomb_tolerance && tension_yield_function_value < tension_tolerance;
}

Vector CoulombWithTensionCutOffImpl::DoReturnMapping(const Properties& rProperties, const Vector& rTrialSigmaTau) const
Vector CoulombWithTensionCutOffImpl::DoReturnMapping(const Properties& rProperties,
const Vector& rTrialSigmaTau,
CoulombYieldSurface::CoulombAveragingType AveragingType) const
{
const auto apex = CalculateApex(ConstitutiveLawUtilities::GetFrictionAngleInRadians(rProperties),
ConstitutiveLawUtilities::GetCohesion(rProperties));
Expand All @@ -131,13 +135,14 @@ Vector CoulombWithTensionCutOffImpl::DoReturnMapping(const Properties& rProperti
}

if (IsStressAtCornerReturnZone(
rTrialSigmaTau, MathUtils<>::DegreesToRadians(rProperties[GEO_DILATANCY_ANGLE]), corner_point)) {
rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
corner_point)) {
return corner_point;
}

// Regular failure region
return ReturnStressAtRegularFailureZone(
rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau),
rTrialSigmaTau, mCoulombYieldSurface.DerivativeOfFlowFunction(rTrialSigmaTau, AveragingType),
ConstitutiveLawUtilities::GetFrictionAngleInRadians(rProperties),
ConstitutiveLawUtilities::GetCohesion(rProperties));
}
Expand Down
6 changes: 4 additions & 2 deletions applications/GeoMechanicsApplication/custom_constitutive/coulomb_with_tension_cut_off_impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,10 @@ class CoulombWithTensionCutOffImpl
double DilatancyAngleInRadians,
double TensileStrength);

[[nodiscard]] bool IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const;
[[nodiscard]] Vector DoReturnMapping(const Properties& rProperties, const Vector& rTrialSigmaTau) const;
[[nodiscard]] bool IsAdmissibleSigmaTau(const Vector& rTrialSigmaTau) const;
[[nodiscard]] Vector DoReturnMapping(const Properties& rProperties,
const Vector& rTrialSigmaTau,
CoulombYieldSurface::CoulombAveragingType AveragingType) const;

private:
CoulombYieldSurface mCoulombYieldSurface;
Expand Down
21 changes: 19 additions & 2 deletions applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,27 @@ double CoulombYieldSurface::YieldFunctionValue(const Vector& rSigmaTau) const
return rSigmaTau[1] + rSigmaTau[0] * std::sin(mFrictionAngle) - mCohesion * std::cos(mFrictionAngle);
}

Vector CoulombYieldSurface::DerivativeOfFlowFunction(const Vector&) const
Vector CoulombYieldSurface::DerivativeOfFlowFunction(const Vector& rSigmaTau) const
{
return DerivativeOfFlowFunction(rSigmaTau, CoulombAveragingType::NO_AVERAGING);
}

Vector CoulombYieldSurface::DerivativeOfFlowFunction(const Vector&, CoulombAveragingType AveragingType) const
{
Vector result(2);
result <<= std::sin(mDilatationAngle), 1.0;
switch (AveragingType) {
case CoulombAveragingType::LOWEST_PRINCIPAL_STRESSES:
result <<= -(1.0 - 3.0 * std::sin(mDilatationAngle)) / 4.0, (3.0 - std::sin(mDilatationAngle)) / 4.0;
break;
case CoulombAveragingType::NO_AVERAGING:
result <<= std::sin(mDilatationAngle), 1.0;
break;
case CoulombAveragingType::HIGHEST_PRINCIPAL_STRESSES:
result <<= (1.0 + 3.0 * std::sin(mDilatationAngle)) / 4.0, (3.0 + std::sin(mDilatationAngle)) / 4.0;
break;
default:
KRATOS_ERROR << "Unsupported Averaging Type: " << static_cast<std::size_t>(AveragingType) << "\n";
}
return result;
}

Expand Down
7 changes: 7 additions & 0 deletions applications/GeoMechanicsApplication/custom_constitutive/coulomb_yield_surface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,19 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) CoulombYieldSurface : public YieldSu
public:
KRATOS_CLASS_POINTER_DEFINITION(CoulombYieldSurface);

enum class CoulombAveragingType {
NO_AVERAGING,
LOWEST_PRINCIPAL_STRESSES,
HIGHEST_PRINCIPAL_STRESSES
};

CoulombYieldSurface() = default;

CoulombYieldSurface(double FrictionAngleInRad, double Cohesion, double DilatationAngleInRad);

[[nodiscard]] double YieldFunctionValue(const Vector& rSigmaTau) const override;
[[nodiscard]] Vector DerivativeOfFlowFunction(const Vector&) const override;
[[nodiscard]] Vector DerivativeOfFlowFunction(const Vector&, CoulombAveragingType AveragingType) const;

private:
friend class Serializer;
Expand Down
10 changes: 6 additions & 4 deletions ...ns/GeoMechanicsApplication/custom_constitutive/interface_coulomb_with_tension_cut_off.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,16 +141,18 @@ void InterfaceCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(Paramete
auto trial_sigma_tau = CalculateTrialTractionVector(rConstitutiveLawParameters.GetStrainVector(),
r_properties[INTERFACE_NORMAL_STIFFNESS],
r_properties[INTERFACE_SHEAR_STIFFNESS]);
auto mapped_sigma_tau = trial_sigma_tau;

const auto negative = std::signbit(trial_sigma_tau[1]);
trial_sigma_tau[1] = std::abs(trial_sigma_tau[1]);

if (!mCoulombWithTensionCutOffImpl.IsAdmissibleSigmaTau(trial_sigma_tau)) {
trial_sigma_tau = mCoulombWithTensionCutOffImpl.DoReturnMapping(r_properties, trial_sigma_tau);
mapped_sigma_tau = mCoulombWithTensionCutOffImpl.DoReturnMapping(
r_properties, trial_sigma_tau, CoulombYieldSurface::CoulombAveragingType::NO_AVERAGING);
if (negative) mapped_sigma_tau[1] *= -1.0;
}

if (negative) trial_sigma_tau[1] *= -1.0;

mTractionVector = trial_sigma_tau;
mTractionVector = mapped_sigma_tau;
rConstitutiveLawParameters.GetStressVector() = mTractionVector;
}

Expand Down
82 changes: 70 additions & 12 deletions ...ications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,49 @@

#include <cmath>

namespace
{
using namespace Kratos;

std::size_t AveragingTypeToArrayIndex(CoulombYieldSurface::CoulombAveragingType AveragingType)
{
switch (AveragingType) {
case CoulombYieldSurface::CoulombAveragingType::LOWEST_PRINCIPAL_STRESSES:
return 0;
case CoulombYieldSurface::CoulombAveragingType::HIGHEST_PRINCIPAL_STRESSES:
return 2;
default:
return 1;
}
}

Vector AveragePrincipalStressComponents(const Vector& rPrincipalStressVector,
CoulombYieldSurface::CoulombAveragingType AveragingType)
{
auto result = rPrincipalStressVector;
if (AveragingType == CoulombYieldSurface::CoulombAveragingType::LOWEST_PRINCIPAL_STRESSES) {
std::fill(result.begin(), result.begin() + 1,
(rPrincipalStressVector[0] + rPrincipalStressVector[1]) * 0.5);
} else if (AveragingType == CoulombYieldSurface::CoulombAveragingType::HIGHEST_PRINCIPAL_STRESSES) {
std::fill(result.begin() + 1, result.begin() + 2,
(rPrincipalStressVector[1] + rPrincipalStressVector[2]) * 0.5);
}
return result;
}

CoulombYieldSurface::CoulombAveragingType FindAveragingType(const Vector& rMappedPrincipalStressVector)
{
if (rMappedPrincipalStressVector[0] < rMappedPrincipalStressVector[1]) {
return CoulombYieldSurface::CoulombAveragingType::LOWEST_PRINCIPAL_STRESSES;
}
if (rMappedPrincipalStressVector[1] < rMappedPrincipalStressVector[2]) {
return CoulombYieldSurface::CoulombAveragingType::HIGHEST_PRINCIPAL_STRESSES;
}
return CoulombYieldSurface::CoulombAveragingType::NO_AVERAGING;
}

} // namespace

namespace Kratos
{

Expand Down Expand Up @@ -149,28 +192,43 @@ void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveL
mpConstitutiveDimension->CalculateElasticMatrix(r_prop[YOUNG_MODULUS], r_prop[POISSON_RATIO]);
}

const auto trail_stress_vector = CalculateTrialStressVector(
const auto trial_stress_vector = CalculateTrialStressVector(
rParameters.GetStrainVector(), r_prop[YOUNG_MODULUS], r_prop[POISSON_RATIO]);

Vector principal_trial_stress_vector;
Matrix rotation_matrix;
StressStrainUtilities::CalculatePrincipalStresses(
trail_stress_vector, principal_trial_stress_vector, rotation_matrix);
trial_stress_vector, principal_trial_stress_vector, rotation_matrix);

auto trial_sigma_tau =
StressStrainUtilities::TransformPrincipalStressesToSigmaTau(principal_trial_stress_vector);
while (!mCoulombWithTensionCutOffImpl.IsAdmissibleSigmaTau(trial_sigma_tau)) {
trial_sigma_tau = mCoulombWithTensionCutOffImpl.DoReturnMapping(r_prop, trial_sigma_tau);
principal_trial_stress_vector = StressStrainUtilities::TransformSigmaTauToPrincipalStresses(
trial_sigma_tau, principal_trial_stress_vector);

StressStrainUtilities::ReorderEigenValuesAndVectors(principal_trial_stress_vector, rotation_matrix);
trial_sigma_tau =
StressStrainUtilities::TransformPrincipalStressesToSigmaTau(principal_trial_stress_vector);

if (mCoulombWithTensionCutOffImpl.IsAdmissibleSigmaTau(trial_sigma_tau)) {
mStressVector = trial_stress_vector;
} else {
auto mapped_sigma_tau = mCoulombWithTensionCutOffImpl.DoReturnMapping(
r_prop, trial_sigma_tau, CoulombYieldSurface::CoulombAveragingType::NO_AVERAGING);
auto mapped_principal_stress_vector = StressStrainUtilities::TransformSigmaTauToPrincipalStresses(
mapped_sigma_tau, principal_trial_stress_vector);

// for interchanging principal stresses, retry mapping with averaged principal stresses.
const auto averaging_type = FindAveragingType(mapped_principal_stress_vector);
if (averaging_type != CoulombYieldSurface::CoulombAveragingType::NO_AVERAGING) {
const auto averaged_principal_trial_stress_vector =
AveragePrincipalStressComponents(principal_trial_stress_vector, averaging_type);
trial_sigma_tau = StressStrainUtilities::TransformPrincipalStressesToSigmaTau(
averaged_principal_trial_stress_vector);
mapped_sigma_tau =
mCoulombWithTensionCutOffImpl.DoReturnMapping(r_prop, trial_sigma_tau, averaging_type);
mapped_principal_stress_vector = StressStrainUtilities::TransformSigmaTauToPrincipalStresses(
Copy link
Copy Markdown
Contributor

@rfaasse rfaasse Jul 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor

@WPK4FEM WPK4FEM Jul 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor Author

@mnabideltares mnabideltares Jul 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

mapped_sigma_tau, averaged_principal_trial_stress_vector);
mapped_principal_stress_vector[1] =
mapped_principal_stress_vector[AveragingTypeToArrayIndex(averaging_type)];
}
mStressVector = StressStrainUtilities::RotatePrincipalStresses(
mapped_principal_stress_vector, rotation_matrix, mpConstitutiveDimension->GetStrainSize());
}

mStressVector = StressStrainUtilities::RotatePrincipalStresses(
principal_trial_stress_vector, rotation_matrix, mpConstitutiveDimension->GetStrainSize());
rParameters.GetStressVector() = mStressVector;
}

Expand Down
1 change: 1 addition & 0 deletions applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) MohrCoulombWithTensionCutOff : publi
[[nodiscard]] Vector CalculateTrialStressVector(const Vector& rStrainVector,
double YoungsModulus,
double PoissonsRatio) const;

friend class Serializer;
void save(Serializer& rSerializer) const override;
void load(Serializer& rSerializer) override;
Expand Down
22 changes: 14 additions & 8 deletions ...Application/tests/cpp_tests/custom_constitutive/test_mohr_coulomb_with_tension_cutoff.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,9 +203,10 @@ KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponse
expected_cauchy_stress_vector, Defaults::absolute_tolerance);

cauchy_stress_vector <<= 12.0, 10.0, -16.0, 0.0;
expected_cauchy_stress_vector <<= 7.338673315592010089, 7.90609951999166506, -9.24477283558367515, 0.0;
expected_cauchy_stress_vector <<= 7.806379130008, 7.806379130008, -9.61275826001616129, 0.0;
constexpr auto tolerance = 2.0e-10;
KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
expected_cauchy_stress_vector, Defaults::absolute_tolerance);
expected_cauchy_stress_vector, tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponseCauchyAtCornerReturnZone,
Expand Down Expand Up @@ -239,8 +240,9 @@ KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponse

cauchy_stress_vector <<= 24.0, 22.0, -8.0, 0.0;
expected_cauchy_stress_vector <<= 10.0, 10.0, -1.5179192179966735, 0.0;
constexpr double tolerance = 1.0e-10;
KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
expected_cauchy_stress_vector, Defaults::absolute_tolerance);
expected_cauchy_stress_vector, tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponseCauchyAtTensionCutoffReturnZone,
Expand Down Expand Up @@ -274,8 +276,9 @@ KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponse

cauchy_stress_vector <<= 14.0, 12.0, 6.0, 0.0;
expected_cauchy_stress_vector <<= 10.0, 10.0, 6.0, 0.0;
constexpr double tolerance = 1.0e-10;
KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
expected_cauchy_stress_vector, Defaults::absolute_tolerance);
expected_cauchy_stress_vector, tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponseCauchyAtTensionApexReturnZone,
Expand All @@ -293,19 +296,20 @@ KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponse
const auto dummy_element_geometry = Geometry<Node>{};
const auto dummy_shape_function_values = Vector{};
law.InitializeMaterial(properties, dummy_element_geometry, dummy_shape_function_values);
constexpr double tolerance = 1.0e-10;

// Act and Assert
Vector cauchy_stress_vector(4);
cauchy_stress_vector <<= 19.0, 12.0, 11.0, 0.0;
Vector expected_cauchy_stress_vector(4);
expected_cauchy_stress_vector <<= 10.0, 10.0, 10.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
expected_cauchy_stress_vector, Defaults::absolute_tolerance);
expected_cauchy_stress_vector, tolerance);

cauchy_stress_vector <<= 11.5, 10.0, 10.5, 0.0;
expected_cauchy_stress_vector <<= 10.0, 10.0, 10.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
expected_cauchy_stress_vector, Defaults::absolute_tolerance);
expected_cauchy_stress_vector, tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponseCauchyWithLargeTensileStrength,
Expand All @@ -329,8 +333,9 @@ KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponse
cauchy_stress_vector <<= 22.0, 20.0, 18.0, 0.0;
Vector expected_cauchy_stress_vector(4);
expected_cauchy_stress_vector <<= 14.2814800674211450, 14.2814800674211450, 14.2814800674211450, 0.0;
constexpr double tolerance = 1.0e-10;
KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
expected_cauchy_stress_vector, Defaults::absolute_tolerance);
expected_cauchy_stress_vector, tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_Serialization, KratosGeoMechanicsFastSuiteWithoutKernel)
Expand Down Expand Up @@ -627,8 +632,9 @@ KRATOS_TEST_CASE_IN_SUITE(MohrCoulombWithTensionCutOff_CalculateMaterialResponse
cauchy_stress_vector <<= 10.0, 0.0, 0.0, -2.0;
Vector expected_cauchy_stress_vector(4);
expected_cauchy_stress_vector <<= -1.0, -1.0, -1.0, 0.0;
constexpr double tolerance = 1.0e-10;
KRATOS_EXPECT_VECTOR_NEAR(CalculateMappedStressVector(cauchy_stress_vector, parameters, law),
expected_cauchy_stress_vector, Defaults::absolute_tolerance);
expected_cauchy_stress_vector, tolerance);
}

} // namespace Kratos::Testing
4 changes: 2 additions & 2 deletions applications/GeoMechanicsApplication/tests/test_mohr_coulomb_with_tension_cutoff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,8 @@ def test_dirichlet_tension_cutoff_return_zone_3d(self):

def test_column_under_gravity(self):
sig_xx, sig_yy = self.simulate_mohr_coulomb('test_column_under_gravity', 2)
self.assertAlmostEqual(sig_xx, -8400.786160492225)
self.assertAlmostEqual(sig_yy, -22365.87047624526)
self.assertAlmostEqual(sig_xx, -6300.238764425369)
self.assertAlmostEqual(sig_yy, -22364.817908413854)

def test_interface_coulomb_2plus2(self):
test_name = "interface_coulomb_2plus2"
Expand Down
Loading