[GeoMechanicsApplication] add a new oedometer test to the test_element_lab (#14403)

* make a mesh, material and project parameters file for the oedometer test (unit convention = kN)

* insert assertion function and apply it in the triaxial test

* add oedometer test with assertion function

* apply assert function on displacement assertion triaxial test

* formatting of json files for the oedometer test

* remove triaxial material parameters from the json format exception list

* minor changes to triaxial test readme

* add readme for the oedometer test incl schematic

* changes in mesh, project and material files according to comments and changes in undrained branch

* add stage name drained and add functions in triaxial test

* remove old code

* change to geo drainage type

* change nr of checked ip's in triaxial test

Co-authored-by: Copilot <copilot@github.com>

* make code use general helper functions

Co-authored-by: Copilot <copilot@github.com>

* change output file type

* change the poisson ratio

* change output file name

* remove value for the bulk modulus fluid

* apply stepwise loading + remove unused material parameters and set to si units (N instead of kN)

* si units in assertion and stepwise loading

* fix fixed base constraint

* use axisymmetric element type and loading

* Update readme

Co-authored-by: Copilot <copilot@github.com>

* apply review comments

---------

Co-authored-by: Copilot <copilot@github.com>

Author: Sabine-van-Dijk

applications/GeoMechanicsApplication/exceptions_for_json_formatting.txt

@@ -67,7 +67,6 @@ tests/test_constant_strip_load_2d/MaterialParameters_stage_1.json
 tests/test_constant_strip_load_2d/ProjectParameters.json
 tests/test_darcy_law.gid/MaterialParameters.json
 tests/test_darcy_law.gid/ProjectParameters.json
-tests/test_element_lab/test_triaxial/MaterialParameters.json
 tests/test_inclinded_phreatic_line_smaller_line.gid/MaterialParameters.json
 tests/test_inclinded_phreatic_line_smaller_line.gid/ProjectParameters.json
 tests/test_inclinded_phreatic_surface.gid/MaterialParameters.json

applications/GeoMechanicsApplication/tests/test_element_lab.py

@@ -1,3 +1,4 @@
+from csv import reader
 import os
 
 import KratosMultiphysics.KratosUnittest as KratosUnittest
@@ -7,18 +8,15 @@
 
 class KratosGeoMechanicsLabElementTests(KratosUnittest.TestCase):
     """
-    This class contains some element tests, such as triaxial and oedometer tests
+    This class contains some element tests, such as triaxial and oedometer tests.
     """
-    def test_triaxial(self):
-        """
-        Regression test for the triaxial experiment.
-        """
+    def test_triaxial_drained(self):
+        """Regression test for the triaxial experiment on a mohr coulomb model with a constant pore water pressure."""
         expected_disp = [[0.0, -0.2, 0.0], [0.0527776, -0.2, 0.0], [0.0, -0.100033, 0.0], [0.0524025, -0.0996931, 0.0], [0.0, 0.0, 0.0], [0.105197, -0.2, 0.0], [0.105114, -0.100049, 0.0], [0.0524406, 0.0, 0.0], [0.104632, 0.0, 0.0]]
-        expected_stress = [[-99.9808, -252.622, -99.9806, 0.193199, 0.0, 0.0], [-99.9991, -252.668, -99.9991, 0.00846584, 0, 0]]
-        expected_strain = [[0.104863, -0.19973, 0.104946, 0.000440186, 0.0, 0.0], [0.1055, -0.200303, 0.104922, 3.84218e-05, 0, 0]]
+        expected_stress = [[-99.9808, -252.622, -99.9806, 0.193199, 0.0, 0.0], [-99.8828, -252.381, -99.8827, 0.0448209, 0, 0], [-100.045, -252.78, -100.045, 0.0314238, 0, 0], [-99.9991, -252.668, -99.9991, 0.00846584, 0, 0], [-99.967, -252.589, -99.967, -0.0103576, 0, 0], [-100.064, -252.827, -100.064, 0.137162, 0, 0]]
+        expected_strain = [[0.104863, -0.19973, 0.104946, 0.000440186, 0.0, 0.0], [0.104939, -0.199861, 0.105022, 0.000100528, 0, 0], [0.10444, -0.199182, 0.104773, 6.24562e-05, 0, 0], [0.1055, -0.200303, 0.104922, 3.84218e-05, 0, 0], [0.104915, -0.200106, 0.105298, -4.97694e-06, 0, 0], [0.105631, -0.200818, 0.105335, 0.000287088, 0, 0]]
         self._run_triaxial_regression_test('drained', 'triaxial_test_output.post.res', expected_disp, expected_stress,
-                                           expected_strain, 4)
-
+                                           expected_strain, 4, assert_all_integration_points=True)
 
     def test_triaxial_undrained(self):
         """
@@ -30,10 +28,10 @@ def test_triaxial_undrained(self):
         self._run_triaxial_regression_test('undrained', 'triaxial_undrained_test_output.post.res', expected_disp,
                                            expected_stress, expected_strain, 4)
 
-    def _run_triaxial_regression_test(self, stage_name, output_file_name, expected_disp, expected_stress,
-                                      expected_strain, precision_places):
-        test_name = os.path.join('test_triaxial', stage_name)
-        file_path = test_helper.get_file_path(os.path.join('test_element_lab', test_name))
+    def _run_triaxial_regression_test(self, stage_name, output_file_name, expected_displacement, expected_stress,
+                                      expected_strain, precision_places, assert_all_integration_points=False):
+        test_name = 'test_triaxial'
+        file_path = test_helper.get_file_path(os.path.join('test_element_lab', test_name, stage_name))
         test_helper.run_kratos(file_path)
 
         # Read the output files from the simulation for comparison.
@@ -42,59 +40,86 @@ def _run_triaxial_regression_test(self, stage_name, output_file_name, expected_d
         time = 1.0
 
         # Assert the displacement in all nodes in all directions.
-        for node_id, expected_node_displacement in enumerate(expected_disp, start=1):
+        for node_id, expected_node_displacement in enumerate(expected_displacement, start=1):
             node_displacement = reader.nodal_values_at_time("DISPLACEMENT", time, result, [node_id])[0]
             self.assertVectorAlmostEqual(node_displacement, expected_node_displacement, precision_places)
 
-        self._assert_first_ip_tensor(reader, result, "CAUCHY_STRESS_TENSOR", expected_stress, precision_places, time)
-        self._assert_first_ip_tensor(reader, result, "ENGINEERING_STRAIN_TENSOR", expected_strain, precision_places, time)
+        if assert_all_integration_points:
+            self._assert_integration_point_tensors(reader, result, "CAUCHY_STRESS_TENSOR", expected_stress,
+                                                  precision_places, time)
+            self._assert_integration_point_tensors(reader, result, "ENGINEERING_STRAIN_TENSOR", expected_strain,
+                                                  precision_places, time)
+        else:
+            self._assert_first_ip_tensor(reader, result, "CAUCHY_STRESS_TENSOR", expected_stress,
+                                         precision_places, time)
+            self._assert_first_ip_tensor(reader, result, "ENGINEERING_STRAIN_TENSOR", expected_strain,
+                                         precision_places, time)
+            
+    def test_oedometer_drained(self):
+        """Regression test for the oedometer experiment on a linear elastic model with constant pore water pressure."""
+        expected_stress_per_ip = [-1e+06/3.0, -1e+06, -1e+06/3.0, 0.0, 0.0, 0.0]
+        expected_stress = [expected_stress_per_ip] * 6
+        expected_y_displacements = [-0.0208, -0.0417, -0.0625, -0.0833]
+        displacement_times = [0.25, 0.5, 0.75, 1.0]
+        top_nodes = [7, 8, 9]
+        self._run_oedometer_regression_test('drained', 'test_oedometer_output.post.res', expected_stress,
+                                            expected_y_displacements, displacement_times, top_nodes, 0, 4)
+        
+    def _run_oedometer_regression_test(self, stage_name, output_file_name, expected_stress,
+                                       expected_y_displacements, displacement_times, top_nodes, precision_places_stress, precision_places_displacement=4):
+        """Run the oedometer regression test and validate stress tensors and y-displacements in time."""
+        test_name = 'test_oedometer'
+        file_path = test_helper.get_file_path(os.path.join('test_element_lab', test_name, stage_name))
+        test_helper.run_kratos(file_path)
 
-    def _assert_first_ip_tensor(self, reader, result, variable_name, expected_values, precision_places, time=1.0):
-        for element_id, expected_tensor in enumerate(expected_values, start=1):
-            tensor = reader.element_integration_point_values_at_time(variable_name, time, result, [element_id], [0])[0]
-            self._assert_integration_point_tensor_results(tensor, expected_tensor, precision_places, variable_name)
+        reader = GiDOutputFileReader()
+        result = reader.read_output_from(os.path.join(file_path, output_file_name))
+
+        self._assert_integration_point_tensors(reader, result, "CAUCHY_STRESS_TENSOR", expected_stress,
+                                              precision_places_stress, time=1.0)
+
+        for time, expected_y in zip(displacement_times, expected_y_displacements):
+            self._assert_y_displacements_at_time(reader, result, top_nodes, expected_y, precision_places_displacement, time)
 
     def test_triaxial_comp_6n(self):
         """
         Drained compression triaxial test on Mohr-Coulomb model with axisymmetric 2D6N elements
-        It consistes of two calculation phases:
+        It consists of two calculation phases:
         1) apply confining stress of -100 kPa
         2) apply deviatoric stress of -200 kPa
         """
-        project_path = test_helper.get_file_path(os.path.join('test_element_lab', 'triaxial_comp_6n'))
-
+        test_name = 'triaxial_comp_6n'
+        project_path = test_helper.get_file_path(os.path.join('test_element_lab', test_name))
         n_stages = 2
         run_multiple_stages.run_stages(project_path, n_stages)
-
         reader = GiDOutputFileReader()
-        self._assert_triaxial_comp_6n_stage(reader, project_path, "triaxial_comp_6n_stage1.post.res", 1.0,
-                                            [-100.0, -100.0, -100.0], 3, 3)
-        self._assert_triaxial_comp_6n_stage(reader, project_path, "triaxial_comp_6n_stage2.post.res", 1.25,
-                                            [-100.0, -300.0, -100.0], 3, 2)
 
+        result = reader.read_output_from(os.path.join(project_path, "triaxial_comp_6n_stage1.post.res"))
+        expected_stress = [[-100.0, -100.0, -100.0, 0.0, 0.0, 0.0]] * 6
+        self._assert_integration_point_tensors(reader, result, "CAUCHY_STRESS_TENSOR", expected_stress, 3, time=1.0)
+
+        result = reader.read_output_from(os.path.join(project_path, "triaxial_comp_6n_stage2.post.res"))
+        expected_stress = [[-100.0, -300.0, -100.0, 0.0, 0.0, 0.0]] * 6
+        self._assert_integration_point_tensors(reader, result, "CAUCHY_STRESS_TENSOR", expected_stress, 2, time=1.25)
 
     def test_oedometer_ULFEM(self):
         """
         Oedometer test on a linear elastic model with 2D6N elements
         """
         test_name = 'oedometer_ULFEM'
-        project_path = test_helper.get_file_path(os.path.join('test_element_lab', test_name))
-        simulation = test_helper.run_kratos(project_path)
+        file_path = test_helper.get_file_path(os.path.join('test_element_lab', test_name))
+        simulation = test_helper.run_kratos(file_path)
         effective_stresses = test_helper.get_cauchy_stress_tensor(simulation)
-        self._assert_oedometer_effective_stresses(effective_stresses, -1000000.0, 2)
+        self._assert_oedometer_effective_stresses(effective_stresses, -1e+06, 2)
 
         top_node_nbrs = [1, 2]
-        output_file_path = os.path.join(project_path, test_name+'.post.res')
-        output_reader = GiDOutputFileReader()
-        output_data = output_reader.read_output_from(output_file_path)
-        displacements = GiDOutputFileReader.nodal_values_at_time("DISPLACEMENT", 0.1, output_data, node_ids=top_node_nbrs)
-        self._assert_y_displacements(displacements, -0.00990099, 6)
-
-        displacements = GiDOutputFileReader.nodal_values_at_time("DISPLACEMENT", 0.7, output_data, node_ids=top_node_nbrs)
-        self._assert_y_displacements(displacements, -0.0654206, 6)
+        output_file = os.path.join(file_path, test_name+'.post.res')
+        reader = GiDOutputFileReader()
+        result = reader.read_output_from(output_file)
 
-        displacements = GiDOutputFileReader.nodal_values_at_time("DISPLACEMENT", 1.0, output_data, node_ids=top_node_nbrs)
-        self._assert_y_displacements(displacements, -0.0909090909516868, 6)
+        self._assert_y_displacements_at_time(reader, result, top_node_nbrs, -0.0099, 4, 0.1)
+        self._assert_y_displacements_at_time(reader, result, top_node_nbrs, -0.0654, 4, 0.7)
+        self._assert_y_displacements_at_time(reader, result, top_node_nbrs, -0.0909, 4, 1.0)
 
     def test_oedometer_ULFEM_diff_order(self):
         """
@@ -104,23 +129,13 @@ def test_oedometer_ULFEM_diff_order(self):
         project_path = test_helper.get_file_path(os.path.join('test_element_lab', test_name))
         simulation = test_helper.run_kratos(project_path)
         effective_stresses = test_helper.get_cauchy_stress_tensor(simulation)
-        displacements = test_helper.get_displacement(simulation)
-        self._assert_oedometer_effective_stresses(effective_stresses, -1000.0, 3)
+        self._assert_oedometer_effective_stresses(effective_stresses, -1e+03, 3)
 
-        y_displacements = [displacement[1] for displacement in displacements]
+        output_file = os.path.join(project_path, test_name + '.post.res')
+        reader = GiDOutputFileReader()
+        result = reader.read_output_from(output_file)
         top_node_nbrs = [1]
-        for top_node_nbr in top_node_nbrs:
-            self.assertAlmostEqual(-1e-04, y_displacements[top_node_nbr], 6)
-
-    def _assert_triaxial_comp_6n_stage(self, reader, project_path, output_file_name, time, expected_stress_vector,
-                                       number_of_integration_points_per_element, places):
-        output_data = reader.read_output_from(os.path.join(project_path, output_file_name))
-        stress_vectors_per_element = reader.element_integration_point_values_at_time("CAUCHY_STRESS_TENSOR", time, output_data)
-        self.assertEqual(2, len(stress_vectors_per_element))
-
-        for element_stress_vectors in stress_vectors_per_element:
-            self.assertEqual(number_of_integration_points_per_element, len(element_stress_vectors))
-            self._assert_integration_point_tensor_results(element_stress_vectors, expected_stress_vector, places, "CAUCHY_STRESS_TENSOR")
+        self._assert_y_displacements_at_time(reader, result, top_node_nbrs, -1e-4, 6, time=1.0)
 
     def _assert_oedometer_effective_stresses(self, effective_stresses, expected_yy, places_yy):
         for element in effective_stresses:
@@ -129,14 +144,32 @@ def _assert_oedometer_effective_stresses(self, effective_stresses, expected_yy,
                 self.assertAlmostEqual(expected_yy, integration_point[1,1], places_yy)
                 self.assertAlmostEqual(0.0, integration_point[2,2], 3)
 
-    def _assert_y_displacements(self, displacements, expected_y, places):
+    def _assert_y_displacements_at_time(self, reader, result, node_ids, expected_y, places, time=1.0):
+        displacements = reader.nodal_values_at_time("DISPLACEMENT", time, result, node_ids=node_ids)
         for displacement in displacements:
             self.assertAlmostEqual(expected_y, displacement[1], places)
 
     def _assert_integration_point_tensor_results(self, integration_point_tensors, expected_integration_point_tensor, places, result_name):
         for idx, ip_tensor in enumerate(integration_point_tensors):
-            self.assertVectorAlmostEqual(expected_integration_point_tensor[:2], ip_tensor[:2], places, msg = f"{result_name} component xx, yy, zz at integration point {idx}")
+            self.assertVectorAlmostEqual(expected_integration_point_tensor, ip_tensor, places, msg = f"{result_name} components at integration point {idx}")
+
+    def _assert_integration_point_tensors(self, reader, result, variable_name, expected_tensors, places, time):
+        """Assert tensor values for all integration points across both elements."""
+        for flat_index, expected_tensor in enumerate(expected_tensors):
+            element_id, ip_index = self._get_element_id_and_ip_index(flat_index)
+            tensor = reader.element_integration_point_values_at_time(variable_name, time, result, [element_id], [ip_index])[0]
+            self._assert_integration_point_tensor_results(tensor, expected_tensor, places, variable_name)
 
+    def _get_element_id_and_ip_index(self, flat_index):
+        """Map a flat integration-point index to the corresponding element and local integration point index."""
+        if flat_index < 3:
+            return 1, flat_index
+        return 2, flat_index - 3
+
+    def _assert_first_ip_tensor(self, reader, result, variable_name, expected_values, precision_places, time=1.0):
+        for element_id, expected_tensor in enumerate(expected_values, start=1):
+            tensor = reader.element_integration_point_values_at_time(variable_name, time, result, [element_id], [0])[0]
+            self._assert_integration_point_tensor_results(tensor, expected_tensor, precision_places, variable_name)
 
 if __name__ == '__main__':
     KratosUnittest.main()

applications/GeoMechanicsApplication/tests/test_element_lab/oedometer_ULFEM_diff_order/MaterialParameters.json

@@ -14,7 +14,6 @@
               "DENSITY_WATER"            :  1.0e3,
               "POROSITY"                 :  0.3,
               "BULK_MODULUS_SOLID"       :  1.0e12,
-              "BULK_MODULUS_FLUID"       :  2.0e-30,
               "PERMEABILITY_XX"          :  4.5e-30,
               "PERMEABILITY_YY"          :  4.5e-30,
               "PERMEABILITY_XY"          :  0.0,

applications/GeoMechanicsApplication/tests/test_element_lab/oedometer_ULFEM_diff_order/ProjectParameters.json

@@ -73,13 +73,13 @@
             "process_name": "GiDOutputProcess",
             "Parameters":    {
                 "model_part_name": "PorousDomain.porous_computational_model_part",
-                "output_name": "Oedometer_ULFEM_diff_order",
+                "output_name": "oedometer_ULFEM_diff_order",
                 "postprocess_parameters": {
                     "result_file_configuration": {
                         "gidpost_flags":       {
                             "WriteDeformedMeshFlag": "WriteUndeformed",
                             "WriteConditionsFlag":   "WriteElementsOnly",
-                            "GiDPostMode":           "GiD_PostBinary",
+                            "GiDPostMode":           "GiD_PostAscii",
                             "MultiFileFlag":         "SingleFile"
                         },
                         "file_label":          "step",

applications/GeoMechanicsApplication/tests/test_element_lab/test_oedometer/README.md

@@ -0,0 +1,30 @@
+# Oedometer Test
+
+This test simulates an oedometer test with a prescribed stepwise top loading process on a linear elastic soil model. It replicates a laboratory test used to determine soil properties, such as consolidation characteristics.
+In the laboratory, this is performed on a cylindrical volume of soil, where an increasing pressure is applied on the top of the cylinder. In the model test, the soil sample is emulated by two axisymmetric differential order elements.
+
+A schematic overview of the model is displayed in the figure below. Here the nodes are displayed as black dots. The two elements are numbered in red and are separated by the red line. The symmetry axis is displayed as a red dashed line.
+
+![MeshStructure](drawing.svg)
+
+## Setup
+
+The test is performed under the following conditions:
+
+- **Constraints**:
+    - The bottom nodes (1, 2, 3) are fixed in the X and Y directions.
+    - The side nodes and the symmetry axis (3, 6, 9 and 1, 4, 7) are fixed in the X direction.
+- **Material**:
+    - The material is described by the linear elastic model with the following parameters:
+        - Poisson's ratio = 0.25
+        - Young's modulus = 1e+07 Pa
+- **Loading Conditions**:
+    - A top load is applied in 4 equal steps:
+        - t = 0.01 to 0.25: load = -2.5e+05 Pa
+        - t = 0.26 to 0.50: load = -5e+05 Pa
+        - t = 0.51 to 0.75: load = -7.5e+05 Pa
+        - t = 0.76 to 1.0: load = -1e+06 Pa
+
+## Assertions
+For this test, the stresses at t = 1.0 are verified.
+Additionally, the displacement in the Y-direction of the top nodes is verified at four time steps: t = 0.25, 0.5, 0.75, and 1.0.

applications/GeoMechanicsApplication/tests/test_element_lab/test_oedometer/common/mesh.mdpa

@@ -0,0 +1,89 @@
+Begin Table 1 TIME LINE_NORMAL_LOAD_Y
+  0.0 0.0
+  0.01 2.5e+05
+  0.25 2.5e+05
+  0.26 5e+05
+  0.50 5e+05
+  0.51 7.5e+05
+  0.75 7.5e+05
+  0.76 1e+06
+  1.0 1e+06
+End Table
+
+Begin Properties 0
+End Properties
+
+Begin Properties 1
+End Properties
+
+Begin Nodes
+  1  0.0 0.0 0.0
+  2  0.5 0.0 0.0
+  3  1.0 0.0 0.0
+  4  0.0 0.5 0.0
+  5  0.5 0.5 0.0
+  6  1.0 0.5 0.0
+  7  0.0 1.0 0.0
+  8  0.5 1.0 0.0
+  9  1.0 1.0 0.0
+End Nodes
+
+Begin Elements SmallStrainUPwDiffOrderAxisymmetricElement2D6N
+  1   1 7 1 3 4 2 5
+  2   1 3 9 7 6 8 5
+End Elements
+
+Begin Conditions AxisymmetricLineNormalLoadDiffOrderCondition2D3N
+  1 0 9 7 8 
+End Conditions
+
+Begin SubModelPart Soil
+  Begin SubModelPartNodes
+    1
+    2
+    3
+    4
+    5
+    6
+    7
+    8
+    9
+  End SubModelPartNodes
+  Begin SubModelPartElements
+    1
+    2
+  End SubModelPartElements
+End SubModelPart
+
+Begin SubModelPart Top_load
+  Begin SubModelPartTables
+    1
+  End SubModelPartTables
+  Begin SubModelPartNodes
+    7
+    8
+    9
+  End SubModelPartNodes
+  Begin SubModelPartConditions
+    1
+  End SubModelPartConditions
+End SubModelPart
+
+Begin SubModelPart Fixed_base
+  Begin SubModelPartNodes
+    1
+    2
+    3
+  End SubModelPartNodes
+End SubModelPart
+
+Begin SubModelPart Sides_rollers
+  Begin SubModelPartNodes
+    3
+    6
+    9
+    1
+    4
+    7
+  End SubModelPartNodes
+End SubModelPart