diff --git a/applications/GeoMechanicsApplication/custom_processes/README.md b/applications/GeoMechanicsApplication/custom_processes/README.md index 98e8835105a0..ed53e26d83ab 100644 --- a/applications/GeoMechanicsApplication/custom_processes/README.md +++ b/applications/GeoMechanicsApplication/custom_processes/README.md @@ -85,7 +85,7 @@ When set to true, the material used for the modelpart used for the $K_0$ procedu When the stress computation is completed, the original constitutive law is restored. Depending on the given input parameters, the following scheme is adapted for computation of the $K_0$ value. -$K_0^{nc}$ is gotten from either "K0_NC" the material input file or by computation from input of "INDEX_OF_UMAT_PHI_PARAMETER" and "UMAT_PARAMETERS": +$K_0^{nc}$ is gotten from either "K0_NC" the material input file or by computation from input of "GEO_FRICTION_ANGLE" or "INDEX_OF_UMAT_PHI_PARAMETER" and "UMAT_PARAMETERS": $$K_0^{nc} = 1.0 - \sin \phi$$ diff --git a/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.cpp b/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.cpp index 8a3b12fc8c3b..575b01cc37eb 100644 --- a/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.cpp +++ b/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.cpp @@ -13,7 +13,6 @@ #include "apply_k0_procedure_process.h" #include -#include #include "containers/model.h" #include "custom_constitutive/linear_elastic_law.h" @@ -32,9 +31,10 @@ using namespace std::string_literals; void SetConsiderDiagonalEntriesOnlyAndNoShear(ModelPart::ElementsContainerType& rElements, bool Whether) { block_for_each(rElements, [Whether](Element& rElement) { - auto pLinearElasticLaw = + const auto p_linear_elastic_law = dynamic_cast(rElement.GetProperties().GetValue(CONSTITUTIVE_LAW).get()); - if (pLinearElasticLaw) pLinearElasticLaw->SetConsiderDiagonalEntriesOnlyAndNoShear(Whether); + if (p_linear_elastic_law) + p_linear_elastic_law->SetConsiderDiagonalEntriesOnlyAndNoShear(Whether); }); } @@ -72,7 +72,7 @@ int ApplyK0ProcedureProcess::Check() { for (const auto& r_model_part : mrModelParts) { KRATOS_ERROR_IF(r_model_part.get().Elements().empty()) - << "ApplyK0ProcedureProces has no elements in modelpart " << r_model_part.get().Name() + << "ApplyK0ProcedureProcess has no elements in modelpart " << r_model_part.get().Name() << std::endl; } @@ -83,8 +83,10 @@ int ApplyK0ProcedureProcess::Check() CheckSufficientMaterialParameters(r_properties, rElement.Id()); CheckOCRorPOP(r_properties, rElement.Id()); CheckPoissonUnloadingReloading(r_properties, rElement.Id()); - CheckPhi(r_properties, rElement.Id()); CheckK0(r_properties, rElement.Id()); + if (!(r_properties.Has(K0_NC) || (r_properties.Has(K0_VALUE_XX) && r_properties.Has(K0_VALUE_YY) && + r_properties.Has(K0_VALUE_ZZ)))) + ConstitutiveLawUtilities::ValidateFrictionAngle(r_properties, rElement.Id()); }); } @@ -97,15 +99,19 @@ void ApplyK0ProcedureProcess::CheckK0MainDirection(const Properties& rProperties << "K0_MAIN_DIRECTION is not defined for element " << ElementId << "." << std::endl; const auto dimension = rProperties.GetValue(CONSTITUTIVE_LAW).get()->WorkingSpaceDimension(); - if (dimension == 2) { - KRATOS_ERROR_IF(rProperties[K0_MAIN_DIRECTION] < 0 || rProperties[K0_MAIN_DIRECTION] > 1) - << "K0_MAIN_DIRECTION should be 0 or 1 for element " << ElementId << "." << std::endl; - } else if (dimension == 3) { - KRATOS_ERROR_IF(rProperties[K0_MAIN_DIRECTION] < 0 || rProperties[K0_MAIN_DIRECTION] > 2) - << "K0_MAIN_DIRECTION should be 0, 1 or 2 for element " << ElementId << "." << std::endl; - } else { - KRATOS_ERROR << "dimension should be 2 or 3 for element " << ElementId << "." << std::endl; + KRATOS_ERROR_IF(dimension != 2 && dimension != 3) + << "dimension should be 2 or 3 for element " << ElementId << "." << std::endl; + std::ostringstream valid_values; + for (std::size_t i = 0; i < dimension; ++i) { + if (i > 0) { + if (i == dimension - 1) valid_values << " or "; + else valid_values << ", "; + } + valid_values << i; } + KRATOS_ERROR_IF(rProperties[K0_MAIN_DIRECTION] < 0 || rProperties[K0_MAIN_DIRECTION] > static_cast(dimension)) + << "K0_MAIN_DIRECTION should be " << valid_values.str() << " for element " << ElementId + << "." << std::endl; } void ApplyK0ProcedureProcess::CheckK0(const Properties& rProperties, IndexType ElementId) @@ -122,37 +128,21 @@ void ApplyK0ProcedureProcess::CheckK0(const Properties& rProperties, IndexType E } } -void ApplyK0ProcedureProcess::CheckPhi(const Properties& rProperties, IndexType ElementId) -{ - if (rProperties.Has(INDEX_OF_UMAT_PHI_PARAMETER) && rProperties.Has(UMAT_PARAMETERS)) { - const auto phi_index = rProperties[INDEX_OF_UMAT_PHI_PARAMETER]; - const auto number_of_umat_parameters = static_cast(rProperties[UMAT_PARAMETERS].size()); - - KRATOS_ERROR_IF(phi_index < 1 || phi_index > number_of_umat_parameters) - << "INDEX_OF_UMAT_PHI_PARAMETER (" << phi_index - << ") is not in range 1, size of UMAT_PARAMETERS for element " << ElementId << "." << std::endl; - - const double phi = rProperties[UMAT_PARAMETERS][phi_index - 1]; - KRATOS_ERROR_IF(phi < 0.0 || phi > 90.0) - << "Phi (" << phi << ") should be between 0 and 90 degrees for element " << ElementId - << "." << std::endl; - } -} - void ApplyK0ProcedureProcess::CheckOCRorPOP(const Properties& rProperties, IndexType ElementId) { - if (rProperties.Has(K0_NC) || - (rProperties.Has(INDEX_OF_UMAT_PHI_PARAMETER) && rProperties.Has(UMAT_PARAMETERS))) { + if (rProperties.Has(K0_NC) || ConstitutiveLawUtilities::HasFrictionAngle(rProperties)) { if (rProperties.Has(OCR)) { const auto ocr = rProperties[OCR]; - KRATOS_ERROR_IF(ocr < 1.0) << "OCR (" << ocr << ") should be in the range [1.0,-> for element " - << ElementId << "." << std::endl; + KRATOS_ERROR_IF(ocr < 1.0) + << "OCR (" << ocr << ") should be in the range [1.0,-> for property Id of " + << rProperties.Id() << " for element " << ElementId << "." << std::endl; } if (rProperties.Has(POP)) { const auto pop = rProperties[POP]; - KRATOS_ERROR_IF(pop < 0.0) << "POP (" << pop << ") should be in the range [0.0,-> for element " - << ElementId << "." << std::endl; + KRATOS_ERROR_IF(pop < 0.0) + << "POP (" << pop << ") should be in the range [0.0,-> for property Id of " + << rProperties.Id() << " element " << ElementId << "." << std::endl; } } } @@ -163,12 +153,14 @@ void ApplyK0ProcedureProcess::CheckPoissonUnloadingReloading(const Properties& r (rProperties[POISSON_UNLOADING_RELOADING] < -1.0 || rProperties[POISSON_UNLOADING_RELOADING] >= 0.5)) << "POISSON_UNLOADING_RELOADING (" << rProperties[POISSON_UNLOADING_RELOADING] - << ") is not in range [-1.0, 0.5> for element " << ElementId << "." << std::endl; + << ") is not in range [-1.0, 0.5> for property Id of " << rProperties.Id() + << " for element " << ElementId << "." << std::endl; if (rProperties.Has(K0_VALUE_XX)) { KRATOS_ERROR_IF(rProperties.Has(POISSON_UNLOADING_RELOADING) || rProperties.Has(OCR) || rProperties.Has(POP)) - << "Insufficient material data for K0 procedure process for element " + << "Insufficient material data for K0 procedure process for property Id of " + << rProperties.Id() << " for element " << ElementId << ". Poisson unloading-reloading, OCR and POP functionality cannot be combined with K0_VALUE_XX, _YY and _ZZ." << std::endl; } @@ -177,10 +169,10 @@ void ApplyK0ProcedureProcess::CheckPoissonUnloadingReloading(const Properties& r void ApplyK0ProcedureProcess::CheckSufficientMaterialParameters(const Properties& rProperties, IndexType ElementId) { KRATOS_ERROR_IF_NOT( - rProperties.Has(K0_NC) || - (rProperties.Has(INDEX_OF_UMAT_PHI_PARAMETER) && rProperties.Has(UMAT_PARAMETERS)) || + rProperties.Has(K0_NC) || ConstitutiveLawUtilities::HasFrictionAngle(rProperties) || (rProperties.Has(K0_VALUE_XX) && rProperties.Has(K0_VALUE_YY) && rProperties.Has(K0_VALUE_ZZ))) - << "Insufficient material data for K0 procedure process for element " << ElementId << ". No K0_NC, " + << "Insufficient material data for K0 procedure process for property Id of " + << rProperties.Id() << " for element " << ElementId << ". No K0_NC, " << "(INDEX_OF_UMAT_PHI_PARAMETER and UMAT_PARAMETERS) or (K0_VALUE_XX, _YY and _ZZ found)." << std::endl; } @@ -206,19 +198,19 @@ bool ApplyK0ProcedureProcess::UseStandardProcedure() const return !mSettings.Has(setting_name) || mSettings[setting_name].GetBool(); } -array_1d ApplyK0ProcedureProcess::CreateK0Vector(const Element::PropertiesType& rProp) +array_1d ApplyK0ProcedureProcess::CreateK0Vector(const Element::PropertiesType& rProperties) { // Check for alternative K0 specifications array_1d k0_vector; - if (rProp.Has(K0_NC)) { - std::fill(k0_vector.begin(), k0_vector.end(), rProp[K0_NC]); - } else if (rProp.Has(INDEX_OF_UMAT_PHI_PARAMETER) && rProp.Has(UMAT_PARAMETERS)) { + if (rProperties.Has(K0_NC)) { + std::ranges::fill(k0_vector, rProperties[K0_NC]); + } else if (ConstitutiveLawUtilities::HasFrictionAngle(rProperties)) { std::ranges::fill(k0_vector, ConstitutiveLawUtilities::CalculateK0NCFromFrictionAngleInRadians( - ConstitutiveLawUtilities::GetFrictionAngleInRadians(rProp))); + ConstitutiveLawUtilities::GetFrictionAngleInRadians(rProperties))); } else { - k0_vector[0] = rProp[K0_VALUE_XX]; - k0_vector[1] = rProp[K0_VALUE_YY]; - k0_vector[2] = rProp[K0_VALUE_ZZ]; + k0_vector[0] = rProperties[K0_VALUE_XX]; + k0_vector[1] = rProperties[K0_VALUE_YY]; + k0_vector[2] = rProperties[K0_VALUE_ZZ]; } return k0_vector; @@ -227,45 +219,46 @@ array_1d ApplyK0ProcedureProcess::CreateK0Vector(const Element::Prope void ApplyK0ProcedureProcess::CalculateK0Stresses(Element& rElement, const ProcessInfo& rProcessInfo) { // Get K0 material parameters of this element ( probably there is something more efficient ) - const Element::PropertiesType& rProp = rElement.GetProperties(); - const auto k0_main_direction = rProp[K0_MAIN_DIRECTION]; + const Element::PropertiesType& r_properties = rElement.GetProperties(); + const auto k0_main_direction = r_properties[K0_MAIN_DIRECTION]; - auto k0_vector = CreateK0Vector(rProp); + auto k0_vector = CreateK0Vector(r_properties); // Corrections on k0_vector by OCR or POP - const auto PoissonUR = rProp.Has(POISSON_UNLOADING_RELOADING) ? rProp[POISSON_UNLOADING_RELOADING] : 0.; - const auto PoissonURfactor = PoissonUR / (1. - PoissonUR); + const auto poisson_u_r = + r_properties.Has(POISSON_UNLOADING_RELOADING) ? r_properties[POISSON_UNLOADING_RELOADING] : 0.; + const auto poisson_u_r_factor = poisson_u_r / (1. - poisson_u_r); double POP_value = 0.0; - if (rProp.Has(K0_NC) || rProp.Has(INDEX_OF_UMAT_PHI_PARAMETER)) { - if (rProp.Has(OCR)) { + if (r_properties.Has(K0_NC) || ConstitutiveLawUtilities::HasFrictionAngle(r_properties)) { + if (r_properties.Has(OCR)) { // Determine OCR dependent K0 values ( constant per element! ) - k0_vector *= rProp[OCR]; - const array_1d correction(3, PoissonURfactor * (rProp[OCR] - 1.0)); + k0_vector *= r_properties[OCR]; + const array_1d correction(3, poisson_u_r_factor * (r_properties[OCR] - 1.0)); k0_vector -= correction; - } else if (rProp.Has(POP)) { + } else if (r_properties.Has(POP)) { // POP is entered as positive value, convention here is compression negative. - POP_value = -rProp[POP]; + POP_value = -r_properties[POP]; } } // Get element stress vectors - std::vector rStressVectors; - rElement.CalculateOnIntegrationPoints(CAUCHY_STRESS_VECTOR, rStressVectors, rProcessInfo); + std::vector stress_vectors; + rElement.CalculateOnIntegrationPoints(CAUCHY_STRESS_VECTOR, stress_vectors, rProcessInfo); // Loop over integration point stress vectors - for (auto& rStressVector : rStressVectors) { + for (auto& r_stress_vector : stress_vectors) { // Apply K0 procedure for (int i_dir = 0; i_dir <= 2; ++i_dir) { if (i_dir != k0_main_direction) { - rStressVector[i_dir] = k0_vector[i_dir] * (rStressVector[k0_main_direction] + POP_value) - - PoissonURfactor * POP_value; + r_stress_vector[i_dir] = k0_vector[i_dir] * (r_stress_vector[k0_main_direction] + POP_value) - + poisson_u_r_factor * POP_value; } } // Erase shear stresses - std::fill(rStressVector.begin() + 3, rStressVector.end(), 0.0); + std::fill(r_stress_vector.begin() + 3, r_stress_vector.end(), 0.0); } // Set element integration point stress tensors - rElement.SetValuesOnIntegrationPoints(CAUCHY_STRESS_VECTOR, rStressVectors, rProcessInfo); + rElement.SetValuesOnIntegrationPoints(CAUCHY_STRESS_VECTOR, stress_vectors, rProcessInfo); } } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.h b/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.h index 036402216da5..1729852b65ad 100644 --- a/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.h +++ b/applications/GeoMechanicsApplication/custom_processes/apply_k0_procedure_process.h @@ -40,13 +40,12 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) ApplyK0ProcedureProcess : public Pro [[nodiscard]] std::string Info() const override; private: - [[nodiscard]] bool UseStandardProcedure() const; - [[nodiscard]] static array_1d CreateK0Vector(const Element::PropertiesType& rProp); + [[nodiscard]] bool UseStandardProcedure() const; + [[nodiscard]] static array_1d CreateK0Vector(const Element::PropertiesType& rProperties); static void CalculateK0Stresses(Element& rElement, const ProcessInfo& rProcessInfo); static void CheckK0(const Properties& rProperties, IndexType ElementId); static void CheckK0MainDirection(const Properties& rProperties, IndexType ElementId); static void CheckOCRorPOP(const Properties& rProperties, IndexType ElementId); - static void CheckPhi(const Properties& rProperties, IndexType ElementId); static void CheckPoissonUnloadingReloading(const Properties& rProperties, IndexType ElementId); static void CheckSufficientMaterialParameters(const Properties& rProperties, IndexType ElementId); diff --git a/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.cpp b/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.cpp index 51cc0042fe4b..39e3972ee6d7 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.cpp +++ b/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.cpp @@ -83,6 +83,45 @@ double ConstitutiveLawUtilities::GetCohesion(const Properties& rProperties) } } +bool ConstitutiveLawUtilities::HasFrictionAngle(const Properties& rProperties) +{ + // Friction angle can be supplied either directly via GEO_FRICTION_ANGLE + // or via UMAT parameters + INDEX_OF_UMAT_PHI_PARAMETER. + return rProperties.Has(GEO_FRICTION_ANGLE) || + (rProperties.Has(INDEX_OF_UMAT_PHI_PARAMETER) && rProperties.Has(UMAT_PARAMETERS)); +} + +void ConstitutiveLawUtilities::ValidateFrictionAngle(const Properties& rProperties, IndexType ElementId) +{ + double phi = 0.0; + std::string phi_name; + + if (rProperties.Has(GEO_FRICTION_ANGLE)) { + phi = rProperties[GEO_FRICTION_ANGLE]; + phi_name = "GEO_FRICTION_ANGLE"; + } else if (rProperties.Has(INDEX_OF_UMAT_PHI_PARAMETER) && rProperties.Has(UMAT_PARAMETERS)) { + const auto phi_index = rProperties[INDEX_OF_UMAT_PHI_PARAMETER]; + const auto number_of_umat_parameters = static_cast(rProperties[UMAT_PARAMETERS].size()); + + KRATOS_ERROR_IF(phi_index < 1 || phi_index > number_of_umat_parameters) + << "Properties ( " << rProperties.Id() << ") of element ( " << ElementId + << "): INDEX_OF_UMAT_PHI_PARAMETER (" << phi_index + << ") is not in range [1, size of UMAT_PARAMETERS]." << std::endl; + + phi = rProperties[UMAT_PARAMETERS][phi_index - 1]; + phi_name = "Phi"; + } + + if (phi_name == "") { + KRATOS_ERROR << "Properties ( " << rProperties.Id() << ") of element ( " << ElementId + << ") does not have GEO_FRICTION_ANGLE nor INDEX_OF_UMAT_PHI_PARAMETER." << std::endl; + } + + KRATOS_ERROR_IF(phi < 0.0 || phi > 90.0) + << "Properties ( " << rProperties.Id() << ") of element ( " << ElementId << "): " << phi_name + << " (" << phi << " degrees) should be between 0 and 90 degrees." << std::endl; +} + double ConstitutiveLawUtilities::GetFrictionAngleInDegrees(const Properties& rProperties) { if (rProperties.Has(GEO_FRICTION_ANGLE)) { diff --git a/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.h b/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.h index 266f02a0c5d2..4f4f8cc2d213 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.h +++ b/applications/GeoMechanicsApplication/custom_utilities/constitutive_law_utilities.h @@ -35,6 +35,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) ConstitutiveLawUtilities double detF); static double GetCohesion(const Properties& rProperties); + static bool HasFrictionAngle(const Properties& rProperties); + static void ValidateFrictionAngle(const Properties& rProperties, IndexType ElementId); static double GetFrictionAngleInDegrees(const Properties& rProperties); static double GetFrictionAngleInRadians(const Properties& rProperties); diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_k0_procedure_process.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_k0_procedure_process.cpp index ba320e76d1fc..8c4451b44c6d 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_k0_procedure_process.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_processes/test_apply_k0_procedure_process.cpp @@ -488,6 +488,24 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandPOPandNu_UR_3 KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance); } +KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_Values_3D, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + // Arrange + auto p_properties = std::make_shared(); + p_properties->SetValue(K0_VALUE_XX, 0.5); + p_properties->SetValue(K0_VALUE_YY, 0.5); + p_properties->SetValue(K0_VALUE_ZZ, 1.0); + p_properties->SetValue(K0_MAIN_DIRECTION, 2); + const auto initial_stress_vector = UblasUtilities::CreateVector({0.0, -10.0, -10.0, 27.0, 10.0, 5.0}); + + // Act + const auto actual_stress_vector = ApplyK0ProcedureOnStubElement(p_properties, initial_stress_vector); + + // Assert + const auto expected_stress_vector = UblasUtilities::CreateVector({-5, -5, -10, 0, 0, 0}); + KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance); +} + KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfProcessHasCorrectMaterialData, KratosGeoMechanicsFastSuiteWithoutKernel) { // Arrange @@ -525,9 +543,9 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfProcessHasCorrectMaterialData, Krat p_element->GetProperties().SetValue(K0_MAIN_DIRECTION, 1); KRATOS_EXPECT_EXCEPTION_IS_THROWN( - process.Check(), - "Insufficient material data for K0 procedure process for element 1. No K0_NC, " - "(INDEX_OF_UMAT_PHI_PARAMETER and UMAT_PARAMETERS) or (K0_VALUE_XX, _YY and _ZZ found).") + process.Check(), " Insufficient material data for K0 procedure process for property Id of " + "0 for element 1. No K0_NC, (INDEX_OF_UMAT_PHI_PARAMETER and " + "UMAT_PARAMETERS) or (K0_VALUE_XX, _YY and _ZZ found).") p_element->GetProperties().SetValue(K0_VALUE_XX, -0.5); p_element->GetProperties().SetValue(K0_VALUE_YY, -0.5); p_element->GetProperties().SetValue(K0_VALUE_ZZ, -0.5); @@ -544,15 +562,15 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfProcessHasCorrectMaterialData, Krat p_element->GetProperties().SetValue(K0_VALUE_ZZ, 0.5); p_element->GetProperties().SetValue(POISSON_UNLOADING_RELOADING, 0.75); - KRATOS_EXPECT_EXCEPTION_IS_THROWN( - process.Check(), - "POISSON_UNLOADING_RELOADING (0.75) is not in range [-1.0, 0.5> for element 1.") + KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(), + "POISSON_UNLOADING_RELOADING (0.75) is not in range [-1.0, " + "0.5> for property Id of 0 for element 1.") p_element->GetProperties().SetValue(POISSON_UNLOADING_RELOADING, 0.25); KRATOS_EXPECT_EXCEPTION_IS_THROWN( - process.Check(), "Insufficient material data for K0 procedure process for " - "element 1. Poisson unloading-reloading, OCR and POP functionality cannot " - "be combined with K0_VALUE_XX, _YY and _ZZ.") + process.Check(), "Insufficient material data for K0 procedure process for property Id of 0 " + "for element 1. Poisson unloading-reloading, OCR and POP functionality " + "cannot be combined with K0_VALUE_XX, _YY and _ZZ.") p_element->GetProperties().Erase(K0_VALUE_XX); p_element->GetProperties().Erase(K0_VALUE_YY); @@ -563,23 +581,26 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfProcessHasCorrectMaterialData, Krat umat_parameters[0] = -30.0; p_element->GetProperties().SetValue(UMAT_PARAMETERS, umat_parameters); KRATOS_EXPECT_EXCEPTION_IS_THROWN( - process.Check(), - "INDEX_OF_UMAT_PHI_PARAMETER (2) is not in range 1, size of UMAT_PARAMETERS for element 1") + process.Check(), "Properties ( 0) of element ( 1): INDEX_OF_UMAT_PHI_PARAMETER (2) is not " + "in range [1, size of UMAT_PARAMETERS].") p_element->GetProperties().SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 1); - KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(), - "Phi (-30) should be between 0 and 90 degrees for element 1.") + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + process.Check(), + "Properties ( 0) of element ( 1): Phi (-30 degrees) should be between 0 and 90 degrees.") umat_parameters[0] = 30.0; p_element->GetProperties().SetValue(UMAT_PARAMETERS, umat_parameters); p_element->GetProperties().SetValue(OCR, 0.5); - KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(), - "OCR (0.5) should be in the range [1.0,-> for element 1.") + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + process.Check(), + "Error: OCR (0.5) should be in the range [1.0,-> for property Id of 0 for element 1.") p_element->GetProperties().Erase(OCR); p_element->GetProperties().SetValue(POP, -100.0); - KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(), - "POP (-100) should be in the range [0.0,-> for element 1.") + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + process.Check(), + "POP (-100) should be in the range [0.0,-> for property Id of 0 element 1.") p_element->GetProperties().Erase(POP); p_element->GetProperties().Erase(INDEX_OF_UMAT_PHI_PARAMETER); @@ -589,6 +610,7 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfProcessHasCorrectMaterialData, Krat "K0_NC (-0.5) should be in the range [0.0,-> for element 1.") p_element->GetProperties().SetValue(K0_NC, 0.5); + p_element->GetProperties().SetValue(GEO_FRICTION_ANGLE, 35.0); KRATOS_EXPECT_EQ(process.Check(), 0); } @@ -601,7 +623,7 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfModelPartHasElements, KratosGeoMech k0_settings.AddString("model_part_name", "dummy"); KRATOS_EXPECT_EXCEPTION_IS_THROWN((ApplyK0ProcedureProcess{model, k0_settings}.Check()), - "ApplyK0ProcedureProces has no elements in modelpart dummy") + "ApplyK0ProcedureProcess has no elements in modelpart dummy") } KRATOS_TEST_CASE_IN_SUITE(CheckInfoK0ProcedureProcess, KratosGeoMechanicsFastSuiteWithoutKernel) diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_utilities/test_constitutive_law_utilities.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_utilities/test_constitutive_law_utilities.cpp index e050c5ae7ea0..a6a5aa2bfb57 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_utilities/test_constitutive_law_utilities.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_utilities/test_constitutive_law_utilities.cpp @@ -183,4 +183,75 @@ KRATOS_TEST_CASE_IN_SUITE(ConstitutiveLawUtilities_CalculateK0NCFromFrictionAngl 1.0 - 0.5 * std::numbers::sqrt3, Defaults::absolute_tolerance); } +KRATOS_TEST_CASE_IN_SUITE(ConstitutiveLawUtilities_HasFrictionAngle, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + // 1) GEO_FRICTION_ANGLE provided + auto properties = Properties{}; + properties.SetValue(GEO_FRICTION_ANGLE, 30.0); + KRATOS_EXPECT_EQ(ConstitutiveLawUtilities::HasFrictionAngle(properties), true); + + // 2) UMAT_PARAMETERS + INDEX_OF_UMAT_PHI_PARAMETER provided + properties = Properties{}; + auto umat_parameters = UblasUtilities::CreateVector({2.0, 30.0}); + properties.SetValue(UMAT_PARAMETERS, umat_parameters); + properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 2); + KRATOS_EXPECT_EQ(ConstitutiveLawUtilities::HasFrictionAngle(properties), true); + + // 3) only INDEX_OF_UMAT_PHI_PARAMETER provided -> should be false + properties = Properties{}; + properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 1); + KRATOS_EXPECT_EQ(ConstitutiveLawUtilities::HasFrictionAngle(properties), false); + + // 4) only UMAT_PARAMETERS provided -> should be false + properties = Properties{}; + umat_parameters = UblasUtilities::CreateVector({2.0, 30.0}); + properties.SetValue(UMAT_PARAMETERS, umat_parameters); + KRATOS_EXPECT_EQ(ConstitutiveLawUtilities::HasFrictionAngle(properties), false); + + // 5) neither provided -> false + properties = Properties{}; + KRATOS_EXPECT_EQ(ConstitutiveLawUtilities::HasFrictionAngle(properties), false); +} + +KRATOS_TEST_CASE_IN_SUITE(ConstitutiveLawUtilities_ValidateFrictionAngle, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + constexpr std::size_t element_id = 1; + + // Valid: GEO_FRICTION_ANGLE provided + auto properties = Properties{}; + properties.SetValue(GEO_FRICTION_ANGLE, 30.0); + EXPECT_NO_THROW(ConstitutiveLawUtilities::ValidateFrictionAngle(properties, element_id)); + + // Valid: UMAT_PARAMETERS + INDEX_OF_UMAT_PHI_PARAMETER provided + properties = Properties{}; + auto umat_parameters = UblasUtilities::CreateVector({2.0, 30.0}); + properties.SetValue(UMAT_PARAMETERS, umat_parameters); + properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 2); + EXPECT_NO_THROW(ConstitutiveLawUtilities::ValidateFrictionAngle(properties, element_id)); + + // Missing both -> error + properties = Properties{}; + KRATOS_EXPECT_EXCEPTION_IS_THROWN(ConstitutiveLawUtilities::ValidateFrictionAngle(properties, element_id), "Properties ( 0) of element ( 1) does not have GEO_FRICTION_ANGLE nor INDEX_OF_UMAT_PHI_PARAMETER."); + + // GEO_FRICTION_ANGLE out of range -> error + properties = Properties{}; + properties.SetValue(GEO_FRICTION_ANGLE, -1.0); + KRATOS_EXPECT_EXCEPTION_IS_THROWN(ConstitutiveLawUtilities::ValidateFrictionAngle(properties, element_id), " Properties ( 0) of element ( 1): GEO_FRICTION_ANGLE (-1 degrees) should be between 0 and 90 degrees."); + + // UMAT phi value out of range -> error + properties = Properties{}; + umat_parameters = UblasUtilities::CreateVector({2.0, -5.0}); + properties.SetValue(UMAT_PARAMETERS, umat_parameters); + properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 2); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + ConstitutiveLawUtilities::ValidateFrictionAngle(properties, element_id), + " Properties ( 0) of element ( 1): Phi (-5 degrees) should be between 0 and 90 degrees."); + + // INDEX_OF_UMAT_PHI_PARAMETER out of bounds -> error + properties = Properties{}; + umat_parameters = UblasUtilities::CreateVector({2.0, 30.0}); + properties.SetValue(UMAT_PARAMETERS, umat_parameters); + properties.SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 3); + KRATOS_EXPECT_EXCEPTION_IS_THROWN(ConstitutiveLawUtilities::ValidateFrictionAngle(properties, element_id), "Properties ( 0) of element ( 1): INDEX_OF_UMAT_PHI_PARAMETER (3) is not in range [1, size of UMAT_PARAMETERS]."); +} } // namespace Kratos::Testing