Skip to content

[GeoMechanicsApplication] Extend ApplyCPhiReductionProcess to use Mohr-Coulomb constitutive model #14085

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

Open
markelov208 wants to merge 1 commit into
base: master
Choose a base branch
from
+58 −44
Open

[GeoMechanicsApplication] Extend ApplyCPhiReductionProcess to use Mohr-Coulomb constitutive model #14085

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
56 changes: 31 additions & 25 deletions applications/GeoMechanicsApplication/custom_processes/apply_c_phi_reduction_process.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,24 @@ int ApplyCPhiReductionProcess::Check()
KRATOS_ERROR_IF(std::ranges::all_of(mrModelParts, [](const auto& r_model_part) {
return r_model_part.get().Elements().empty();
})) << "None of the provided model parts contains at least one element. A c-phi reduction analysis requires at least one element.\n";

for (const auto& r_model_part : mrModelParts) {
for (const auto& r_element : r_model_part.get().Elements()) {
auto r_properties = r_element.GetProperties();
const CheckProperties check_properties(r_properties, "model part property",
CheckProperties::Bounds::AllInclusive);
if (r_properties.Has(GEO_COHESION)) {
check_properties.Check(GEO_COHESION);
check_properties.Check(GEO_FRICTION_ANGLE);
} else {
check_properties.CheckAvailability(UMAT_PARAMETERS);
check_properties.Check(INDEX_OF_UMAT_PHI_PARAMETER, 1,
static_cast<int>(r_properties[UMAT_PARAMETERS].size()));
check_properties.Check(INDEX_OF_UMAT_C_PARAMETER, 1,
static_cast<int>(r_properties[UMAT_PARAMETERS].size()));
}
}
}
return 0;
}

Expand All @@ -99,12 +117,6 @@ double ApplyCPhiReductionProcess::GetAndCheckPhi(const Properties& rModelPartPro
// properties object with reduced c and phi for each and every element. Those reduced
// properties objects are not linked to the original ones.

const CheckProperties check_properties(rModelPartProperties, "model part property",
CheckProperties::Bounds::AllInclusive);
check_properties.CheckAvailability(UMAT_PARAMETERS);
check_properties.Check(INDEX_OF_UMAT_PHI_PARAMETER, 1,
static_cast<int>(rModelPartProperties[UMAT_PARAMETERS].size()));

const auto phi = ConstitutiveLawUtilities::GetFrictionAngleInDegrees(rModelPartProperties);
KRATOS_ERROR_IF(phi < 0. || phi > 90.)
<< "Friction angle Phi in the model part property with Id " << ElementPropertyId
Expand All @@ -126,32 +138,26 @@ double ApplyCPhiReductionProcess::GetAndCheckC(const Properties& rModelPartPrope
// properties object with reduced c and phi for each and every element. Those reduced
// properties objects are not linked to the original ones.

KRATOS_ERROR_IF_NOT(rModelPartProperties.Has(UMAT_PARAMETERS))
<< "Missing required item UMAT_PARAMETERS" << std::endl;
KRATOS_ERROR_IF_NOT(rModelPartProperties.Has(INDEX_OF_UMAT_C_PARAMETER))
<< "Missing required item INDEX_OF_UMAT_C_PARAMETER" << std::endl;

KRATOS_ERROR_IF(rModelPartProperties[INDEX_OF_UMAT_C_PARAMETER] < 1 ||
rModelPartProperties[INDEX_OF_UMAT_C_PARAMETER] >
static_cast<int>(rModelPartProperties[UMAT_PARAMETERS].size()))
<< "invalid INDEX_OF_UMAT_C_PARAMETER: " << rModelPartProperties[INDEX_OF_UMAT_C_PARAMETER]
<< " (out-of-bounds index)" << std::endl;
const auto c = ConstitutiveLawUtilities::GetCohesion(rModelPartProperties);
KRATOS_ERROR_IF(c < 0.) << "Cohesion C out of range: " << c << std::endl;
return c;
}

void ApplyCPhiReductionProcess::SetCPhiAtElement(Element& rElement, double ReducedPhi, double ReducedC)
{
// Get C/Phi material properties of this element
const auto& r_prop = rElement.GetProperties();

// Overwrite C and Phi in the UMAT_PARAMETERS
auto Umat_parameters = r_prop[UMAT_PARAMETERS];
Umat_parameters[r_prop[INDEX_OF_UMAT_PHI_PARAMETER] - 1] = ReducedPhi;
Umat_parameters[r_prop[INDEX_OF_UMAT_C_PARAMETER] - 1] = ReducedC;

SetValueAtElement(rElement, UMAT_PARAMETERS, Umat_parameters);
auto& r_properties = rElement.GetProperties();

if (r_properties.Has(GEO_FRICTION_ANGLE)) {
r_properties.SetValue(GEO_FRICTION_ANGLE, ReducedPhi);
r_properties.SetValue(GEO_COHESION, ReducedC);
} else {
// Overwrite C and Phi in the UMAT_PARAMETERS
auto Umat_parameters = r_properties[UMAT_PARAMETERS];
Umat_parameters[r_properties[INDEX_OF_UMAT_PHI_PARAMETER] - 1] = ReducedPhi;
Umat_parameters[r_properties[INDEX_OF_UMAT_C_PARAMETER] - 1] = ReducedC;

SetValueAtElement(rElement, UMAT_PARAMETERS, Umat_parameters);
}
}

void ApplyCPhiReductionProcess::SetValueAtElement(Element& rElement, const Variable<Vector>& rVariable, const Vector& rValue)
Expand Down
46 changes: 27 additions & 19 deletions ...hanicsApplication/tests/cpp_tests/custom_processes/test_apply_c_phi_reduction_process.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,25 +115,46 @@ KRATOS_TEST_CASE_IN_SUITE(CheckFailureUmatInputsApplyCPhiReductionProcess, Krato
const auto parameters = Parameters{R"({"model_part_name" : "dummy"})"};

KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()),
(ApplyCPhiReductionProcess{model, parameters}.Check()),
"UMAT_PARAMETERS does not exist in the model part property with Id 0.")

Vector umat_parameters(6);
umat_parameters <<= 10000000, 0.2, 10.0, 25.0, 25.0, 1000;
r_model_part_properties.SetValue(UMAT_PARAMETERS, umat_parameters);

// checking of Phi
// checking settings for Phi
KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()),
(ApplyCPhiReductionProcess{model, parameters}.Check()),
"INDEX_OF_UMAT_PHI_PARAMETER does not exist in the model part property with Id 0.")

r_model_part_properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 0);
KRATOS_EXPECT_EXCEPTION_IS_THROWN((ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()), "INDEX_OF_UMAT_PHI_PARAMETER in the model part property with Id 0 has an invalid value: 0 is out of the range [1, 6].")
KRATOS_EXPECT_EXCEPTION_IS_THROWN((ApplyCPhiReductionProcess{model, parameters}.Check()),
"INDEX_OF_UMAT_PHI_PARAMETER in the model part property with "
"Id 0 has an invalid value: 0 is out of the range [1, 6].")

r_model_part_properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 7);
KRATOS_EXPECT_EXCEPTION_IS_THROWN((ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()), "INDEX_OF_UMAT_PHI_PARAMETER in the model part property with Id 0 has an invalid value: 7 is out of the range [1, 6].")
KRATOS_EXPECT_EXCEPTION_IS_THROWN((ApplyCPhiReductionProcess{model, parameters}.Check()),
"INDEX_OF_UMAT_PHI_PARAMETER in the model part property with "
"Id 0 has an invalid value: 7 is out of the range [1, 6].")

r_model_part_properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 4);

// checking settings for c
KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.Check()),
"INDEX_OF_UMAT_C_PARAMETER does not exist in the model part property with Id 0.")

r_model_part_properties.SetValue(INDEX_OF_UMAT_C_PARAMETER, 0);
KRATOS_EXPECT_EXCEPTION_IS_THROWN((ApplyCPhiReductionProcess{model, parameters}.Check()),
"INDEX_OF_UMAT_C_PARAMETER in the model part property with "
"Id 0 has an invalid value: 0 is out of the range [1, 6].")

r_model_part_properties.SetValue(INDEX_OF_UMAT_C_PARAMETER, 7);
KRATOS_EXPECT_EXCEPTION_IS_THROWN((ApplyCPhiReductionProcess{model, parameters}.Check()),
"INDEX_OF_UMAT_C_PARAMETER in the model part property with "
"Id 0 has an invalid value: 7 is out of the range [1, 6].")

// checking Phi value
umat_parameters(3) = -0.0001;
r_model_part_properties.SetValue(UMAT_PARAMETERS, umat_parameters);

Expand All @@ -145,21 +166,8 @@ KRATOS_TEST_CASE_IN_SUITE(CheckFailureUmatInputsApplyCPhiReductionProcess, Krato

umat_parameters(3) = 25.0;
r_model_part_properties.SetValue(UMAT_PARAMETERS, umat_parameters);
// checking of c
KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()),
"Missing required item INDEX_OF_UMAT_C_PARAMETER")

r_model_part_properties.SetValue(INDEX_OF_UMAT_C_PARAMETER, 0);
KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()),
"invalid INDEX_OF_UMAT_C_PARAMETER: 0 (out-of-bounds index)")

r_model_part_properties.SetValue(INDEX_OF_UMAT_C_PARAMETER, 7);
KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()),
"invalid INDEX_OF_UMAT_C_PARAMETER: 7 (out-of-bounds index)")

// checking c value
r_model_part_properties.SetValue(INDEX_OF_UMAT_C_PARAMETER, 3);
umat_parameters(2) = -0.00001;
r_model_part_properties.SetValue(UMAT_PARAMETERS, umat_parameters);
Expand Down
Loading