Skip to content

Commit

Permalink
[GeoMechanicsApplication] Standard K0 procedure for 3D (#12850)
Browse files Browse the repository at this point in the history
* clean setup and added 3D tests

* used K0_MAIN_DIRECTION 2

* check for K0_MAIN_DIRECTION > 0

* used the patch

* format

* added integration test for 3D

* added dimension into CheckK0MainDirection

* format

* used mock

* updated integration test

* format

* removed commented line
  • Loading branch information
markelov208 authored Nov 19, 2024
1 parent 5845ac2 commit 6015a90
Show file tree
Hide file tree
Showing 8 changed files with 831 additions and 33 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,17 @@ void ApplyK0ProcedureProcess::CheckK0MainDirection(const Properties& rProperties
{
KRATOS_ERROR_IF(!rProperties.Has(K0_MAIN_DIRECTION))
<< "K0_MAIN_DIRECTION is not defined for element " << ElementId << "." << std::endl;
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;

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;
}
}

void ApplyK0ProcedureProcess::CheckK0(const Properties& rProperties, IndexType ElementId)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -13,21 +13,25 @@
#include "custom_processes/apply_k0_procedure_process.h"
#include "geo_mechanics_application_variables.h"
#include "geo_mechanics_fast_suite.h"
#include "geometries/quadrilateral_2d_4.h"
#include "includes/element.h"
#include "processes/structured_mesh_generator_process.h"
#include "stub_linear_elastic_law.h"
#include "test_utilities.h"
#include <custom_constitutive/incremental_linear_elastic_law.h>

#include <boost/algorithm/cxx11/all_of.hpp>
#include <boost/algorithm/cxx11/none_of.hpp>
#include <boost/numeric/ublas/assignment.hpp>
#include <gmock/gmock.h>

namespace
{

using namespace Kratos;

class MockConstitutiveLaw : public GeoIncrementalLinearElasticLaw
{
public:
MOCK_METHOD(std::size_t, WorkingSpaceDimension, (), (override));
};

class StubElement : public Element
{
public:
Expand Down Expand Up @@ -89,26 +93,14 @@ namespace

ModelPart& PrepareTestModelPart(Model& rModel)
{
auto& result = rModel.CreateModelPart("dummy");

// Set up the test model part mesh
auto p_point_1 = Kratos::make_intrusive<Node>(1, 0.0, 0.0, 0.0);
auto p_point_2 = Kratos::make_intrusive<Node>(2, 0.0, 1.0, 0.0);
auto p_point_3 = Kratos::make_intrusive<Node>(3, 1.0, 1.0, 0.0);
auto p_point_4 = Kratos::make_intrusive<Node>(4, 1.0, 0.0, 0.0);
Quadrilateral2D4<Node> domain_geometry(p_point_1, p_point_2, p_point_3, p_point_4);

Parameters mesher_parameters(R"({
"number_of_divisions": 2,
"element_name": "Element2D3N",
"condition_name": "LineCondition",
"create_skin_sub_model_part": true
})");
StructuredMeshGeneratorProcess(domain_geometry, result, mesher_parameters).Execute();

auto& result = rModel.CreateModelPart("dummy");
auto p_dummy_law = std::make_shared<Testing::StubLinearElasticLaw>();
auto& r_model_part_properties = result.GetProperties(0);
r_model_part_properties.SetValue(CONSTITUTIVE_LAW, p_dummy_law);
auto p_model_part_properties = result.pGetProperties(0);
p_model_part_properties->SetValue(CONSTITUTIVE_LAW, p_dummy_law);

auto p_element = make_intrusive<StubElement>();
p_element->SetProperties(p_model_part_properties);
result.AddElement(p_element);

return result;
}
Expand All @@ -126,7 +118,7 @@ namespace Kratos::Testing
{

KRATOS_TEST_CASE_IN_SUITE(AllElementsConsiderDiagonalEntriesOnlyAndNoShearWhenUseStandardProcedureFlagIsNotDefined,
KratosGeoMechanicsFastSuite)
KratosGeoMechanicsFastSuiteWithoutKernel)
{
Model model;
auto& r_model_part = PrepareTestModelPart(model);
Expand All @@ -136,11 +128,11 @@ KRATOS_TEST_CASE_IN_SUITE(AllElementsConsiderDiagonalEntriesOnlyAndNoShearWhenUs
ApplyK0ProcedureProcess process{r_model_part, k0_settings};
process.ExecuteInitialize();

KRATOS_EXPECT_TRUE(boost::algorithm::all_of(r_model_part.Elements(), ElementConsidersDiagonalEntriesOnlyAndNoShear))
KRATOS_EXPECT_TRUE(ElementConsidersDiagonalEntriesOnlyAndNoShear(r_model_part.Elements()[0]))
}

KRATOS_TEST_CASE_IN_SUITE(AllElementsConsiderDiagonalEntriesOnlyAndNoShearWhenUsingStandardProcedure,
KratosGeoMechanicsFastSuite)
KratosGeoMechanicsFastSuiteWithoutKernel)
{
Model model;
auto& r_model_part = PrepareTestModelPart(model);
Expand All @@ -150,11 +142,11 @@ KRATOS_TEST_CASE_IN_SUITE(AllElementsConsiderDiagonalEntriesOnlyAndNoShearWhenUs
ApplyK0ProcedureProcess process{r_model_part, k0_settings};
process.ExecuteInitialize();

KRATOS_EXPECT_TRUE(boost::algorithm::all_of(r_model_part.Elements(), ElementConsidersDiagonalEntriesOnlyAndNoShear))
KRATOS_EXPECT_TRUE(ElementConsidersDiagonalEntriesOnlyAndNoShear(r_model_part.Elements()[0]))
}

KRATOS_TEST_CASE_IN_SUITE(NoneOfElementsConsiderDiagonalEntriesOnlyAndNoShearWhenNotUsingStandardProcedure,
KratosGeoMechanicsFastSuite)
KratosGeoMechanicsFastSuiteWithoutKernel)
{
Model model;
auto& r_model_part = PrepareTestModelPart(model);
Expand All @@ -164,10 +156,11 @@ KRATOS_TEST_CASE_IN_SUITE(NoneOfElementsConsiderDiagonalEntriesOnlyAndNoShearWhe
ApplyK0ProcedureProcess process{r_model_part, k0_settings};
process.ExecuteInitialize();

KRATOS_EXPECT_TRUE(boost::algorithm::none_of(r_model_part.Elements(), ElementConsidersDiagonalEntriesOnlyAndNoShear))
KRATOS_EXPECT_FALSE(ElementConsidersDiagonalEntriesOnlyAndNoShear(r_model_part.Elements()[0]))
}

KRATOS_TEST_CASE_IN_SUITE(UseStandardProcedureFlagIsInEffectDuringProcessExecutionOnly, KratosGeoMechanicsFastSuite)
KRATOS_TEST_CASE_IN_SUITE(UseStandardProcedureFlagIsInEffectDuringProcessExecutionOnly,
KratosGeoMechanicsFastSuiteWithoutKernel)
{
Model model;
auto& r_model_part = PrepareTestModelPart(model);
Expand All @@ -178,7 +171,7 @@ KRATOS_TEST_CASE_IN_SUITE(UseStandardProcedureFlagIsInEffectDuringProcessExecuti
process.ExecuteInitialize(); // start considering diagonal entries only and no shear
process.ExecuteFinalize(); // stop considering diagonal entries only and no shear

KRATOS_EXPECT_TRUE(boost::algorithm::none_of(r_model_part.Elements(), ElementConsidersDiagonalEntriesOnlyAndNoShear))
KRATOS_EXPECT_FALSE(ElementConsidersDiagonalEntriesOnlyAndNoShear(r_model_part.Elements()[0]))
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NC, KratosGeoMechanicsFastSuiteWithoutKernel)
Expand All @@ -199,6 +192,24 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NC, KratosGeoMecha
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NC_3D, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
auto p_properties = std::make_shared<Properties>();
p_properties->SetValue(K0_NC, 0.5);
p_properties->SetValue(K0_MAIN_DIRECTION, 2);
Vector initial_stress_vector{6};
initial_stress_vector <<= 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
Vector expected_stress_vector{6};
expected_stress_vector <<= -5.0, -5.0, -10.0, 0.0, 0.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithPhi, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
Expand All @@ -222,6 +233,29 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithPhi, KratosGeoMechani
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithPhi_3D, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
auto p_properties = std::make_shared<Properties>();
p_properties->SetValue(INDEX_OF_UMAT_PHI_PARAMETER, 1);
p_properties->SetValue(NUMBER_OF_UMAT_PARAMETERS, 1);
Vector umat_parameters{1};
umat_parameters[0] = 30.0;
p_properties->SetValue(UMAT_PARAMETERS, umat_parameters);
p_properties->SetValue(K0_MAIN_DIRECTION, 2);

Vector initial_stress_vector{6};
initial_stress_vector <<= 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
Vector expected_stress_vector{6};
expected_stress_vector <<= -5.0, -5.0, -10.0, 0.0, 0.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandOCR, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
Expand All @@ -241,6 +275,25 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandOCR, KratosGe
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandOCR_3D, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
auto p_properties = std::make_shared<Properties>();
p_properties->SetValue(K0_NC, 0.5);
p_properties->SetValue(K0_MAIN_DIRECTION, 2);
p_properties->SetValue(OCR, 1.5);
Vector initial_stress_vector{6};
initial_stress_vector <<= 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
Vector expected_stress_vector{6};
expected_stress_vector <<= -7.5, -7.5, -10.0, 0.0, 0.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandOCRandNu_UR, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
Expand All @@ -261,6 +314,26 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandOCRandNu_UR,
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandOCRandNu_UR_3D, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
auto p_properties = std::make_shared<Properties>();
p_properties->SetValue(K0_NC, 0.5);
p_properties->SetValue(K0_MAIN_DIRECTION, 2);
p_properties->SetValue(OCR, 1.5);
p_properties->SetValue(POISSON_UNLOADING_RELOADING, 0.25);
Vector initial_stress_vector{6};
initial_stress_vector <<= 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
Vector expected_stress_vector{6};
expected_stress_vector <<= -7 * 10.0 / 12.0, -7 * 10.0 / 12.0, -10.0, 0.0, 0.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandPOP, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
Expand All @@ -280,6 +353,25 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandPOP, KratosGe
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandPOP_3D, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
auto p_properties = std::make_shared<Properties>();
p_properties->SetValue(K0_NC, 0.5);
p_properties->SetValue(K0_MAIN_DIRECTION, 2);
p_properties->SetValue(POP, 50.0);
Vector initial_stress_vector{6};
initial_stress_vector <<= 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
Vector expected_stress_vector{6};
expected_stress_vector <<= -30.0, -30.0, -10.0, 0.0, 0.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandPOPandNu_UR, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
Expand All @@ -300,6 +392,26 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandPOPandNu_UR,
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureIsAppliedCorrectlyWithK0_NCandPOPandNu_UR_3D, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
auto p_properties = std::make_shared<Properties>();
p_properties->SetValue(K0_NC, 0.5);
p_properties->SetValue(K0_MAIN_DIRECTION, 2);
p_properties->SetValue(POP, 50.0);
p_properties->SetValue(POISSON_UNLOADING_RELOADING, 0.25);
Vector initial_stress_vector{6};
initial_stress_vector <<= 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
Vector expected_stress_vector{6};
expected_stress_vector <<= -30.0 + 50.0 / 3.0, -30.0 + 50.0 / 3.0, -10.0, 0.0, 0.0, 0.0;
KRATOS_EXPECT_VECTOR_NEAR(actual_stress_vector, expected_stress_vector, Defaults::absolute_tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfProcessHasCorrectMaterialData, KratosGeoMechanicsFastSuiteWithoutKernel)
{
// Arrange
Expand All @@ -313,14 +425,27 @@ KRATOS_TEST_CASE_IN_SUITE(K0ProcedureChecksIfProcessHasCorrectMaterialData, Krat
const auto k0_settings = Parameters{};
ApplyK0ProcedureProcess process{r_model_part, k0_settings};

auto mock_constitutive_law = std::make_shared<MockConstitutiveLaw>();
p_element->GetProperties().SetValue(CONSTITUTIVE_LAW, mock_constitutive_law);
EXPECT_CALL(*mock_constitutive_law, WorkingSpaceDimension()).WillRepeatedly(testing::Return(2));

// Act & Assert
KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(),
"K0_MAIN_DIRECTION is not defined for element 1.");

p_element->GetProperties().SetValue(K0_MAIN_DIRECTION, 4);

EXPECT_CALL(*mock_constitutive_law, WorkingSpaceDimension()).WillRepeatedly(testing::Return(1));
KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(), "dimension should be 2 or 3 for element 1.")

EXPECT_CALL(*mock_constitutive_law, WorkingSpaceDimension()).WillRepeatedly(testing::Return(2));
KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(),
"K0_MAIN_DIRECTION should be 0 or 1 for element 1.")

EXPECT_CALL(*mock_constitutive_law, WorkingSpaceDimension()).WillRepeatedly(testing::Return(3));
KRATOS_EXPECT_EXCEPTION_IS_THROWN(process.Check(),
"K0_MAIN_DIRECTION should be 0, 1 or 2 for element 1.")

p_element->GetProperties().SetValue(K0_MAIN_DIRECTION, 1);
KRATOS_EXPECT_EXCEPTION_IS_THROWN(
process.Check(),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,17 @@ def assert_stresses_at_integration_points(self, cauchy_stress_tensors, integrati
for integration_point in integration_points:
self.assert_stresses_at_integration_point(cauchy_stress_tensors, integration_point, expected_horizontal_stress, expected_vertical_stress, rel_tol)

def compare_stresses_at_integration_point(self, cauchy_stress_tensors, integration_point, k0_nc, rel_tol):
"""
Verifies whether the computed stresses are (nearly) equal to some expected values at the
given integration point. Note that this function assumes there are no shear stresses!
"""
element_id, integration_point_index = integration_point
stress_tensor = cauchy_stress_tensors[element_id-1][integration_point_index]
z_stress = stress_tensor[2, 2]
self.assertIsClose(stress_tensor[0, 0], z_stress*k0_nc, rel_tol=rel_tol, msg=f"X stress at integration point {integration_point_index} of element {element_id}")
self.assertIsClose(stress_tensor[1, 1], z_stress*k0_nc, rel_tol=rel_tol, msg=f"Y stress at integration point {integration_point_index} of element {element_id}")

def test_k0_procedure_k0_nc(self):
"""
Test to check if CAUCHY_STRESS_XX is correctly derived from CAUCHY_STRESS_YY using K0_NC
Expand Down Expand Up @@ -472,5 +483,26 @@ def test_k0_procedure_for_tilted_layers(self):
integration_point = (2, 0) # far right
self.assert_stresses_at_integration_point(cauchy_stress_tensors, integration_point, expected_vertical_stress=-22084, expected_horizontal_stress=-11042, rel_tol=0.02)

def test_k0_procedure_for_3d(self):
"""
Test to check whether the effective stress distribution is in line with regression data.
To this end, we test the horizontal, vertical and shear stresses at a selection
of integration points (defined as pairs of element IDs and integration point indices).
The settings are taken from lysmer_column3d_hexa_in_Z test.
"""
test_path = test_helper.get_file_path(os.path.join("test_k0_procedure_process", "test_k0_procedure_k0_nc_3D"))
simulation = test_helper.run_kratos(test_path)

cauchy_stress_tensors = test_helper.get_on_integration_points(simulation, Kratos.CAUCHY_STRESS_TENSOR)

k0_nc = 0.5
# Check the stresses at a few integration points
integration_point = (1, 1) # bottom
self.compare_stresses_at_integration_point(cauchy_stress_tensors, integration_point, k0_nc, rel_tol=0.02)
integration_point = (10, 1) # middle
self.compare_stresses_at_integration_point(cauchy_stress_tensors, integration_point, k0_nc, rel_tol=0.02)
integration_point = (20, 1) # top
self.compare_stresses_at_integration_point(cauchy_stress_tensors, integration_point, k0_nc, rel_tol=0.02)

if __name__ == '__main__':
KratosUnittest.main()
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
This folder contains the test cases:

- [K<sub>0</sub> procedure normal consolidation](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc/README.md)
- [K<sub>0</sub> procedure normal consolidation in 3D](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_3D/README.md)
- [K<sub>0</sub> procedure normal consolidation with layers](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_layers/README.md)
- [K<sub>0</sub> procedure normal consolidation with OCR](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_ocr/README.md)
- [K<sub>0</sub> procedure normal consolidation with OCR_field](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_k0_procedure_process/test_k0_procedure_k0_nc_ocr_field/README.md)
Expand Down
Loading

0 comments on commit 6015a90

Please sign in to comment.