diff --git a/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py b/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py index 162a749848c2..2737ea9d9609 100644 --- a/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py +++ b/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py @@ -36,6 +36,7 @@ from test_normal_load_on_1d_element import KratosGeoMechanicsNormalLoad1DTests from test_k0_procedure_process import KratosGeoMechanicsK0ProcedureProcessTests from test_geomechanics_solver import KratosGeoMechanicsSolverTests +from test_column_changing_waterlevel import KratosGeoMechanicsChangingWaterLevelTests def AssembleTestSuites(): ''' Populates the test suites to run. @@ -70,7 +71,8 @@ def AssembleTestSuites(): KratosGeoMechanicsParameterFieldTests, KratosGeoMechanicsNormalLoad1DTests, KratosGeoMechanicsK0ProcedureProcessTests, - KratosGeoMechanicsSolverTests + KratosGeoMechanicsSolverTests, + KratosGeoMechanicsChangingWaterLevelTests ] # Create an array with the selected tests diff --git a/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel.py b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel.py new file mode 100644 index 000000000000..4c456b0db0c0 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel.py @@ -0,0 +1,207 @@ +import os +import KratosMultiphysics.KratosUnittest as KratosUnittest +import test_helper + + +class KratosGeoMechanicsChangingWaterLevelTests(KratosUnittest.TestCase): + """ + This class contains tests which check the result of a changing water level process. + The test contains a column of 10m x 50m with a water level set to -25m at the first step. + """ + + def parameters(self): + """ + Set parameters for the test. + rho_wet: density of the soil in Mg/m3 + porosity: porosity of the soil + h_soil: height of the soil column in m + rho_water: density of the water in Mg/m3 + g: gravity in m/s2 + h_water_list: list of water levels in m + + Args: + - None + + Returns: + - None + """ + + self.rho_wet = 2.65 + self.porosity = 0.3 + self.h_soil = 0.0 + self.rho_water = 1.0 + self.g = -9.81 + self.h_water_list = [-25, -30, -35, -40, -35, -30, -25, -20, -15, -10] + + def calculate_water_pressure(self, rho_water, g, y, h_water): + """ + Calculate water pressure. + + Args: + - rho_water (float): density of water + - g (float): gravity + - y (float): y-coordinate of the gauss point + - h_water (float): water level + + Returns: + - float: water pressure + """ + + return rho_water * g * (h_water - y) if y < h_water else 0.0 + + def calculate_analytical_total_stress_yy(self, rho_water, g, y, h_water): + """ + Analytical solution for calculating the total stress in the column. + + Args: + - rho_water (float): density of water + - g (float): gravity + - y (float): y-coordinate of the gauss point + - h_water (float): water level + + Returns: + - list: analytical total stress at the y coordinate of gauss point + """ + + return ((1.0 - self.porosity) * self.rho_wet * g * (self.h_soil - y) + + self.porosity * self.calculate_water_pressure(rho_water, g, y, h_water)) + + def calculate_analytical_effective_stress_yy(self, rho_water, g, y, h_water): + """ + Analytical solution for calculating the effective stress in the column. + + Args: + - rho_water (float): density of water + - g (float): gravity + - y (float): y-coordinate of the gauss point + - h_water (float): water level + + Returns: + - list: analytical effective stress at the y coordinate of gauss point + """ + + return ((1.0 - self.porosity) * self.rho_wet * g * (self.h_soil - y) + + (self.porosity - 1.0) * self.calculate_water_pressure(rho_water, g, y, h_water)) + + def assert_almost_equal_values(self, values1, values2): + """ + Compare analytical vs numerical solution. + + Args: + - expected_values (list): analytical solution + - calculated_values_yy (list): numerical solution + + Returns: + - None + """ + + for value1, value2 in zip(values1, values2): + self.assertAlmostEqual(value1, value2, places=3) + + def get_y_coordinates_of_gauss_points(self, simulation): + """ + Get y-coordinates of the gauss points. + + Args: + - simulation (object): simulation object + + Returns: + - list: y-coordinates of the gauss points + """ + + y_coordinates = [] + for all_gauss_points_in_element in test_helper.get_gauss_coordinates(simulation): + for position in all_gauss_points_in_element: + y_coordinates.append(position[1]) + return y_coordinates + + def compare_effective_stresses(self, simulation_output, simulation): + """ + Compare analytical vs numerical solution for the effective stresses. + + Args: + - simulation (object): simulation object + - file_path (str): path to the file + + Returns: + - None + """ + + for h_water, effective_stresses in zip(self.h_water_list, simulation_output["results"]["CAUCHY_STRESS_TENSOR"]): + analytical_effective_stresses_yy =\ + [self.calculate_analytical_effective_stress_yy(self.rho_water, self.g, y, h_water) + for y in self.get_y_coordinates_of_gauss_points(simulation)] + + numerical_effective_stresses_yy = [] + for value in effective_stresses["values"]: + numerical_effective_stresses_yy.extend([stress[1] for stress in value["value"]]) + self.assert_almost_equal_values(analytical_effective_stresses_yy, numerical_effective_stresses_yy) + + def compare_total_stresses(self, simulation_output, simulation): + """ + Compare analytical vs numerical solution for the total stresses. + + Args: + - simulation_output (object): simulation output + - simulation (object): simulation object + + Returns: + - None + """ + + for h_water, total_stresses in zip(self.h_water_list, simulation_output["results"]["TOTAL_STRESS_TENSOR"]): + analytical_total_stresses_yy =\ + [self.calculate_analytical_total_stress_yy(self.rho_water, self.g, y, h_water) + for y in self.get_y_coordinates_of_gauss_points(simulation)] + + numerical_total_stresses_yy = [] + for value in total_stresses["values"]: + numerical_total_stresses_yy.extend([stress[1] for stress in value["value"]]) + self.assert_almost_equal_values(analytical_total_stresses_yy, numerical_total_stresses_yy) + + def compare_water_pressures(self, simulation_output, simulation): + """ + Compare analytical vs numerical solution for the water pressures. + + Args: + - simulation (object): simulation object + + Returns: + - None + """ + + y_coordinates = [node.Y0 for node in simulation._list_of_output_processes[0].model_part.Nodes] + for h_water, water_pressures in zip(self.h_water_list, simulation_output["results"]["WATER_PRESSURE"]): + analytical_water_pressures = [self.calculate_water_pressure(self.rho_water, self.g, y, h_water) + for y in y_coordinates] + numerical_water_pressures = [item["value"] for item in water_pressures["values"]] + self.assert_almost_equal_values(analytical_water_pressures, numerical_water_pressures) + + def test_side_pressure_boundaries(self): + """ + Test to check if CAUCHY_STRESS_YY, TOTAL_STRESS_YY and WATER_PRESSURE are consistent + with the analytical solution. + + Args: + - None + + Returns: + - None + """ + + test_name = os.path.join("test_column_changing_waterlevel", "test_side_pressure_boundaries") + file_path = test_helper.get_file_path(test_name) + # run simulation + simulation = test_helper.run_kratos(file_path) + # read results + result_file_path = os.path.join(file_path, "test_phreatic.post.res") + simulation_output = test_helper.read_numerical_results_from_post_res(result_file_path) + # compare results + self.parameters() + self.compare_effective_stresses(simulation_output, simulation) + self.compare_total_stresses(simulation_output, simulation) + self.compare_water_pressures(simulation_output, simulation) + + +if __name__ == '__main__': + KratosUnittest.main() diff --git a/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/MaterialParameters.json b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/MaterialParameters.json new file mode 100644 index 000000000000..d1279895729d --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/MaterialParameters.json @@ -0,0 +1,35 @@ +{ + "properties": [{ + "model_part_name": "PorousDomain.Soil_column", + "properties_id": 1, + "Material": { + "constitutive_law": { + "name" : "GeoLinearElasticPlaneStrain2DLaw" + }, + "Variables": { + "IGNORE_UNDRAINED" : false, + "YOUNG_MODULUS" : 10000, + "POISSON_RATIO" : 0.2, + "DENSITY_SOLID" : 2.65, + "DENSITY_WATER" : 1.0, + "POROSITY" : 0.3, + "BULK_MODULUS_SOLID" : 1.0e9, + "BULK_MODULUS_FLUID" : 2.0e6, + "PERMEABILITY_XX" : 4.5e+13, + "PERMEABILITY_YY" : 4.5e+13, + "PERMEABILITY_XY" : 0.0, + "DYNAMIC_VISCOSITY" : 8.90e-17, + "THICKNESS" : 1.0, + "BIOT_COEFFICIENT" : 1.0, + "RETENTION_LAW" : "SaturatedBelowPhreaticLevelLaw", + "SATURATED_SATURATION" : 1.0, + "RESIDUAL_SATURATION" : 1e-10, + "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE" : 2.561, + "VAN_GENUCHTEN_GN" : 1.377, + "VAN_GENUCHTEN_GL" : 1.25, + "MINIMUM_RELATIVE_PERMEABILITY" : 0.0001 + }, + "Tables": {} + } + }] +} diff --git a/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/ProjectParameters.json b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/ProjectParameters.json new file mode 100644 index 000000000000..1135ac7fd704 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/ProjectParameters.json @@ -0,0 +1,226 @@ +{ + "problem_data": { + "problem_name": "test_phreatic", + "start_time": 0.0, + "end_time": 10.0, + "echo_level": 1, + "parallel_type": "OpenMP", + "number_of_threads": 1 + }, + "solver_settings": { + "solver_type": "U_Pw", + "model_part_name": "PorousDomain", + "domain_size": 2, + "model_import_settings": { + "input_type": "mdpa", + "input_filename": "test_phreatic" + }, + "material_import_settings": { + "materials_filename": "MaterialParameters.json" + }, + "time_stepping": { + "time_step": 1, + "max_delta_time_factor": 1000 + }, + "buffer_size": 2, + "echo_level": 1, + "clear_storage": false, + "compute_reactions": true, + "move_mesh_flag": false, + "reform_dofs_at_each_step": false, + "nodal_smoothing": false, + "block_builder": true, + "solution_type": "Quasi-static", + "scheme_type": "Newmark", + "reset_displacements": true, + "newmark_beta": 0.25, + "newmark_gamma": 0.5, + "newmark_theta": 0.5, + "rayleigh_m": 0.0, + "rayleigh_k": 0.0, + "strategy_type": "newton_raphson", + "convergence_criterion": "water_pressure_criterion", + "displacement_relative_tolerance": 0.0001, + "displacement_absolute_tolerance": 1e-09, + "residual_relative_tolerance": 0.0001, + "residual_absolute_tolerance": 1e-09, + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1e-09, + "min_iterations": 6, + "max_iterations": 15, + "number_cycles": 100, + "reduction_factor": 1, + "increase_factor": 1, + "desired_iterations": 4, + "max_radius_factor": 10.0, + "min_radius_factor": 0.1, + "calculate_reactions": true, + "max_line_search_iterations": 5, + "first_alpha_value": 0.5, + "second_alpha_value": 1.0, + "min_alpha": 0.1, + "max_alpha": 2.0, + "line_search_tolerance": 0.5, + "rotation_dofs": true, + "linear_solver_settings": { + "solver_type": "LinearSolversApplication.sparse_lu", + "scaling": true + }, + "problem_domain_sub_model_part_list": [ + "Soil_column" + ], + "processes_sub_model_part_list": [ + "Side_rollers", + "Base_Fixed", + "Initial_field_pressure", + "Fixed_base_pressure", + "Gravity_load" + ], + "body_domain_sub_model_part_list": [ + "Soil_column" + ] + }, + "output_processes": { + "gid_output": [ + { + "python_module": "gid_output_process", + "kratos_module": "KratosMultiphysics", + "process_name": "GiDOutputProcess", + "Parameters": { + "model_part_name": "PorousDomain.porous_computational_model_part", + "output_name": "test_phreatic", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag": "WriteElementsOnly", + "GiDPostMode": "GiD_PostAscii", + "MultiFileFlag": "SingleFile" + }, + "file_label": "step", + "output_control_type": "step", + "output_interval": 1, + "body_output": true, + "node_output": false, + "skin_output": false, + "plane_output": [], + "nodal_results": [ + "WATER_PRESSURE" + ], + "gauss_point_results": [ + "CAUCHY_STRESS_TENSOR", + "TOTAL_STRESS_TENSOR" + ] + }, + "point_data_configuration": [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "apply_vector_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name": "PorousDomain.Side_rollers", + "variable_name": "DISPLACEMENT", + "active": [ + true, + true, + true + ], + "is_fixed": [ + true, + false, + true + ], + "value": [ + 0.0, + 0.0, + 0.0 + ], + "table": [ + 0, + 0, + 0 + ] + } + }, + { + "python_module": "apply_vector_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name": "PorousDomain.Base_Fixed", + "variable_name": "DISPLACEMENT", + "active": [ + true, + true, + true + ], + "is_fixed": [ + true, + true, + true + ], + "value": [ + 0.0, + 0.0, + 0.0 + ], + "table": [ + 0, + 0, + 0 + ] + } + }, + { + "python_module": "apply_scalar_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name": "PorousDomain.Side_pressure", + "variable_name": "WATER_PRESSURE", + "is_fixed": true, + "fluid_pressure_type": "Hydrostatic", + "gravity_direction": 1, + "reference_coordinate": -25.0, + "table": 1, + "pressure_tension_cut_off": 0.0, + "specific_weight": 9.81 + } + } + ], + "loads_process_list": [ + { + "python_module": "apply_vector_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name": "PorousDomain.Gravity_load", + "variable_name": "VOLUME_ACCELERATION", + "active": [ + false, + true, + false + ], + "value": [ + 0.0, + -9.81, + 0.0 + ], + "table": [ + 0, + 0, + 0 + ] + } + } + ], + "auxiliar_process_list": [] + } +} \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/test_phreatic.mdpa b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/test_phreatic.mdpa new file mode 100644 index 000000000000..515e18b01fee --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_column_changing_waterlevel/test_side_pressure_boundaries/test_phreatic.mdpa @@ -0,0 +1,455 @@ +Begin Table 1 TIME WATER_PRESSURE + 0.0 0.0 + 1.0 0.0 + 2.0 -5.0 + 3.0 -10.0 + 4.0 -15.0 + 5.0 -10.0 + 6.0 -5.0 + 7.0 0.0 + 8.0 5.0 + 9.0 10.0 + 10.0 15.0 + 11.0 20.0 +End Table + +Begin Table 2 TIME VOLUME_ACCELERATION + 0.0 1.0 + 11.0 1.0 +End Table + +Begin Properties 1 +End Properties + + +Begin Nodes + 1 0.0000000000 0.0000000000 0.0000000000 + 2 0.0000000000 -2.5000000000 0.0000000000 + 3 0.0000000000 -5.0000000000 0.0000000000 + 4 5.0000000000 0.0000000000 0.0000000000 + 5 5.0000000000 -5.0000000000 0.0000000000 + 6 0.0000000000 -7.5000000000 0.0000000000 + 7 0.0000000000 -10.0000000000 0.0000000000 + 8 10.0000000000 0.0000000000 0.0000000000 + 9 10.0000000000 -2.5000000000 0.0000000000 + 10 5.0000000000 -10.0000000000 0.0000000000 + 11 10.0000000000 -5.0000000000 0.0000000000 + 12 10.0000000000 -7.5000000000 0.0000000000 + 13 0.0000000000 -12.5000000000 0.0000000000 + 14 10.0000000000 -10.0000000000 0.0000000000 + 15 0.0000000000 -15.0000000000 0.0000000000 + 16 5.0000000000 -15.0000000000 0.0000000000 + 17 10.0000000000 -12.5000000000 0.0000000000 + 18 0.0000000000 -17.5000000000 0.0000000000 + 19 10.0000000000 -15.0000000000 0.0000000000 + 20 0.0000000000 -20.0000000000 0.0000000000 + 21 10.0000000000 -17.5000000000 0.0000000000 + 22 5.0000000000 -20.0000000000 0.0000000000 + 23 10.0000000000 -20.0000000000 0.0000000000 + 24 0.0000000000 -22.5000000000 0.0000000000 + 25 10.0000000000 -22.5000000000 0.0000000000 + 26 0.0000000000 -25.0000000000 0.0000000000 + 27 5.0000000000 -25.0000000000 0.0000000000 + 28 10.0000000000 -25.0000000000 0.0000000000 + 29 0.0000000000 -27.5000000000 0.0000000000 + 30 10.0000000000 -27.5000000000 0.0000000000 + 31 0.0000000000 -30.0000000000 0.0000000000 + 32 5.0000000000 -30.0000000000 0.0000000000 + 33 10.0000000000 -30.0000000000 0.0000000000 + 34 0.0000000000 -32.5000000000 0.0000000000 + 35 10.0000000000 -32.5000000000 0.0000000000 + 36 0.0000000000 -35.0000000000 0.0000000000 + 37 5.0000000000 -35.0000000000 0.0000000000 + 38 10.0000000000 -35.0000000000 0.0000000000 + 39 0.0000000000 -37.5000000000 0.0000000000 + 40 10.0000000000 -37.5000000000 0.0000000000 + 41 0.0000000000 -40.0000000000 0.0000000000 + 42 5.0000000000 -40.0000000000 0.0000000000 + 43 10.0000000000 -40.0000000000 0.0000000000 + 44 0.0000000000 -42.5000000000 0.0000000000 + 45 10.0000000000 -42.5000000000 0.0000000000 + 46 0.0000000000 -45.0000000000 0.0000000000 + 47 5.0000000000 -45.0000000000 0.0000000000 + 48 10.0000000000 -45.0000000000 0.0000000000 + 49 0.0000000000 -47.5000000000 0.0000000000 + 50 10.0000000000 -47.5000000000 0.0000000000 + 51 0.0000000000 -50.0000000000 0.0000000000 + 52 5.0000000000 -50.0000000000 0.0000000000 + 53 10.0000000000 -50.0000000000 0.0000000000 +End Nodes + + +Begin Elements SmallStrainUPwDiffOrderElement2D8N + 1 1 53 48 46 51 50 47 49 52 + 2 1 48 43 41 46 45 42 44 47 + 3 1 43 38 36 41 40 37 39 42 + 4 1 38 33 31 36 35 32 34 37 + 5 1 33 28 26 31 30 27 29 32 + 6 1 28 23 20 26 25 22 24 27 + 7 1 23 19 15 20 21 16 18 22 + 8 1 19 14 7 15 17 10 13 16 + 9 1 14 11 3 7 12 5 6 10 + 10 1 11 8 1 3 9 4 2 5 +End Elements + + + +Begin SubModelPart Soil_column + Begin SubModelPartTables + End SubModelPartTables + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart Side_rollers + Begin SubModelPartTables + End SubModelPartTables + Begin SubModelPartNodes + 1 + 2 + 3 + 6 + 7 + 8 + 9 + 11 + 12 + 13 + 14 + 15 + 17 + 18 + 19 + 20 + 21 + 23 + 24 + 25 + 26 + 28 + 29 + 30 + 31 + 33 + 34 + 35 + 36 + 38 + 39 + 40 + 41 + 43 + 44 + 45 + 46 + 48 + 49 + 50 + 51 + 53 + End SubModelPartNodes + Begin SubModelPartElements + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart Side_pressure + Begin SubModelPartTables + 1 + End SubModelPartTables + Begin SubModelPartNodes + 1 + 2 + 3 + 6 + 7 + 8 + 9 + 11 + 12 + 13 + 14 + 15 + 17 + 18 + 19 + 20 + 21 + 23 + 24 + 25 + 26 + 28 + 29 + 30 + 31 + 33 + 34 + 35 + 36 + 38 + 39 + 40 + 41 + 43 + 44 + 45 + 46 + 48 + 49 + 50 + 51 + 53 + End SubModelPartNodes + Begin SubModelPartElements + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart Base_Fixed + Begin SubModelPartTables + End SubModelPartTables + Begin SubModelPartNodes + 51 + 52 + 53 + End SubModelPartNodes + Begin SubModelPartElements + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart Initial_field_pressure + Begin SubModelPartTables + 1 + End SubModelPartTables + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart Fixed_base_pressure + Begin SubModelPartTables + 1 + End SubModelPartTables + Begin SubModelPartNodes + 51 + 52 + 53 + End SubModelPartNodes + Begin SubModelPartElements + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart Gravity_load + Begin SubModelPartTables + 2 + End SubModelPartTables + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart diff --git a/applications/GeoMechanicsApplication/tests/test_helper.py b/applications/GeoMechanicsApplication/tests/test_helper.py index 2252fb6c413a..d9c6b5aa5222 100644 --- a/applications/GeoMechanicsApplication/tests/test_helper.py +++ b/applications/GeoMechanicsApplication/tests/test_helper.py @@ -11,7 +11,7 @@ def get_file_path(fileName): import os - return os.path.dirname(__file__) + "/" + fileName + return os.path.join(os.path.dirname(__file__), fileName) def run_kratos(file_path, model=None): @@ -412,3 +412,94 @@ def find_closest_index_greater_than_value(input_list, value): if value < list_value: return index return None + + +def read_numerical_results_from_post_res(project_path): + """ + Reads the numerical results from a Kratos "post.res" output file + :param project_path: path to the result file + :return: Puts the nodal and gauss point results into a dictionary sorted by result name for every time step + + """ + output_data = {"GaussPoints": {}, + "results": {}} + with open(project_path, "r") as result_file: + gauss_points_name = None + result_name = None + result_type = None + current_block_name = None + result_location = None + current_integration_point = 0 + for line in result_file: + line = line.strip() + if line.startswith("GaussPoints"): + current_block_name = "GaussPoints" + words = line.split() + gauss_points_name = words[1][1:-1] # strip off enclosing double quotes + output_data["GaussPoints"][gauss_points_name] = {} + continue + + if line == "End GaussPoints": + current_block_name = None + gauss_points_name = None + continue + + if line.startswith("Result"): + words = line.split() + result_name = words[1][1:-1] # strip off enclosing double quotes + if result_name not in output_data["results"]: + output_data["results"][result_name] = [] + result_type = words[4] + result_location = words[5] + this_result = {"time": float(words[3]), + "location": result_location, + "values": []} + output_data["results"][result_name].append(this_result) + if result_location == "OnGaussPoints": + current_integration_point = 0 + gauss_points_name = words[6][1:-1] # strip off enclosing double quotes + continue + + if line == "Values": + current_block_name = "Values" + continue + + if line == "End Values": + current_block_name = None + continue + + if current_block_name == "GaussPoints": + if line.startswith("Number Of Gauss Points:"): + pos = line.index(":") + num_gauss_points = int(line[pos+1:].strip()) + output_data["GaussPoints"][gauss_points_name]["size"] = num_gauss_points + continue + + if current_block_name == "Values": + words = line.split() + + if result_location == "OnNodes": + value = {"node": int(words[0])} + if result_type == "Scalar": + value["value"] = float(words[1]) + elif result_type == "Vector": + value["value"] = [float(x) for x in words[1:]] + + output_data["results"][result_name][-1]["values"].append(value) + elif result_location == "OnGaussPoints": + current_integration_point %= output_data["GaussPoints"][gauss_points_name]["size"] + current_integration_point += 1 + if current_integration_point == 1: + value = {"element": int(words[0]), + "value": []} + output_data["results"][result_name][-1]["values"].append(value) + words.pop(0) + + value = output_data["results"][result_name][-1]["values"][-1]["value"] + if result_type == "Scalar": + value.append(float(words[0])) + elif result_type == "Matrix": + value.append([float(v) for v in words]) + + continue + return output_data