diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab.py b/applications/GeoMechanicsApplication/tests/test_element_lab.py
index 555fe87e98df..572b26cadae0 100644
--- a/applications/GeoMechanicsApplication/tests/test_element_lab.py
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab.py
@@ -13,52 +13,52 @@ class KratosGeoMechanicsLabElementTests(KratosGeoUnittest.TestCase):
     """
     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.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, assert_all_integration_points=True)
+        file_path = test_helper.get_file_path(os.path.join('test_element_lab', 'test_triaxial', 'drained'))
+        expected_disp = test_helper.get_values_from_csv_as_vectors(
+            Path(file_path) / "expected_disp.csv",
+            ["node_id"],
+            ["disp_x", "disp_y", "disp_z"])
+        expected_stress = test_helper.get_values_from_csv_as_vectors(
+            Path(file_path) / "expected_stress.csv",
+            ["element_id", "ip_index"],
+            ["stress_xx", "stress_yy", "stress_zz", "stress_xy", "stress_yz", "stress_xz"])
+        expected_strain = test_helper.get_values_from_csv_as_vectors(
+            Path(file_path) / "expected_strain.csv",
+            ["element_id", "ip_index"],
+            ["strain_xx", "strain_yy", "strain_zz", "strain_xy", "strain_yz", "strain_xz"])
+        self._run_triaxial_regression_test(file_path, 'triaxial_test_output.post.res', expected_disp, expected_stress,
+                                           expected_strain, 4)
 
     def test_triaxial_undrained(self):
         """
         Regression test for the undrained triaxial experiment.
         """
-        expected_disp = [[0.0, -0.2, 0.0], [0.045, -0.2, 0.0], [0.0, -0.1, 0.0], [ 0.045, -0.1, 0.0],[0.0, 0.0, 0.0], [0.09, -0.2, 0.0], [0.09, -0.1, 0.0], [0.045, 0.0, 0.0], [0.09, 0.0, 0.0]]
-        expected_stress = [[-100.0, -4740.0, -100.0, 0.0, 0.0, 0.0], [-100.0, -4740.0, -100.0, 0.0, 0.0, 0.0]]
-        expected_strain = [[0.09, -0.2, 0.09, 0.0, 0.0, 0.0], [0.09, -0.2, 0.09, 0.0, 0.0, 0.0]]
-        self._run_triaxial_regression_test('undrained', 'triaxial_undrained_test_output.post.res', expected_disp,
+        file_path = test_helper.get_file_path(os.path.join('test_element_lab', 'test_triaxial', 'undrained'))
+        expected_disp = test_helper.get_values_from_csv_as_vectors(
+            Path(file_path) / "expected_disp.csv",
+            ["node_id"],
+            ["disp_x", "disp_y", "disp_z"])
+        expected_stress = test_helper.get_values_from_csv_as_vectors(
+            Path(file_path) / "expected_stress.csv",
+            ["element_id", "ip_index"],
+            ["stress_xx", "stress_yy", "stress_zz", "stress_xy", "stress_yz", "stress_xz"])
+        expected_strain = test_helper.get_values_from_csv_as_vectors(
+            Path(file_path) / "expected_strain.csv",
+            ["element_id", "ip_index"],
+            ["strain_xx", "strain_yy", "strain_zz", "strain_xy", "strain_yz", "strain_xz"])
+        self._run_triaxial_regression_test(file_path, '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_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))
+    def _run_triaxial_regression_test(self, file_path, output_file_name, expected_displacement, expected_stress, expected_strain, precision_places):
         test_helper.run_kratos(file_path)
-
-        # Read the output files from the simulation for comparison.
         reader = GiDOutputFileReader()
         result = reader.read_output_from(os.path.join(file_path, output_file_name))
         time = 1.0
-
-        # Assert the displacement in all nodes in all directions.
-        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)
-
-        if assert_all_integration_points:
-            self.assert_integration_point_tensors(result, "CAUCHY_STRESS_TENSOR",
-                self._make_integration_point_tensor_entries(expected_stress, num_elements=2, num_integration_points_per_element=3),
-                time ,precision_places)
-            self.assert_integration_point_tensors(result, "ENGINEERING_STRAIN_TENSOR",
-                self._make_integration_point_tensor_entries(expected_strain, num_elements=2, num_integration_points_per_element=3),
-                time, precision_places)
-        else:
-            self.assert_integration_point_tensors(result, "CAUCHY_STRESS_TENSOR", 
-                self._make_integration_point_tensor_entries(expected_stress, num_elements=2, num_integration_points_per_element=1),
-                time, precision_places)
-            self.assert_integration_point_tensors(result, "ENGINEERING_STRAIN_TENSOR", 
-                self._make_integration_point_tensor_entries(expected_strain, num_elements=2, num_integration_points_per_element=1),
-                time, precision_places)
+        node_ids = sorted(expected_displacement)
+        expected_disp = [expected_displacement[node_id] for node_id in node_ids]
+        self.assert_nodal_values_at_time(result, "DISPLACEMENT", node_ids, expected_disp, time, precision_places)
+        self.assert_integration_point_tensors(result, "CAUCHY_STRESS_TENSOR", expected_stress, time, precision_places)
+        self.assert_integration_point_tensors(result, "ENGINEERING_STRAIN_TENSOR", expected_strain, time, precision_places)
             
     def test_oedometer_drained(self):
         """Regression test for the oedometer experiment on a linear elastic model with constant pore water pressure."""
@@ -122,38 +122,34 @@ def test_crs_drained(self):
         """Regression test for the CRS experiment with constant pore water pressure."""
         stage_name = 'drained'
         nr_of_phases = 5
-        nr_of_elements = 2
-        nr_of_integration_points_per_element = 3
-        total_integration_points = nr_of_elements * nr_of_integration_points_per_element
-
-        expected_strains = [
-            self._repeat_tensor([0.0, -0.1, 0.0, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([0.0, -0.05, 0.0, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([0.0, -0.25, 0.0, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([0.0, -0.25, 0.0, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([0.0, -0.40, 0.0, 0.0, 0.0, 0.0], total_integration_points),
-        ]
         expected_water_pressures = [[0.0] * 9 for _ in range(nr_of_phases)]
 
-        expected_stresses = [
-            self._repeat_tensor([-4e+05, -12e+05, -4e+05, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-2e+05, -6e+05, -2e+05, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-1e+06, -3e+06, -1e+06, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-1e+06, -3e+06, -1e+06, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-1.6e+06, -4.8e+06, -1.6e+06, 0.0, 0.0, 0.0], total_integration_points),
-        ]
+        linear_elastic_file_path = test_helper.get_file_path(os.path.join('test_element_lab', 'test_crs', stage_name, 'linear_elastic'))
+        linear_elastic_stresses = test_helper.get_values_from_csv_grouped(
+            Path(linear_elastic_file_path) / "expected_stress.csv",
+            ['phase'], ['element_id', 'ip_index'],
+            ["stress_xx", "stress_yy", "stress_zz", "stress_yz", "stress_xz", "stress_xy"])
+        linear_elastic_strains = test_helper.get_values_from_csv_grouped(
+            Path(linear_elastic_file_path) / "expected_strain.csv",
+            ['phase'], ['element_id', 'ip_index'],
+            ["strain_xx", "strain_yy", "strain_zz", "strain_yz", "strain_xz", "strain_xy"])
+        expected_stresses = [linear_elastic_stresses[phase] for phase in range(1, nr_of_phases + 1)]
+        expected_strains = [linear_elastic_strains[phase] for phase in range(1, nr_of_phases + 1)]
 
         self._run_crs_test_for_model(stage_name, "linear_elastic")
         self._check_crs_results(stage_name, "linear_elastic", nr_of_phases, expected_stresses, expected_strains, expected_water_pressures, delta_stress=1.0)
 
-        # mohr coulomb expectations per phase
-        expected_stresses = [
-            self._repeat_tensor([-9.01878e+05, -2.22843e+06, -9.01878e+05, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-5.01878e+05, -1.02843e+06, -5.01878e+05, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-2.25697e+06, -5.56725e+06, -2.25697e+06, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-2.25697e+06, -5.56725e+06, -2.25697e+06, 0.0, 0.0, 0.0], total_integration_points),
-            self._repeat_tensor([-3.61205e+06, -8.90607e+06, -3.61205e+06, 0.0, 0.0, 0.0], total_integration_points)
-            ]
+        mohr_coulomb_file_path = test_helper.get_file_path(os.path.join('test_element_lab', 'test_crs', stage_name, 'mohr_coulomb'))
+        mohr_coulomb_stresses = test_helper.get_values_from_csv_grouped(
+            Path(mohr_coulomb_file_path) / "expected_stress.csv",
+            ['phase'], ['element_id', 'ip_index'],
+            ["stress_xx", "stress_yy", "stress_zz", "stress_xy", "stress_yz", "stress_xz"])
+        mohr_coulomb_strains = test_helper.get_values_from_csv_grouped(
+            Path(mohr_coulomb_file_path) / "expected_strain.csv",
+            ['phase'], ['element_id', 'ip_index'],
+            ["strain_xx", "strain_yy", "strain_zz", "strain_xy", "strain_yz", "strain_xz"])
+        expected_stresses = [mohr_coulomb_stresses[phase] for phase in range(1, nr_of_phases + 1)]
+        expected_strains = [mohr_coulomb_strains[phase] for phase in range(1, nr_of_phases + 1)]
 
         self._run_crs_test_for_model(stage_name, "mohr_coulomb")
         self._check_crs_results(stage_name, "mohr_coulomb", nr_of_phases, expected_stresses, expected_strains, expected_water_pressures, delta_stress=1.0)
@@ -171,13 +167,9 @@ def _check_crs_results(self, stage_name, model_name, nr_of_phases, expected_stre
         node_ids = list(range(1, 10))
         for i in range(nr_of_phases):
             self.assert_integration_point_tensors(
-                result, "CAUCHY_STRESS_TENSOR",
-                self._make_integration_point_tensor_entries(expected_stresses[i], num_elements=2, num_integration_points_per_element=3),
-                times[i], places=places_stress, delta=delta_stress)
+                result, "CAUCHY_STRESS_TENSOR", expected_stresses[i], times[i], places=places_stress, delta=delta_stress)
             self.assert_integration_point_tensors(
-                result, "ENGINEERING_STRAIN_TENSOR",
-                self._make_integration_point_tensor_entries(expected_strains[i], num_elements=2, num_integration_points_per_element=3),
-                times[i], places=places_strain)
+                result, "ENGINEERING_STRAIN_TENSOR", expected_strains[i], times[i], places=places_strain)
             self.assert_nodal_values_at_time(result, "WATER_PRESSURE", node_ids, expected_water_pressures[i], times[i], places=places_water_pressure)
 
     def test_triaxial_comp_6n(self):
@@ -249,9 +241,6 @@ 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 _repeat_tensor(self, tensor, count):
-        return [list(tensor) for _ in range(count)]
-
     def _make_integration_point_tensor_entries(self, expected_tensors, num_elements, num_integration_points_per_element):
         result = {}
         idx = 0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/linear_elastic/expected_strain.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/linear_elastic/expected_strain.csv
new file mode 100644
index 000000000000..7021c9fd342c
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/linear_elastic/expected_strain.csv
@@ -0,0 +1,31 @@
+phase,element_id,ip_index,strain_xx,strain_yy,strain_zz,strain_xy,strain_yz,strain_xz
+1,1,0,0.0,-0.1,0.0,0.0,0.0,0.0
+1,1,1,0.0,-0.1,0.0,0.0,0.0,0.0
+1,1,2,0.0,-0.1,0.0,0.0,0.0,0.0
+1,2,0,0.0,-0.1,0.0,0.0,0.0,0.0
+1,2,1,0.0,-0.1,0.0,0.0,0.0,0.0
+1,2,2,0.0,-0.1,0.0,0.0,0.0,0.0
+2,1,0,0.0,-0.05,0.0,0.0,0.0,0.0
+2,1,1,0.0,-0.05,0.0,0.0,0.0,0.0
+2,1,2,0.0,-0.05,0.0,0.0,0.0,0.0
+2,2,0,0.0,-0.05,0.0,0.0,0.0,0.0
+2,2,1,0.0,-0.05,0.0,0.0,0.0,0.0
+2,2,2,0.0,-0.05,0.0,0.0,0.0,0.0
+3,1,0,0.0,-0.25,0.0,0.0,0.0,0.0
+3,1,1,0.0,-0.25,0.0,0.0,0.0,0.0
+3,1,2,0.0,-0.25,0.0,0.0,0.0,0.0
+3,2,0,0.0,-0.25,0.0,0.0,0.0,0.0
+3,2,1,0.0,-0.25,0.0,0.0,0.0,0.0
+3,2,2,0.0,-0.25,0.0,0.0,0.0,0.0
+4,1,0,0.0,-0.25,0.0,0.0,0.0,0.0
+4,1,1,0.0,-0.25,0.0,0.0,0.0,0.0
+4,1,2,0.0,-0.25,0.0,0.0,0.0,0.0
+4,2,0,0.0,-0.25,0.0,0.0,0.0,0.0
+4,2,1,0.0,-0.25,0.0,0.0,0.0,0.0
+4,2,2,0.0,-0.25,0.0,0.0,0.0,0.0
+5,1,0,0.0,-0.4,0.0,0.0,0.0,0.0
+5,1,1,0.0,-0.4,0.0,0.0,0.0,0.0
+5,1,2,0.0,-0.4,0.0,0.0,0.0,0.0
+5,2,0,0.0,-0.4,0.0,0.0,0.0,0.0
+5,2,1,0.0,-0.4,0.0,0.0,0.0,0.0
+5,2,2,0.0,-0.4,0.0,0.0,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/linear_elastic/expected_stress.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/linear_elastic/expected_stress.csv
new file mode 100644
index 000000000000..ab416f355072
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/linear_elastic/expected_stress.csv
@@ -0,0 +1,31 @@
+phase,element_id,ip_index,stress_xx,stress_yy,stress_zz,stress_xy,stress_yz,stress_xz
+1,1,0,-400000,-1200000,-400000,0.0,0.0,0.0
+1,1,1,-400000,-1200000,-400000,0.0,0.0,0.0
+1,1,2,-400000,-1200000,-400000,0.0,0.0,0.0
+1,2,0,-400000,-1200000,-400000,0.0,0.0,0.0
+1,2,1,-400000,-1200000,-400000,0.0,0.0,0.0
+1,2,2,-400000,-1200000,-400000,0.0,0.0,0.0
+2,1,0,-200000,-600000,-200000,0.0,0.0,0.0
+2,1,1,-200000,-600000,-200000,0.0,0.0,0.0
+2,1,2,-200000,-600000,-200000,0.0,0.0,0.0
+2,2,0,-200000,-600000,-200000,0.0,0.0,0.0
+2,2,1,-200000,-600000,-200000,0.0,0.0,0.0
+2,2,2,-200000,-600000,-200000,0.0,0.0,0.0
+3,1,0,-1000000,-3000000,-1000000,0.0,0.0,0.0
+3,1,1,-1000000,-3000000,-1000000,0.0,0.0,0.0
+3,1,2,-1000000,-3000000,-1000000,0.0,0.0,0.0
+3,2,0,-1000000,-3000000,-1000000,0.0,0.0,0.0
+3,2,1,-1000000,-3000000,-1000000,0.0,0.0,0.0
+3,2,2,-1000000,-3000000,-1000000,0.0,0.0,0.0
+4,1,0,-1000000,-3000000,-1000000,0.0,0.0,0.0
+4,1,1,-1000000,-3000000,-1000000,0.0,0.0,0.0
+4,1,2,-1000000,-3000000,-1000000,0.0,0.0,0.0
+4,2,0,-1000000,-3000000,-1000000,0.0,0.0,0.0
+4,2,1,-1000000,-3000000,-1000000,0.0,0.0,0.0
+4,2,2,-1000000,-3000000,-1000000,0.0,0.0,0.0
+5,1,0,-1600000,-4800000,-1600000,0.0,0.0,0.0
+5,1,1,-1600000,-4800000,-1600000,0.0,0.0,0.0
+5,1,2,-1600000,-4800000,-1600000,0.0,0.0,0.0
+5,2,0,-1600000,-4800000,-1600000,0.0,0.0,0.0
+5,2,1,-1600000,-4800000,-1600000,0.0,0.0,0.0
+5,2,2,-1600000,-4800000,-1600000,0.0,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/mohr_coulomb/expected_strain.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/mohr_coulomb/expected_strain.csv
new file mode 100644
index 000000000000..7021c9fd342c
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/mohr_coulomb/expected_strain.csv
@@ -0,0 +1,31 @@
+phase,element_id,ip_index,strain_xx,strain_yy,strain_zz,strain_xy,strain_yz,strain_xz
+1,1,0,0.0,-0.1,0.0,0.0,0.0,0.0
+1,1,1,0.0,-0.1,0.0,0.0,0.0,0.0
+1,1,2,0.0,-0.1,0.0,0.0,0.0,0.0
+1,2,0,0.0,-0.1,0.0,0.0,0.0,0.0
+1,2,1,0.0,-0.1,0.0,0.0,0.0,0.0
+1,2,2,0.0,-0.1,0.0,0.0,0.0,0.0
+2,1,0,0.0,-0.05,0.0,0.0,0.0,0.0
+2,1,1,0.0,-0.05,0.0,0.0,0.0,0.0
+2,1,2,0.0,-0.05,0.0,0.0,0.0,0.0
+2,2,0,0.0,-0.05,0.0,0.0,0.0,0.0
+2,2,1,0.0,-0.05,0.0,0.0,0.0,0.0
+2,2,2,0.0,-0.05,0.0,0.0,0.0,0.0
+3,1,0,0.0,-0.25,0.0,0.0,0.0,0.0
+3,1,1,0.0,-0.25,0.0,0.0,0.0,0.0
+3,1,2,0.0,-0.25,0.0,0.0,0.0,0.0
+3,2,0,0.0,-0.25,0.0,0.0,0.0,0.0
+3,2,1,0.0,-0.25,0.0,0.0,0.0,0.0
+3,2,2,0.0,-0.25,0.0,0.0,0.0,0.0
+4,1,0,0.0,-0.25,0.0,0.0,0.0,0.0
+4,1,1,0.0,-0.25,0.0,0.0,0.0,0.0
+4,1,2,0.0,-0.25,0.0,0.0,0.0,0.0
+4,2,0,0.0,-0.25,0.0,0.0,0.0,0.0
+4,2,1,0.0,-0.25,0.0,0.0,0.0,0.0
+4,2,2,0.0,-0.25,0.0,0.0,0.0,0.0
+5,1,0,0.0,-0.4,0.0,0.0,0.0,0.0
+5,1,1,0.0,-0.4,0.0,0.0,0.0,0.0
+5,1,2,0.0,-0.4,0.0,0.0,0.0,0.0
+5,2,0,0.0,-0.4,0.0,0.0,0.0,0.0
+5,2,1,0.0,-0.4,0.0,0.0,0.0,0.0
+5,2,2,0.0,-0.4,0.0,0.0,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/mohr_coulomb/expected_stress.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/mohr_coulomb/expected_stress.csv
new file mode 100644
index 000000000000..e36b2b8bb38e
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_crs/drained/mohr_coulomb/expected_stress.csv
@@ -0,0 +1,31 @@
+phase,element_id,ip_index,stress_xx,stress_yy,stress_zz,stress_xy,stress_yz,stress_xz
+1,1,0,-901878,-2228430,-901878,0.0,0.0,0.0
+1,1,1,-901878,-2228430,-901878,0.0,0.0,0.0
+1,1,2,-901878,-2228430,-901878,0.0,0.0,0.0
+1,2,0,-901878,-2228430,-901878,0.0,0.0,0.0
+1,2,1,-901878,-2228430,-901878,0.0,0.0,0.0
+1,2,2,-901878,-2228430,-901878,0.0,0.0,0.0
+2,1,0,-501878,-1028430,-501878,0.0,0.0,0.0
+2,1,1,-501878,-1028430,-501878,0.0,0.0,0.0
+2,1,2,-501878,-1028430,-501878,0.0,0.0,0.0
+2,2,0,-501878,-1028430,-501878,0.0,0.0,0.0
+2,2,1,-501878,-1028430,-501878,0.0,0.0,0.0
+2,2,2,-501878,-1028430,-501878,0.0,0.0,0.0
+3,1,0,-2256970,-5567250,-2256970,0.0,0.0,0.0
+3,1,1,-2256970,-5567250,-2256970,0.0,0.0,0.0
+3,1,2,-2256970,-5567250,-2256970,0.0,0.0,0.0
+3,2,0,-2256970,-5567250,-2256970,0.0,0.0,0.0
+3,2,1,-2256970,-5567250,-2256970,0.0,0.0,0.0
+3,2,2,-2256970,-5567250,-2256970,0.0,0.0,0.0
+4,1,0,-2256970,-5567250,-2256970,0.0,0.0,0.0
+4,1,1,-2256970,-5567250,-2256970,0.0,0.0,0.0
+4,1,2,-2256970,-5567250,-2256970,0.0,0.0,0.0
+4,2,0,-2256970,-5567250,-2256970,0.0,0.0,0.0
+4,2,1,-2256970,-5567250,-2256970,0.0,0.0,0.0
+4,2,2,-2256970,-5567250,-2256970,0.0,0.0,0.0
+5,1,0,-3612050,-8906070,-3612050,0.0,0.0,0.0
+5,1,1,-3612050,-8906070,-3612050,0.0,0.0,0.0
+5,1,2,-3612050,-8906070,-3612050,0.0,0.0,0.0
+5,2,0,-3612050,-8906070,-3612050,0.0,0.0,0.0
+5,2,1,-3612050,-8906070,-3612050,0.0,0.0,0.0
+5,2,2,-3612050,-8906070,-3612050,0.0,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_disp.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_disp.csv
new file mode 100644
index 000000000000..45c283d7e2e5
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_disp.csv
@@ -0,0 +1,10 @@
+node_id,disp_x,disp_y,disp_z
+1,0.0,-0.2,0.0
+2,0.0527776,-0.2,0.0
+3,0.0,-0.100033,0.0
+4,0.0524025,-0.0996931,0.0
+5,0.0,0.0,0.0
+6,0.105197,-0.2,0.0
+7,0.105114,-0.100049,0.0
+8,0.0524406,0.0,0.0
+9,0.104632,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_strain.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_strain.csv
new file mode 100644
index 000000000000..efdfdbe142e6
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_strain.csv
@@ -0,0 +1,7 @@
+element_id,ip_index,strain_xx,strain_yy,strain_zz,strain_xy,strain_yz,strain_xz
+1,0,0.104863,-0.19973,0.104946,0.000440186,0.0,0.0
+1,1,0.104939,-0.199861,0.105022,0.000100528,0.0,0.0
+1,2,0.10444,-0.199182,0.104773,6.24562e-05,0.0,0.0
+2,0,0.1055,-0.200303,0.104922,3.84218e-05,0.0,0.0
+2,1,0.104915,-0.200106,0.105298,-4.97694e-06,0.0,0.0
+2,2,0.105631,-0.200818,0.105335,0.000287088,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_stress.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_stress.csv
new file mode 100644
index 000000000000..5f68683fef68
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/drained/expected_stress.csv
@@ -0,0 +1,7 @@
+element_id,ip_index,stress_xx,stress_yy,stress_zz,stress_xy,stress_yz,stress_xz
+1,0,-99.9808,-252.622,-99.9806,0.193199,0.0,0.0
+1,1,-99.8828,-252.381,-99.8827,0.0448209,0.0,0.0
+1,2,-100.045,-252.78,-100.045,0.0314238,0.0,0.0
+2,0,-99.9991,-252.668,-99.9991,0.00846584,0.0,0.0
+2,1,-99.967,-252.589,-99.967,-0.0103576,0.0,0.0
+2,2,-100.064,-252.827,-100.064,0.137162,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_disp.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_disp.csv
new file mode 100644
index 000000000000..7fcad9839bff
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_disp.csv
@@ -0,0 +1,10 @@
+node_id,disp_x,disp_y,disp_z
+1,0.0,-0.2,0.0
+2,0.045,-0.2,0.0
+3,0.0,-0.1,0.0
+4,0.045,-0.1,0.0
+5,0.0,0.0,0.0
+6,0.09,-0.2,0.0
+7,0.09,-0.1,0.0
+8,0.045,0.0,0.0
+9,0.09,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_strain.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_strain.csv
new file mode 100644
index 000000000000..f6bee494190e
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_strain.csv
@@ -0,0 +1,3 @@
+element_id,ip_index,strain_xx,strain_yy,strain_zz,strain_xy,strain_yz,strain_xz
+1,0,0.09,-0.2,0.09,0.0,0.0,0.0
+2,0,0.09,-0.2,0.09,0.0,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_stress.csv b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_stress.csv
new file mode 100644
index 000000000000..d59890480aa0
--- /dev/null
+++ b/applications/GeoMechanicsApplication/tests/test_element_lab/test_triaxial/undrained/expected_stress.csv
@@ -0,0 +1,3 @@
+element_id,ip_index,stress_xx,stress_yy,stress_zz,stress_xy,stress_yz,stress_xz
+1,0,-100.0,-4740.0,-100.0,0.0,0.0,0.0
+2,0,-100.0,-4740.0,-100.0,0.0,0.0,0.0
diff --git a/applications/GeoMechanicsApplication/tests/test_helper.py b/applications/GeoMechanicsApplication/tests/test_helper.py
index 4cb33ff11a61..fbe7481c671e 100644
--- a/applications/GeoMechanicsApplication/tests/test_helper.py
+++ b/applications/GeoMechanicsApplication/tests/test_helper.py
@@ -1,6 +1,7 @@
 from typing import Dict, Any
 import sys,os
 import math
+import csv
 
 sys.path.append(os.path.join('..','..','..'))
 
@@ -46,6 +47,34 @@ def run_kratos(file_path, model=None):
     os.chdir(cwd)
     return simulation
 
+def get_values_from_csv(csv_filepath, keys, values):
+    result = {}
+    with open(csv_filepath, newline = "") as csv_file:
+        reader = csv.DictReader(csv_file)
+        for row in reader:           
+            if len(keys) == 1:
+                key = int(row[keys[0]])
+            else:
+                key = tuple(int(row[key]) for key in keys)
+            expected_values = {value: float(row[value]) for value in values}
+            result[key] = expected_values
+    return result
+
+def get_values_from_csv_as_vectors(csv_filepath, keys, values):
+    return {key: [row[value] for value in values]
+            for key, row in get_values_from_csv(csv_filepath, keys, values).items()}
+
+def get_values_from_csv_grouped(csv_filepath, grouping_keys, entry_keys, values):
+    grouped_results = {}
+    combined = get_values_from_csv(csv_filepath, grouping_keys + entry_keys, values)
+    for key, values_dict in combined.items():
+        if len(grouping_keys) == 1:
+            group_key = key[0] if isinstance(key, tuple) else key
+        else:
+            group_key = tuple(key[:len(grouping_keys)])
+        entry_key = key[len(grouping_keys):] if isinstance(key, tuple) else ()
+        grouped_results.setdefault(group_key, {})[entry_key] = [values_dict[v] for v in values]
+    return grouped_results
 
 def get_displacement(simulation):
     """