From aa921daa506d9a6a15852c17d71d250f947bafae Mon Sep 17 00:00:00 2001 From: Chengshun Shang Date: Mon, 7 Oct 2024 17:02:31 +0200 Subject: [PATCH 1/6] add a measurement for density with periodic condition --- .../python_scripts/DEM_analysis_stage.py | 99 ++++++++++++++++++- 1 file changed, 98 insertions(+), 1 deletion(-) diff --git a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py index 94f647c3c0e6..931fa25d1f32 100644 --- a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py +++ b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py @@ -742,7 +742,7 @@ def OutputSingleTimeLoop(self): self.PrintResultsForGid(self.time) self.time_old_print = self.time - def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, center_z, type): + def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, center_z, type, domain_size=[1,1,1]): ''' This is a function to establish a sphere to measure local packing properties The type could be "porosity", "averaged_coordination_number", "fabric_tensor", "stress_tensor" or "strain" @@ -781,6 +781,103 @@ def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, c return measured_porosity + if type == "porosity_periodic": + + measure_sphere_volume = 4.0 / 3.0 * math.pi * radius * radius * radius + sphere_volume_inside_range = 0.0 + measured_porosity = 0.0 + + center_list = [] + center_list.append([center_x, center_y, center_z]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y + domain_size[1], center_z]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x + domain_size[0], center_y - domain_size[1], center_z]) + if center_y + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z + domain_size[2]]) + if center_y - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z - domain_size[2]]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x - domain_size[0], center_y + domain_size[1], center_z]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y - domain_size[1], center_z]) + if center_y + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z + domain_size[2]]) + if center_y - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z - domain_size[2]]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y + domain_size[1], center_z]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y + domain_size[1], center_z]) + if center_z + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z + domain_size[2]]) + if center_z - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z - domain_size[2]]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y - domain_size[1], center_z]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y - domain_size[1], center_z]) + if center_z + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z + domain_size[2]]) + if center_z - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z - domain_size[2]]) + if center_z + 0.5 * domain_size[2] < 2 * radius: + center_list.append([center_x, center_y, center_z + domain_size[2]]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z + domain_size[2]]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z + domain_size[2]]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z + domain_size[2]]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z + domain_size[2]]) + if center_z - 0.5 * domain_size[2] > -2 * radius: + center_list.append([center_x, center_y, center_z - domain_size[2]]) + if center_x + 0.5 * domain_size[0] < 2 * radius: + center_list.append([center_x + domain_size[0], center_y, center_z - domain_size[2]]) + if center_x - 0.5 * domain_size[0] > -2 * radius: + center_list.append([center_x - domain_size[0], center_y, center_z - domain_size[2]]) + if center_y + 0.5 * domain_size[1] < 2 * radius: + center_list.append([center_x, center_y + domain_size[1], center_z - domain_size[2]]) + if center_y - 0.5 * domain_size[1] > -2 * radius: + center_list.append([center_x, center_y - domain_size[1], center_z - domain_size[2]]) + + for center in center_list: + for node in self.spheres_model_part.Nodes: + + r = node.GetSolutionStepValue(RADIUS) + x = node.X + y = node.Y + z = node.Z + + center_to_sphere_distance = ((x - center[0])**2 + (y - center[1])**2 + (z - center[2])**2)**0.5 + + if center_to_sphere_distance < (radius - r): + + sphere_volume_inside_range += 4/3 * math.pi * r * r * r + + elif center_to_sphere_distance <= (radius + r): + + other_part_d = radius - (radius * radius + center_to_sphere_distance * center_to_sphere_distance - r * r) / (center_to_sphere_distance * 2) + + my_part_d = r - (r * r + center_to_sphere_distance * center_to_sphere_distance - radius * radius) / (center_to_sphere_distance * 2) + + cross_volume = math.pi * other_part_d * other_part_d * (radius - 1/3 * other_part_d) + math.pi * my_part_d * my_part_d * (r - 1/3 * my_part_d) + + sphere_volume_inside_range += cross_volume + + measured_porosity = 1.0 - (sphere_volume_inside_range / measure_sphere_volume) + + return measured_porosity + if type == "porosity_kdtree": measure_sphere_volume = 4.0 / 3.0 * math.pi * radius * radius * radius From 1b617096a24bbe778eea4fb2b448f05d48bb85c4 Mon Sep 17 00:00:00 2001 From: Chengshun Shang Date: Tue, 8 Oct 2024 23:20:43 +0200 Subject: [PATCH 2/6] first servo loading - not check yet --- .../DEMApplication/DEM_application.cpp | 2 + .../DEM_application_variables.h | 1 + .../custom_elements/spheric_particle.cpp | 5 + .../python_scripts/DEM_analysis_stage.py | 105 ++++++++++++++++-- .../dem_default_input_parameters.py | 7 +- .../python_scripts/sphere_strategy.py | 9 ++ 6 files changed, 116 insertions(+), 13 deletions(-) diff --git a/applications/DEMApplication/DEM_application.cpp b/applications/DEMApplication/DEM_application.cpp index 6cb1ef5dafba..55699529cbc6 100644 --- a/applications/DEMApplication/DEM_application.cpp +++ b/applications/DEMApplication/DEM_application.cpp @@ -153,6 +153,7 @@ KRATOS_CREATE_VARIABLE(bool, DELTA_OPTION) KRATOS_CREATE_VARIABLE(double, SKIN_SPHERE) KRATOS_CREATE_VARIABLE(int, PROPERTIES_ID) KRATOS_CREATE_VARIABLE(int, CONTACT_MESH_OPTION) +KRATOS_CREATE_VARIABLE(int, BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_CREATE_VARIABLE(double, MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_CREATE_VARIABLE(int, FAILURE_CRITERION_OPTION) KRATOS_CREATE_VARIABLE(int, CONCRETE_TEST_OPTION) @@ -660,6 +661,7 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(SKIN_SPHERE) KRATOS_REGISTER_VARIABLE(PROPERTIES_ID) KRATOS_REGISTER_VARIABLE(CONTACT_MESH_OPTION) + KRATOS_REGISTER_VARIABLE(BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_REGISTER_VARIABLE(MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_REGISTER_VARIABLE(FAILURE_CRITERION_OPTION) KRATOS_REGISTER_VARIABLE(CONCRETE_TEST_OPTION) diff --git a/applications/DEMApplication/DEM_application_variables.h b/applications/DEMApplication/DEM_application_variables.h index 2cf219a098cb..84ebdd5b358f 100644 --- a/applications/DEMApplication/DEM_application_variables.h +++ b/applications/DEMApplication/DEM_application_variables.h @@ -125,6 +125,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, PROPERTIES_ID) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CONTACT_MESH_OPTION) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, FAILURE_CRITERION_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CONCRETE_TEST_OPTION) diff --git a/applications/DEMApplication/custom_elements/spheric_particle.cpp b/applications/DEMApplication/custom_elements/spheric_particle.cpp index cbb3004337dd..abf006f168d9 100644 --- a/applications/DEMApplication/custom_elements/spheric_particle.cpp +++ b/applications/DEMApplication/custom_elements/spheric_particle.cpp @@ -894,6 +894,11 @@ void SphericParticle::ComputeBallToBallContactForceAndMoment(SphericParticle::Pa if ((i < (int)mNeighbourElements.size()) && this->Id() < neighbour_iterator_id) { CalculateOnContactElements(i, LocalContactForce, GlobalContactForce); } + } else if (r_process_info[BOUNDING_BOX_SERVO_LOADING_OPTION] == 1 && r_process_info[CONTACT_MESH_OPTION] == 1) { + unsigned int neighbour_iterator_id = data_buffer.mpOtherParticle->Id(); + if ((i < (int)mNeighbourElements.size()) && this->Id() < neighbour_iterator_id) { + CalculateOnContactElements(i, LocalContactForce, GlobalContactForce); + } } if (r_process_info[ENERGY_CALCULATION_OPTION]){ diff --git a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py index 931fa25d1f32..d6ac60c82122 100644 --- a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py +++ b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py @@ -530,32 +530,71 @@ def InitializeSolutionStep(self): if self.spheres_model_part.ProcessInfo[IMPOSED_Z_STRAIN_OPTION]: self.spheres_model_part.ProcessInfo.SetValue(IMPOSED_Z_STRAIN_VALUE, eval(self.DEM_parameters["ZStrainValue"].GetString())) + bounding_box_servo_loading_option = False + if "BoundingBoxServoLoadingOption" in self.DEM_parameters.keys(): + if self.DEM_parameters["BoundingBoxServoLoadingOption"].GetBool(): + bounding_box_servo_loading_option = True + if "BoundingBoxMoveOption" in self.DEM_parameters.keys(): if self.DEM_parameters["BoundingBoxMoveOption"].GetBool(): time_step = self.spheres_model_part.ProcessInfo[TIME_STEPS] - NStepSearch = self.DEM_parameters["NeighbourSearchFrequency"].GetInt() - if (time_step + 1) % NStepSearch == 0 and (time_step > 0): - self.UpdateSearchStartegyAndCPlusPlusStrategy() - self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor) + if bounding_box_servo_loading_option: + self.move_velocity = [0.0, 0.0, 0.0] + NStepSearch = self.DEM_parameters["BoundingBoxServoLoadingFrequency"].GetInt() + if (time_step + 1) % NStepSearch == 0 and (time_step > 0): + measured_global_stress = self.MeasureSphereForGettingGlobalStressTensor() + self.CalculateBoundingBoxMoveVelocity(measured_global_stress) + self.UpdateSearchStartegyAndCPlusPlusStrategy(self.move_velocity) + self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor) + else: + NStepSearch = self.DEM_parameters["NeighbourSearchFrequency"].GetInt() + if (time_step + 1) % NStepSearch == 0 and (time_step > 0): + move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetVector() + self.UpdateSearchStartegyAndCPlusPlusStrategy(move_velocity) + self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor) + + def CalculateBoundingBoxMoveVelocity(self, measured_global_stress): - def UpdateSearchStartegyAndCPlusPlusStrategy(self): + move_acceleration = [0.0, 0.0, 0.0] + servo_loading_stress = self.DEM_parameters["BoundingBoxServoLoadingStress"].GetVector() + servo_loading_factor = self.DEM_parameters["BoundingBoxServoLoadingFactor"].GetDouble() + + move_acceleration[0] = (servo_loading_stress[0] - measured_global_stress[0]) * servo_loading_factor + move_acceleration[1] = (servo_loading_stress[1] - measured_global_stress[1]) * servo_loading_factor + move_acceleration[2] = (servo_loading_stress[2] - measured_global_stress[2]) * servo_loading_factor + delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) - move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetDouble() + self.move_velocity[0] += delta_time * move_acceleration[0] + self.move_velocity[1] += delta_time * move_acceleration[1] + self.move_velocity[2] += delta_time * move_acceleration[2] + + bonuding_box_move_velocity_max = self.DEM_parameters["BoundingBoxMoveVelocityMax"].GetDouble() + if abs(self.move_velocity[0]) > bonuding_box_move_velocity_max: + self.move_velocity[0] = bonuding_box_move_velocity_max * np.sign(self.move_velocity[0]) + if abs(self.move_velocity[1]) > bonuding_box_move_velocity_max: + self.move_velocity[1] = bonuding_box_move_velocity_max * np.sign(self.move_velocity[1]) + if abs(self.move_velocity[2]) > bonuding_box_move_velocity_max: + self.move_velocity[2] = bonuding_box_move_velocity_max * np.sign(self.move_velocity[2]) + def UpdateSearchStartegyAndCPlusPlusStrategy(self, move_velocity): + + delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) + + # note: control_bool_vector = [MoveMinX, MoveMinY, MoveMinZ, MoveMaxX, MoveMaxY, MoveMaxZ] control_bool_vector = self.DEM_parameters["BoundingBoxMoveOptionDetail"].GetVector() if control_bool_vector[0]: - self.BoundingBoxMinX_update += delta_time * move_velocity + self.BoundingBoxMinX_update += delta_time * move_velocity[0] if control_bool_vector[1]: - self.BoundingBoxMinY_update += delta_time * move_velocity + self.BoundingBoxMinY_update += delta_time * move_velocity[1] if control_bool_vector[2]: - self.BoundingBoxMinZ_update += delta_time * move_velocity + self.BoundingBoxMinZ_update += delta_time * move_velocity[2] if control_bool_vector[3]: - self.BoundingBoxMaxX_update -= delta_time * move_velocity + self.BoundingBoxMaxX_update -= delta_time * move_velocity[0] if control_bool_vector[4]: - self.BoundingBoxMaxY_update -= delta_time * move_velocity + self.BoundingBoxMaxY_update -= delta_time * move_velocity[1] if control_bool_vector[5]: - self.BoundingBoxMaxZ_update -= delta_time * move_velocity + self.BoundingBoxMaxZ_update -= delta_time * move_velocity[2] self._GetSolver().search_strategy = OMP_DEMSearch(self.BoundingBoxMinX_update, self.BoundingBoxMinY_update, @@ -1204,6 +1243,48 @@ def MeasureSphereForGettingPackingProperties(self, radius, center_x, center_y, c if type == "strain": pass + def MeasureSphereForGettingGlobalStressTensor(self): + + if self.DEM_parameters["PostStressStrainOption"].GetBool() and self.DEM_parameters["ContactMeshOption"].GetBool(): + + bounding_box_volume = (self.BoundingBoxMaxX_update - self.BoundingBoxMinX_update) * \ + (self.BoundingBoxMaxY_update - self.BoundingBoxMinY_update) * \ + (self.BoundingBoxMaxZ_update - self.BoundingBoxMinZ_update) + + total_tensor = np.empty((3, 3)) + total_tensor[:] = 0.0 + stress_tensor_modulus = 0.0 + + for element in self.contact_model_part.Elements: + + x_0 = element.GetNode(0).X + x_1 = element.GetNode(1).X + y_0 = element.GetNode(0).Y + y_1 = element.GetNode(1).Y + z_0 = element.GetNode(0).Z + z_1 = element.GetNode(1).Z + + if abs(x_0 - x_1) > 0.5 * (self.BoundingBoxMaxX_update - self.BoundingBoxMinX_update) or \ + abs(y_0 - y_1) > 0.5 * (self.BoundingBoxMaxY_update - self.BoundingBoxMinY_update) or \ + abs(z_0 - z_1) > 0.5 * (self.BoundingBoxMaxZ_update - self.BoundingBoxMinZ_update): + continue + + local_contact_force_X = element.GetValue(GLOBAL_CONTACT_FORCE)[0] + local_contact_force_Y = element.GetValue(GLOBAL_CONTACT_FORCE)[1] + local_contact_force_Z = element.GetValue(GLOBAL_CONTACT_FORCE)[2] + contact_force_vector = np.array([local_contact_force_X , local_contact_force_Y, local_contact_force_Z]) + vector_l = np.array([x_0 - x_1 , y_0 - y_1, z_0 - z_1]) + tensor = np.outer(contact_force_vector, vector_l) + total_tensor += tensor + + total_tensor = total_tensor / bounding_box_volume + + return total_tensor + + else: + + raise Exception('The \"PostStressStrainOption\" and \"ContactMeshOption\" in the [ProjectParametersDEM.json] should be [True].') + def MeasureSphereForGettingRadialDistributionFunction(self, radius, center_x, center_y, center_z, delta_r, d_mean): min_reference_particle_to_center_distance = 1e10 diff --git a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py index 52476b1f7d73..b6fb8db50158 100644 --- a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py +++ b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py @@ -20,7 +20,12 @@ def GetDefaultInputParameters(): "BoundingBoxMinZ" : -10.0, "BoundingBoxMoveOption" : false, "BoundingBoxMoveOptionDetail" : [1, 1, 1, 1, 1, 1], - "BoundingBoxMoveVelocity" : 0.001, + "BoundingBoxMoveVelocity" : [0.01, 0.01, 0.01], + "BoundingBoxServoLoadingOption" : false, + "BoundingBoxServoLoadingStress" : [0.0, 0.0, 0.0], + "BoundingBoxServoLoadingFactor" : 1.0, + "BoundingBoxServoLoadingFrequency" : 50, + "BoundingBoxServoLoadingVelocityMax": 0.05, "dem_inlet_option" : true, "dem_inlets_settings" : {}, "seed" : 42, diff --git a/applications/DEMApplication/python_scripts/sphere_strategy.py b/applications/DEMApplication/python_scripts/sphere_strategy.py index d9ccadeebd4e..944298d52c30 100644 --- a/applications/DEMApplication/python_scripts/sphere_strategy.py +++ b/applications/DEMApplication/python_scripts/sphere_strategy.py @@ -81,6 +81,10 @@ def __init__(self, all_model_parts, creator_destructor, dem_fem_search, DEM_para self.contact_mesh_option = DEM_parameters["ContactMeshOption"].GetBool() self.automatic_bounding_box_option = DEM_parameters["AutomaticBoundingBoxOption"].GetBool() + self.bounding_box_servo_loading_option = 0 + if "BoundingBoxServoLoadingOption" in DEM_parameters.keys(): + self.bounding_box_servo_loading_option = DEM_parameters["BoundingBoxServoLoadingOption"].GetBool() + self.delta_option = DEM_parameters["DeltaOption"].GetString() #TODO: this is not an option (bool) let's change the name to something including 'type' self.search_increment = 0.0 @@ -322,6 +326,11 @@ def SetVariablesAndOptions(self): else: self.spheres_model_part.ProcessInfo.SetValue(CONTACT_MESH_OPTION, 0) + if self.servo_loading_option: + self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_SERVO_LOADING_OPTION, 1) + else: + self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_SERVO_LOADING_OPTION, 0) + # PRINTING VARIABLES self.spheres_model_part.ProcessInfo.SetValue(PRINT_EXPORT_ID, self.print_export_id) From 71860d9bc4b8c4dd5a99128367b4cb795219aeff Mon Sep 17 00:00:00 2001 From: Chengshun Shang Date: Thu, 10 Oct 2024 10:31:45 +0200 Subject: [PATCH 3/6] a stable servo loading --- .../DEMApplication/DEM_application.cpp | 2 + .../DEM_application_variables.h | 1 + .../custom_elements/spheric_particle.cpp | 2 +- .../custom_python/DEM_python_application.cpp | 2 + .../strategies/explicit_solver_strategy.cpp | 3 + .../python_scripts/DEM_analysis_stage.py | 100 +++++++++++++----- .../python_scripts/DEM_procedures.py | 17 ++- .../dem_default_input_parameters.py | 13 ++- .../python_scripts/sphere_strategy.py | 2 +- 9 files changed, 98 insertions(+), 44 deletions(-) diff --git a/applications/DEMApplication/DEM_application.cpp b/applications/DEMApplication/DEM_application.cpp index 55699529cbc6..6a05ec219e23 100644 --- a/applications/DEMApplication/DEM_application.cpp +++ b/applications/DEMApplication/DEM_application.cpp @@ -101,6 +101,7 @@ KRATOS_CREATE_VARIABLE(int, ROTATION_OPTION) KRATOS_CREATE_VARIABLE(int, VIRTUAL_MASS_OPTION) KRATOS_CREATE_VARIABLE(int, SEARCH_CONTROL) KRATOS_CREATE_VARIABLE(bool, IS_TIME_TO_PRINT) +KRATOS_CREATE_VARIABLE(bool, IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_CREATE_VARIABLE(double, COORDINATION_NUMBER) KRATOS_CREATE_VARIABLE(double, CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_CREATE_VARIABLE(double, MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) @@ -613,6 +614,7 @@ void KratosDEMApplication::Register() { KRATOS_REGISTER_VARIABLE(VIRTUAL_MASS_OPTION) KRATOS_REGISTER_VARIABLE(SEARCH_CONTROL) KRATOS_REGISTER_VARIABLE(IS_TIME_TO_PRINT) + KRATOS_REGISTER_VARIABLE(IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_REGISTER_VARIABLE(COORDINATION_NUMBER) KRATOS_REGISTER_VARIABLE(CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_REGISTER_VARIABLE(MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) diff --git a/applications/DEMApplication/DEM_application_variables.h b/applications/DEMApplication/DEM_application_variables.h index 84ebdd5b358f..a45f1d4753a1 100644 --- a/applications/DEMApplication/DEM_application_variables.h +++ b/applications/DEMApplication/DEM_application_variables.h @@ -68,6 +68,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, VIRTUAL_MASS_OPTION) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, SEARCH_CONTROL) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, IS_TIME_TO_PRINT) + KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, COORDINATION_NUMBER) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) diff --git a/applications/DEMApplication/custom_elements/spheric_particle.cpp b/applications/DEMApplication/custom_elements/spheric_particle.cpp index abf006f168d9..ce7366a02521 100644 --- a/applications/DEMApplication/custom_elements/spheric_particle.cpp +++ b/applications/DEMApplication/custom_elements/spheric_particle.cpp @@ -894,7 +894,7 @@ void SphericParticle::ComputeBallToBallContactForceAndMoment(SphericParticle::Pa if ((i < (int)mNeighbourElements.size()) && this->Id() < neighbour_iterator_id) { CalculateOnContactElements(i, LocalContactForce, GlobalContactForce); } - } else if (r_process_info[BOUNDING_BOX_SERVO_LOADING_OPTION] == 1 && r_process_info[CONTACT_MESH_OPTION] == 1) { + } else if (r_process_info[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] && r_process_info[CONTACT_MESH_OPTION] == 1) { unsigned int neighbour_iterator_id = data_buffer.mpOtherParticle->Id(); if ((i < (int)mNeighbourElements.size()) && this->Id() < neighbour_iterator_id) { CalculateOnContactElements(i, LocalContactForce, GlobalContactForce); diff --git a/applications/DEMApplication/custom_python/DEM_python_application.cpp b/applications/DEMApplication/custom_python/DEM_python_application.cpp index 44bbc81b2c0f..0c3fa1f75fcc 100644 --- a/applications/DEMApplication/custom_python/DEM_python_application.cpp +++ b/applications/DEMApplication/custom_python/DEM_python_application.cpp @@ -71,6 +71,7 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, AUTOMATIC_SKIN_COMPUTATION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SKIN_FACTOR_RADIUS) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, IS_TIME_TO_PRINT) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, IS_TIME_TO_UPDATE_CONTACT_ELEMENT) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, COORDINATION_NUMBER) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CONTINUUM_SEARCH_RADIUS_AMPLIFICATION_FACTOR) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MAX_AMPLIFICATION_RATIO_OF_THE_SEARCH_RADIUS) @@ -125,6 +126,7 @@ PYBIND11_MODULE(KratosDEMApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, COHESIVE_GROUP) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PROPERTIES_ID) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CONTACT_MESH_OPTION) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, BOUNDING_BOX_SERVO_LOADING_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, MAX_NUMBER_OF_INTACT_BONDS_TO_CONSIDER_A_SPHERE_BROKEN) //KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FAILURE_CRITERION_OPTION) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, CONCRETE_TEST_OPTION) diff --git a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp index 74c115e07a96..154bd71697be 100644 --- a/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp +++ b/applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp @@ -459,6 +459,9 @@ namespace Kratos { if (is_time_to_print_results && r_process_info[CONTACT_MESH_OPTION] == 1) { CreateContactElements(); InitializeContactElements(); + }else if (r_process_info[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] && r_process_info[CONTACT_MESH_OPTION] == 1) { + CreateContactElements(); + InitializeContactElements(); } //RebuildPropertiesProxyPointers(mListOfSphericParticles); diff --git a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py index d6ac60c82122..8614ea26f5e7 100644 --- a/applications/DEMApplication/python_scripts/DEM_analysis_stage.py +++ b/applications/DEMApplication/python_scripts/DEM_analysis_stage.py @@ -245,6 +245,7 @@ def FillAnalyticSubModelPartsWithNewParticles(self): def Initialize(self): self.time = 0.0 self.time_old_print = 0.0 + self.time_old_update_contact_element = 0.0 self.ReadModelParts() @@ -309,6 +310,13 @@ def Initialize(self): self.BoundingBoxMaxY_update = self.DEM_parameters["BoundingBoxMaxY"].GetDouble() self.BoundingBoxMaxZ_update = self.DEM_parameters["BoundingBoxMaxZ"].GetDouble() + self.bounding_box_servo_loading_option = False + if "BoundingBoxServoLoadingOption" in self.DEM_parameters.keys(): + if self.DEM_parameters["BoundingBoxServoLoadingOption"].GetBool(): + self.bounding_box_servo_loading_option = True + + self.bounding_box_move_velocity = [0.0, 0.0, 0.0] + def SetMaterials(self): self.ReadMaterialsFile() @@ -498,6 +506,10 @@ def RunAnalytics(self, time): def IsTimeToPrintPostProcess(self): return self.do_print_results_option and self.DEM_parameters["OutputTimeStep"].GetDouble() - (self.time - self.time_old_print) < 1e-2 * self._GetSolver().dt + def IsTimeToUpdateContactElementForServo(self): + return self.bounding_box_servo_loading_option and self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFrequency"].GetInt() * \ + self._GetSolver().dt - (self.time - self.time_old_update_contact_element) < 1e-2 * self._GetSolver().dt + def PrintResults(self): #### GiD IO ########################################## if self.IsTimeToPrintPostProcess(): @@ -523,6 +535,8 @@ def InitializeSolutionStep(self): self.FillAnalyticSubModelPartsWithNewParticles() if self.DEM_parameters["ContactMeshOption"].GetBool(): self.UpdateIsTimeToPrintInModelParts(self.IsTimeToPrintPostProcess()) + if self.bounding_box_servo_loading_option: + self.UpdateIsTimeToUpdateContactElementForServo(self.IsTimeToUpdateContactElementForServo()) if self.DEM_parameters["Dimension"].GetInt() == 2: self.spheres_model_part.ProcessInfo[IMPOSED_Z_STRAIN_OPTION] = self.DEM_parameters["ImposeZStrainIn2DOption"].GetBool() @@ -539,43 +553,70 @@ def InitializeSolutionStep(self): if self.DEM_parameters["BoundingBoxMoveOption"].GetBool(): time_step = self.spheres_model_part.ProcessInfo[TIME_STEPS] if bounding_box_servo_loading_option: - self.move_velocity = [0.0, 0.0, 0.0] - NStepSearch = self.DEM_parameters["BoundingBoxServoLoadingFrequency"].GetInt() + NStepSearch = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFrequency"].GetInt() if (time_step + 1) % NStepSearch == 0 and (time_step > 0): measured_global_stress = self.MeasureSphereForGettingGlobalStressTensor() self.CalculateBoundingBoxMoveVelocity(measured_global_stress) - self.UpdateSearchStartegyAndCPlusPlusStrategy(self.move_velocity) - self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor) + self.UpdateSearchStartegyAndCPlusPlusStrategy(self.bounding_box_move_velocity) + self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor, self.bounding_box_move_velocity) else: NStepSearch = self.DEM_parameters["NeighbourSearchFrequency"].GetInt() if (time_step + 1) % NStepSearch == 0 and (time_step > 0): - move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetVector() - self.UpdateSearchStartegyAndCPlusPlusStrategy(move_velocity) - self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor) + bounding_box_move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetVector() + self.UpdateSearchStartegyAndCPlusPlusStrategy(bounding_box_move_velocity) + self.procedures.UpdateBoundingBox(self.spheres_model_part, self.creator_destructor, bounding_box_move_velocity) def CalculateBoundingBoxMoveVelocity(self, measured_global_stress): - move_acceleration = [0.0, 0.0, 0.0] + # note: servo_loading_type = "isotropic" or "anisotropic" + servo_loading_type = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingType"].GetString() + + if servo_loading_type == "isotropic": + self.CalculateBoundingBoxMoveVelocityIsotropic(measured_global_stress) + elif servo_loading_type == "anisotropic": + self.CalculateBoundingBoxMoveVelocityAnisotropic(measured_global_stress) + + def CalculateBoundingBoxMoveVelocityIsotropic(self, measured_global_stress): + + servo_loading_stress = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingStress"].GetVector() + servo_loading_factor = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFactor"].GetDouble() + d50 = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["MeanParticleDiameterD50"].GetDouble() + Youngs_modulus = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["ParticleYoungsModulus"].GetDouble() - servo_loading_stress = self.DEM_parameters["BoundingBoxServoLoadingStress"].GetVector() - servo_loading_factor = self.DEM_parameters["BoundingBoxServoLoadingFactor"].GetDouble() + servo_loading_stress_mean = (servo_loading_stress[0] + servo_loading_stress[1] + servo_loading_stress[2]) / 3.0 + measured_global_stress_mean = (measured_global_stress[0][0] + measured_global_stress[1][1] + measured_global_stress[2][2]) / 3.0 - move_acceleration[0] = (servo_loading_stress[0] - measured_global_stress[0]) * servo_loading_factor - move_acceleration[1] = (servo_loading_stress[1] - measured_global_stress[1]) * servo_loading_factor - move_acceleration[2] = (servo_loading_stress[2] - measured_global_stress[2]) * servo_loading_factor + delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) + servo_loading_coefficient = servo_loading_factor * d50 / (delta_time * Youngs_modulus) + self.bounding_box_move_velocity[0] = (servo_loading_stress_mean - measured_global_stress_mean) * servo_loading_coefficient + + bonuding_box_serva_loading_velocity_max = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingVelocityMax"].GetDouble() + if abs(self.bounding_box_move_velocity[0]) > bonuding_box_serva_loading_velocity_max: + self.bounding_box_move_velocity[0] = bonuding_box_serva_loading_velocity_max * np.sign(self.bounding_box_move_velocity[0]) + + self.bounding_box_move_velocity[1] = self.bounding_box_move_velocity[0] + self.bounding_box_move_velocity[2] = self.bounding_box_move_velocity[0] + + def CalculateBoundingBoxMoveVelocityAnisotropic(self, measured_global_stress): + + servo_loading_stress = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingStress"].GetVector() + servo_loading_factor = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingFactor"].GetDouble() + d50 = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["MeanParticleDiameterD50"].GetDouble() + Youngs_modulus = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["ParticleYoungsModulus"].GetDouble() delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) - self.move_velocity[0] += delta_time * move_acceleration[0] - self.move_velocity[1] += delta_time * move_acceleration[1] - self.move_velocity[2] += delta_time * move_acceleration[2] - - bonuding_box_move_velocity_max = self.DEM_parameters["BoundingBoxMoveVelocityMax"].GetDouble() - if abs(self.move_velocity[0]) > bonuding_box_move_velocity_max: - self.move_velocity[0] = bonuding_box_move_velocity_max * np.sign(self.move_velocity[0]) - if abs(self.move_velocity[1]) > bonuding_box_move_velocity_max: - self.move_velocity[1] = bonuding_box_move_velocity_max * np.sign(self.move_velocity[1]) - if abs(self.move_velocity[2]) > bonuding_box_move_velocity_max: - self.move_velocity[2] = bonuding_box_move_velocity_max * np.sign(self.move_velocity[2]) + servo_loading_coefficient = servo_loading_factor * d50 / (delta_time * Youngs_modulus) + self.bounding_box_move_velocity[0] = (servo_loading_stress[0] - measured_global_stress[0][0]) * servo_loading_coefficient + self.bounding_box_move_velocity[1] = (servo_loading_stress[1] - measured_global_stress[1][1]) * servo_loading_coefficient + self.bounding_box_move_velocity[2] = (servo_loading_stress[2] - measured_global_stress[2][2]) * servo_loading_coefficient + + bonuding_box_move_velocity_max = self.DEM_parameters["BoundingBoxServoLoadingSettings"]["BoundingBoxServoLoadingVelocityMax"].GetDouble() + if abs(self.bounding_box_move_velocity[0]) > bonuding_box_move_velocity_max: + self.bounding_box_move_velocity[0] = bonuding_box_move_velocity_max * np.sign(self.bounding_box_move_velocity[0]) + if abs(self.bounding_box_move_velocity[1]) > bonuding_box_move_velocity_max: + self.bounding_box_move_velocity[1] = bonuding_box_move_velocity_max * np.sign(self.bounding_box_move_velocity[1]) + if abs(self.bounding_box_move_velocity[2]) > bonuding_box_move_velocity_max: + self.bounding_box_move_velocity[2] = bonuding_box_move_velocity_max * np.sign(self.bounding_box_move_velocity[2]) def UpdateSearchStartegyAndCPlusPlusStrategy(self, move_velocity): @@ -610,6 +651,13 @@ def UpdateIsTimeToPrintInModelParts(self, is_time_to_print): self.UpdateIsTimeToPrintInOneModelPart(self.dem_inlet_model_part, is_time_to_print) self.UpdateIsTimeToPrintInOneModelPart(self.rigid_face_model_part, is_time_to_print) + def UpdateIsTimeToUpdateContactElementForServo(self, is_time_to_update_contact_element): + self.spheres_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.cluster_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.dem_inlet_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.rigid_face_model_part.ProcessInfo[IS_TIME_TO_UPDATE_CONTACT_ELEMENT] = is_time_to_update_contact_element + self.time_old_update_contact_element = self.time + def UpdateIsTimeToPrintInOneModelPart(self, model_part, is_time_to_print): model_part.ProcessInfo[IS_TIME_TO_PRINT] = is_time_to_print @@ -1253,7 +1301,6 @@ def MeasureSphereForGettingGlobalStressTensor(self): total_tensor = np.empty((3, 3)) total_tensor[:] = 0.0 - stress_tensor_modulus = 0.0 for element in self.contact_model_part.Elements: @@ -1821,9 +1868,6 @@ def MeasureCubicForGettingPackingProperties(self, side_length, center_x, center_ if total_contact_number: averaged_contact_force_modulus_square = total_contact_force_tensor_modulus_square / total_contact_number - print(averaged_total_particle_force_tensor_modulus_square) - print(averaged_contact_force_modulus_square) - if averaged_contact_force_modulus_square: return (averaged_total_particle_force_tensor_modulus_square / averaged_contact_force_modulus_square)**0.5 else: diff --git a/applications/DEMApplication/python_scripts/DEM_procedures.py b/applications/DEMApplication/python_scripts/DEM_procedures.py index ba5e3df9dec8..5201c4266098 100644 --- a/applications/DEMApplication/python_scripts/DEM_procedures.py +++ b/applications/DEMApplication/python_scripts/DEM_procedures.py @@ -826,28 +826,25 @@ def SetBoundingBox(self, spheres_model_part, clusters_model_part, rigid_faces_mo creator_destructor.SetHighNode(b_box_high) creator_destructor.CalculateSurroundingBoundingBox(spheres_model_part, clusters_model_part, rigid_faces_model_part, dem_inlet_model_part, self.bounding_box_enlargement_factor, self.automatic_bounding_box_OPTION) - def UpdateBoundingBox(self, spheres_model_part, creator_destructor): + def UpdateBoundingBox(self, spheres_model_part, creator_destructor, move_velocity): delta_time = spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) - #NStepSearch = self.DEM_parameters["NeighbourSearchFrequency"].GetInt() - #delta_time_real = delta_time * NStepSearch - move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetDouble() b_box_low = Array3() b_box_high = Array3() control_bool_vector = self.DEM_parameters["BoundingBoxMoveOptionDetail"].GetVector() if control_bool_vector[0]: - self.b_box_minX += delta_time * move_velocity + self.b_box_minX += delta_time * move_velocity[0] if control_bool_vector[1]: - self.b_box_minY += delta_time * move_velocity + self.b_box_minY += delta_time * move_velocity[1] if control_bool_vector[2]: - self.b_box_minZ += delta_time * move_velocity + self.b_box_minZ += delta_time * move_velocity[2] if control_bool_vector[3]: - self.b_box_maxX -= delta_time * move_velocity + self.b_box_maxX -= delta_time * move_velocity[0] if control_bool_vector[4]: - self.b_box_maxY -= delta_time * move_velocity + self.b_box_maxY -= delta_time * move_velocity[1] if control_bool_vector[5]: - self.b_box_maxZ -= delta_time * move_velocity + self.b_box_maxZ -= delta_time * move_velocity[2] b_box_low[0] = self.b_box_minX b_box_low[1] = self.b_box_minY b_box_low[2] = self.b_box_minZ diff --git a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py index b6fb8db50158..e35bae680caa 100644 --- a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py +++ b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py @@ -22,10 +22,15 @@ def GetDefaultInputParameters(): "BoundingBoxMoveOptionDetail" : [1, 1, 1, 1, 1, 1], "BoundingBoxMoveVelocity" : [0.01, 0.01, 0.01], "BoundingBoxServoLoadingOption" : false, - "BoundingBoxServoLoadingStress" : [0.0, 0.0, 0.0], - "BoundingBoxServoLoadingFactor" : 1.0, - "BoundingBoxServoLoadingFrequency" : 50, - "BoundingBoxServoLoadingVelocityMax": 0.05, + "BoundingBoxServoLoadingSettings":{ + "BoundingBoxServoLoadingType" : "isotropic", + "BoundingBoxServoLoadingStress" : [0.0, 0.0, 0.0], + "BoundingBoxServoLoadingFactor" : 0.8, + "BoundingBoxServoLoadingFrequency" : 50, + "BoundingBoxServoLoadingVelocityMax": 0.05, + "MeanParticleDiameterD50" : 0.1, + "ParticleYoungsModulus" : 1e9 + }, "dem_inlet_option" : true, "dem_inlets_settings" : {}, "seed" : 42, diff --git a/applications/DEMApplication/python_scripts/sphere_strategy.py b/applications/DEMApplication/python_scripts/sphere_strategy.py index 944298d52c30..3d8599bcbfe5 100644 --- a/applications/DEMApplication/python_scripts/sphere_strategy.py +++ b/applications/DEMApplication/python_scripts/sphere_strategy.py @@ -326,7 +326,7 @@ def SetVariablesAndOptions(self): else: self.spheres_model_part.ProcessInfo.SetValue(CONTACT_MESH_OPTION, 0) - if self.servo_loading_option: + if self.bounding_box_servo_loading_option: self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_SERVO_LOADING_OPTION, 1) else: self.spheres_model_part.ProcessInfo.SetValue(BOUNDING_BOX_SERVO_LOADING_OPTION, 0) From 9f3383313508ace92474b95780f0496adaf31b23 Mon Sep 17 00:00:00 2001 From: Chengshun Shang Date: Thu, 10 Oct 2024 10:37:07 +0200 Subject: [PATCH 4/6] fix the test case of moving boundaries --- .../ProjectParametersDEM.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json b/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json index f3762b36437a..7f94c9e90c2e 100644 --- a/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json +++ b/applications/DEMApplication/tests/moving_periodic_boundary_tests_files/ProjectParametersDEM.json @@ -13,7 +13,7 @@ "BoundingBoxMinY" : 0.0, "BoundingBoxMinZ" : -0.0015, "BoundingBoxMoveOption" : true, - "BoundingBoxMoveVelocity" : 0.001, + "BoundingBoxMoveVelocity" : [0.001, 0.001, 0.001], "dem_inlet_option" : false, "GravityX" : 0.0, "GravityY" : 0.0, From aa0ad5ab98fb8910074c7146532ad2c67bbc0ba3 Mon Sep 17 00:00:00 2001 From: Chengshun Shang Date: Fri, 25 Oct 2024 18:12:11 +0200 Subject: [PATCH 5/6] add test case for servo control --- .../MaterialsDEM.json | 29 +++++++ .../ProjectParametersDEM.json | 87 +++++++++++++++++++ .../servo_control_tests_files/testSCDEM.mdpa | 55 ++++++++++++ .../testSCDEM_FEM_boundary.mdpa | 6 ++ .../testSC_Graphs/EnergyPlot.grf | 2 + .../tests/test_DEMApplication.py | 2 + .../tests/test_servo_control.py | 52 +++++++++++ 7 files changed, 233 insertions(+) create mode 100644 applications/DEMApplication/tests/servo_control_tests_files/MaterialsDEM.json create mode 100644 applications/DEMApplication/tests/servo_control_tests_files/ProjectParametersDEM.json create mode 100644 applications/DEMApplication/tests/servo_control_tests_files/testSCDEM.mdpa create mode 100644 applications/DEMApplication/tests/servo_control_tests_files/testSCDEM_FEM_boundary.mdpa create mode 100644 applications/DEMApplication/tests/servo_control_tests_files/testSC_Graphs/EnergyPlot.grf create mode 100644 applications/DEMApplication/tests/test_servo_control.py diff --git a/applications/DEMApplication/tests/servo_control_tests_files/MaterialsDEM.json b/applications/DEMApplication/tests/servo_control_tests_files/MaterialsDEM.json new file mode 100644 index 000000000000..2573997bcfd9 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/MaterialsDEM.json @@ -0,0 +1,29 @@ +{ + "materials" : [{ + "material_name" : "DEM-DefaultMaterial", + "material_id" : 1, + "Variables" : { + "PARTICLE_DENSITY" : 2650.0, + "YOUNG_MODULUS" : 3.5e7, + "POISSON_RATIO" : 0.2, + "PARTICLE_SPHERICITY" : 1.0 + } + }], + "material_relations" : [{ + "material_names_list" : ["DEM-DefaultMaterial","DEM-DefaultMaterial"], + "material_ids_list" : [1,1], + "Variables" : { + "DEM_DISCONTINUUM_CONSTITUTIVE_LAW_NAME" : "DEM_D_Hertz_viscous_Coulomb", + "PARTICLE_COHESION" : 0.0, + "STATIC_FRICTION" : 0.6, + "DYNAMIC_FRICTION" : 0.6, + "FRICTION_DECAY" : 500, + "COEFFICIENT_OF_RESTITUTION" : 0.01, + "ROLLING_FRICTION" : 0.01, + "ROLLING_FRICTION_WITH_WALLS" : 0.01, + "DEM_ROLLING_FRICTION_MODEL_NAME" : "DEMRollingFrictionModelConstantTorque" + } + }], + "material_assignation_table" : [["SpheresPart.DEMParts_Body","DEM-DefaultMaterial"]] +} + diff --git a/applications/DEMApplication/tests/servo_control_tests_files/ProjectParametersDEM.json b/applications/DEMApplication/tests/servo_control_tests_files/ProjectParametersDEM.json new file mode 100644 index 000000000000..b6a68a9500be --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/ProjectParametersDEM.json @@ -0,0 +1,87 @@ +{ + "Dimension" : 3, + "PeriodicDomainOption" : true, + "BoundingBoxOption" : true, + "AutomaticBoundingBoxOption" : false, + "BoundingBoxEnlargementFactor" : 1.0, + "BoundingBoxStartTime" : 0.0, + "BoundingBoxStopTime" : 1000.0, + "BoundingBoxMaxX" : 0.2, + "BoundingBoxMaxY" : 0.241421356237, + "BoundingBoxMaxZ" : 0.2, + "BoundingBoxMinX" : -0.2, + "BoundingBoxMinY" : -0.241421356237, + "BoundingBoxMinZ" : -0.2, + "BoundingBoxMoveOption" : true, + "BoundingBoxMoveVelocity" : [1.0, 1.0, 1.0], + "BoundingBoxMoveOptionDetail" : [1, 1, 1, 1, 1, 1], + "BoundingBoxServoLoadingOption" : true, + "BoundingBoxServoLoadingSettings":{ + "BoundingBoxServoLoadingType" : "isotropic", + "BoundingBoxServoLoadingStress" : [5e3, 5e3, 5e3], + "BoundingBoxServoLoadingFactor" : 5, + "BoundingBoxServoLoadingFrequency" : 10, + "BoundingBoxServoLoadingVelocityMax": 1.0, + "MeanParticleDiameterD50" : 2.1e-4, + "ParticleYoungsModulus" : 3.5e10 + }, + "dem_inlet_option" : false, + "dem_inlets_settings" : {}, + "GravityX" : 0.0, + "GravityY" : 0.0, + "GravityZ" : 0.0, + "RotationOption" : true, + "CleanIndentationsOption" : false, + "RadiusExpansionOption" : false, + "RadiusExpansionRate" : 1000.0, + "RadiusMultiplierMax" : 2.0, + "RadiusExpansionRateChangeOption": false, + "RadiusExpansionAcceleration" : -4e5, + "RadiusExpansionRateMin" : 100.0, + "EnergyCalculationOption" : true, + "ComputeStressTensorOption" : true, + "solver_settings" : { + "RemoveBallsInitiallyTouchingWalls" : false, + "strategy" : "sphere_strategy", + "material_import_settings" : { + "materials_filename" : "MaterialsDEM.json" + } + }, + "VirtualMassCoefficient" : 1.0, + "RollingFrictionOption" : true, + "GlobalDamping" : 0.0, + "GlobalViscousDamping" : 0.0, + "ContactMeshOption" : true, + "OutputFileType" : "Ascii", + "Multifile" : "multiple_files", + "ElementType" : "SphericPartDEMElement3D", + "TranslationalIntegrationScheme" : "Symplectic_Euler", + "RotationalIntegrationScheme" : "Direct_Integration", + "MaxTimeStep" : 2e-8, + "FinalTime" : 0.00006, + "NeighbourSearchFrequency" : 50, + "SearchTolerance" : 0.0000002, + "GraphExportFreq" : 0.00004, + "VelTrapGraphExportFreq" : 0.00004, + "OutputTimeStep" : 0.00004, + "PostBoundingBox" : true, + "PostLocalContactForce" : false, + "PostDisplacement" : true, + "PostRadius" : true, + "PostVelocity" : true, + "PostAngularVelocity" : false, + "PostElasticForces" : false, + "PostContactForces" : false, + "PostRigidElementForces" : false, + "PostStressStrainOption" : true, + "PostTangentialElasticForces" : false, + "PostTotalForces" : false, + "PostPressure" : false, + "PostShearStress" : false, + "PostSkinSphere" : false, + "PostNonDimensionalVolumeWear" : false, + "PostParticleMoment" : false, + "PostEulerAngles" : false, + "PostRollingResistanceMoment" : false, + "problem_name" : "testSC" +} diff --git a/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM.mdpa b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM.mdpa new file mode 100644 index 000000000000..84ae1326d776 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM.mdpa @@ -0,0 +1,55 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes +1 -0.1 0.0 0.1 +2 0.1 0.0 0.1 +3 -0.1 0.0 -0.1 +4 0.1 0.0 -0.1 +5 0.0 0.141421356237 0.0 +6 0.0 -0.141421356237 0.0 +End Nodes + + +Begin Elements SphericParticle3D// GUI group identifier: Body +1 0 1 +2 0 2 +3 0 3 +4 0 4 +5 0 5 +6 0 6 +End Elements + +Begin NodalData RADIUS // GUI group identifier: Body + 1 0 0.1 + 2 0 0.1 + 3 0 0.1 + 4 0 0.1 + 5 0 0.1 + 6 0 0.1 +End NodalData + +Begin SubModelPart DEMParts_Body // Group Body // Subtree DEMParts +Begin SubModelPartNodes +1 +2 +3 +4 +5 +6 +End SubModelPartNodes +Begin SubModelPartElements +1 +2 +3 +4 +5 +6 +End SubModelPartElements +Begin SubModelPartConditions +End SubModelPartConditions +End SubModelPart diff --git a/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM_FEM_boundary.mdpa b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM_FEM_boundary.mdpa new file mode 100644 index 000000000000..c3ac404e2a68 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/testSCDEM_FEM_boundary.mdpa @@ -0,0 +1,6 @@ +Begin ModelPartData + // VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties diff --git a/applications/DEMApplication/tests/servo_control_tests_files/testSC_Graphs/EnergyPlot.grf b/applications/DEMApplication/tests/servo_control_tests_files/testSC_Graphs/EnergyPlot.grf new file mode 100644 index 000000000000..355c1f541d85 --- /dev/null +++ b/applications/DEMApplication/tests/servo_control_tests_files/testSC_Graphs/EnergyPlot.grf @@ -0,0 +1,2 @@ + Time Trans kinematic energy Rot kinematic energy Kinematic energy Gravitational energy Elastic energy Frictional energy Viscodamping energy Rolling resistance energy Total energy +4.002e-05 6.17221e-19 1.64829e-28 6.17221e-19 0 9.58515e-12 0 1.26916e-21 2.09691e-25 9.58515e-12 diff --git a/applications/DEMApplication/tests/test_DEMApplication.py b/applications/DEMApplication/tests/test_DEMApplication.py index 987099fd52ef..579f6a9077b6 100644 --- a/applications/DEMApplication/tests/test_DEMApplication.py +++ b/applications/DEMApplication/tests/test_DEMApplication.py @@ -35,6 +35,7 @@ import test_dem_3d_parallel_bond_model import test_dem_3d_smooth_joint_model import test_moving_periodic_boundary +import test_servo_control import sys sys.path.append('DEM3D_chung_ooi_tests/test1_data') sys.path.append('DEM3D_chung_ooi_tests/test2_data') @@ -107,6 +108,7 @@ def AssembleTestSuites(): smallSuite.addTest(test_dem_3d_parallel_bond_model.TestParallelBondModel("test_ParallelBondModel_1")) smallSuite.addTest(test_dem_3d_smooth_joint_model.TestSmoothJointModel("test_SmoothJointModel_1")) smallSuite.addTest(test_moving_periodic_boundary.TestMovingPeriodicBoundary("test_MovingPeriodicBoundary")) + smallSuite.addTest(test_servo_control.TestServoControl("test_ServoControl")) # Create a test suit with the selected tests plus all small tests nightSuite = suites['nightly'] diff --git a/applications/DEMApplication/tests/test_servo_control.py b/applications/DEMApplication/tests/test_servo_control.py new file mode 100644 index 000000000000..92156c684b40 --- /dev/null +++ b/applications/DEMApplication/tests/test_servo_control.py @@ -0,0 +1,52 @@ +import os +import KratosMultiphysics as Kratos +from KratosMultiphysics import Logger +Logger.GetDefaultOutput().SetSeverity(Logger.Severity.WARNING) +import KratosMultiphysics.KratosUnittest as KratosUnittest +import KratosMultiphysics.DEMApplication.DEM_analysis_stage + +import auxiliary_functions_for_tests + +this_working_dir_backup = os.getcwd() + +def GetFilePath(fileName): + return os.path.join(os.path.dirname(os.path.realpath(__file__)), fileName) + +class ServoConrolTestSolution(KratosMultiphysics.DEMApplication.DEM_analysis_stage.DEMAnalysisStage, KratosUnittest.TestCase): + + @classmethod + def GetMainPath(self): + return os.path.join(os.path.dirname(os.path.realpath(__file__)), "servo_control_tests_files") + + def GetProblemNameWithPath(self): + return os.path.join(self.main_path, self.DEM_parameters["problem_name"].GetString()) + + def FinalizeSolutionStep(self): + super().FinalizeSolutionStep() + tolerance = 1e-8 + stress_tensor = self.MeasureSphereForGettingGlobalStressTensor() + mean_stress = (stress_tensor[0][0]+stress_tensor[1][1]+stress_tensor[2][2])/3 + + if self.time >= 0.000048 and self.time < 0.000049: + expected_value = 5.007421025614752e-08 + self.assertAlmostEqual(mean_stress, expected_value, delta=tolerance) + + def Finalize(self): + self.procedures.RemoveFoldersWithResults(str(self.main_path), str(self.problem_name), '') + super().Finalize() + +class TestServoControl(KratosUnittest.TestCase): + + def setUp(self): + pass + + @classmethod + def test_ServoControl(self): + path = os.path.join(os.path.dirname(os.path.realpath(__file__)), "servo_control_tests_files") + parameters_file_name = os.path.join(path, "ProjectParametersDEM.json") + model = Kratos.Model() + auxiliary_functions_for_tests.CreateAndRunStageInSelectedNumberOfOpenMPThreads(ServoConrolTestSolution, model, parameters_file_name, 1) + +if __name__ == "__main__": + Kratos.Logger.GetDefaultOutput().SetSeverity(Logger.Severity.WARNING) + KratosUnittest.main() From 5289b98f26bf0c974570daa1adb362053dd0e81e Mon Sep 17 00:00:00 2001 From: markelov208 <146726001+markelov208@users.noreply.github.com> Date: Wed, 30 Oct 2024 16:35:34 +0100 Subject: [PATCH 6/6] [GeoMechanicsApplication] Making unit tests for the erosion mechanism (#12789) * first version of unit tests * corrected second unit test * removed unused definition * formatted * formatetd 2.0 * added this-> to call of BaseClassFinalizeSolutionStep * updated README file * changed GravitationalForce * changed to GravitationalAcceleration * response to all reviews * fixed that * response to review on geo_mechanics_newton_raphson_erosion_process_strategy * removed some include, p_scheme and p_solver --- .../steady_state_Pw_piping_element.cpp | 2 - .../custom_strategies/strategies/README.md | 6 + ...ewton_raphson_erosion_process_strategy.hpp | 35 ++- ...ewton_rapson_erosion_processs_strategy.cpp | 216 ++++++++++++++++++ 4 files changed, 235 insertions(+), 24 deletions(-) create mode 100644 applications/GeoMechanicsApplication/tests/cpp_tests/test_geo_mechanics_newton_rapson_erosion_processs_strategy.cpp diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp index eb6f8f7658ca..2d90908bc669 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp @@ -18,7 +18,6 @@ namespace Kratos { - template Element::Pointer SteadyStatePwPipingElement::Create(IndexType NewId, NodesArrayType const& ThisNodes, @@ -341,5 +340,4 @@ bool SteadyStatePwPipingElement::InEquilibrium(const Properties template class SteadyStatePwPipingElement<2, 4>; template class SteadyStatePwPipingElement<3, 6>; template class SteadyStatePwPipingElement<3, 8>; - } // Namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md b/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md index b819a0ebaa4b..52febd5a1e65 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/README.md @@ -54,6 +54,12 @@ This function performs an iterative process towards certain equilibrium values o the modelling is performed. Then a pipe height is corrected by a given increment for each open piping element. The maximum number of iterations is the user input, "max_piping_iterations". +The typical case is the following. Initially, the equilibrium pipe height is small and quickly the pipe height is +increased above the equilibrium height. At this moment PIPE_EROSION is set to true. Then due to the flow development, +the water pressure gradient decreases, which leads to a higher value of the equilibrium height. As a result the +equilibrium pipe height grows faster than the pipe height. When both of the heights become equal, the iterative process +is stopped. + ![check_pipe_equilibrium.svg](check_pipe_equilibrium.svg) ### Recalculate function diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp index aca6fca4a05d..1a3f8f443e97 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp @@ -31,7 +31,6 @@ namespace Kratos { - template class GeoMechanicsNewtonRaphsonErosionProcessStrategy : public GeoMechanicsNewtonRaphsonStrategy @@ -85,7 +84,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy if (PipeElements.empty()) { KRATOS_INFO_IF("PipingLoop", this->GetEchoLevel() > 0 && rank == 0) << "No Pipe Elements -> Finalizing Solution " << std::endl; - GeoMechanicsNewtonRaphsonStrategy::FinalizeSolutionStep(); + this->BaseClassFinalizeSolutionStep(); return; } // calculate max pipe height and pipe increment @@ -124,13 +123,13 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy save_or_reset_pipe_heights(OpenPipeElements, grow); } // recalculate groundwater flow - bool converged = Recalculate(); + bool converged = this->Recalculate(); // error check KRATOS_ERROR_IF_NOT(converged) << "Groundwater flow calculation failed to converge." << std::endl; } - GeoMechanicsNewtonRaphsonStrategy::FinalizeSolutionStep(); + this->BaseClassFinalizeSolutionStep(); } /** @@ -140,7 +139,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy int Check() override { KRATOS_TRY - BaseType::Check(); this->GetBuilderAndSolver()->Check(BaseType::GetModelPart()); this->GetScheme()->Check(BaseType::GetModelPart()); @@ -296,25 +294,22 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy return max_pipe_height / (n_steps - 1); } - bool Recalculate() + virtual bool Recalculate() { KRATOS_INFO_IF("ResidualBasedNewtonRaphsonStrategy", this->GetEchoLevel() > 0 && rank == 0) << "Recalculating" << std::endl; - // KRATOS_INFO_IF("PipingLoop") << "Recalculating" << std::endl; - // ModelPart& CurrentModelPart = this->GetModelPart(); - // this->Clear(); - - // Reset displacements to the initial (Assumes Water Pressure is the convergence criteria) - /* block_for_each(CurrentModelPart.Nodes(), [&](Node& rNode) { - auto dold = rNode.GetSolutionStepValue(WATER_PRESSURE, 1); - rNode.GetSolutionStepValue(WATER_PRESSURE, 0) = dold; - });*/ GeoMechanicsNewtonRaphsonStrategy::InitializeSolutionStep(); GeoMechanicsNewtonRaphsonStrategy::Predict(); return GeoMechanicsNewtonRaphsonStrategy::SolveSolutionStep(); } + virtual void BaseClassFinalizeSolutionStep() + { + // to override in a unit test + GeoMechanicsNewtonRaphsonStrategy::FinalizeSolutionStep(); + } + bool check_pipe_equilibrium(filtered_elements open_pipe_elements, double amax, unsigned int mPipingIterations) { bool equilibrium = false; @@ -331,7 +326,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy equilibrium = true; // perform a flow calculation and stop growing if the calculation doesn't converge - converged = Recalculate(); + converged = this->Recalculate(); // todo: JDN (20220817) : grow not used. // if (!converged) @@ -344,7 +339,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy equilibrium = true; for (auto OpenPipeElement : open_pipe_elements) { auto pElement = static_cast*>(OpenPipeElement); - // get open pipe element geometry and properties auto& Geom = OpenPipeElement->GetGeometry(); auto& prop = OpenPipeElement->GetProperties(); @@ -373,7 +367,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy equilibrium = true; OpenPipeElement->SetValue(PIPE_HEIGHT, small_pipe_height); } - // calculate difference between equilibrium height and current pipe height OpenPipeElement->SetValue(DIFF_PIPE_HEIGHT, eq_height - current_height); } @@ -437,7 +430,5 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy } } } - -}; // Class GeoMechanicsNewtonRaphsonStrategy - -} // namespace Kratos \ No newline at end of file +}; // Class GeoMechanicsNewtonRaphsonErosionProcessStrategy +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_geo_mechanics_newton_rapson_erosion_processs_strategy.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_geo_mechanics_newton_rapson_erosion_processs_strategy.cpp new file mode 100644 index 000000000000..bbe6e65e0697 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_geo_mechanics_newton_rapson_erosion_processs_strategy.cpp @@ -0,0 +1,216 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Gennady Markelov +// + +// Project includes +#include "custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp" +#include "geo_mechanics_fast_suite.h" +#include "test_utilities.h" + +#include +#include +#include + +namespace Kratos +{ +template +class MockGeoMechanicsNewtonRaphsonErosionProcessStrategy + : public GeoMechanicsNewtonRaphsonErosionProcessStrategy +{ +public: + using SolvingStrategyType = ImplicitSolvingStrategy; + using TConvergenceCriteriaType = ConvergenceCriteria; + using TBuilderAndSolverType = typename SolvingStrategyType::TBuilderAndSolverType; + using TSchemeType = typename SolvingStrategyType::TSchemeType; + + MockGeoMechanicsNewtonRaphsonErosionProcessStrategy(ModelPart& rModelPart, + typename TSchemeType::Pointer pScheme, + typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, + typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, + Parameters& rParameters) + : GeoMechanicsNewtonRaphsonErosionProcessStrategy( + rModelPart, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters), + mrModelPart(rModelPart) + { + } + +private: + bool Recalculate() override + { + // The function provides a pressure change over iterations that leads to a non-linear growth of an equilibrium pipe height. + // As a result the height curve crosses twice a pipe height growth curve in the search algorithm implemented + // in GeoMechanicsNewtonRaphsonErosionProcessStrategy. The search stops at the second point. + for (auto& node : mrModelPart.Nodes()) { + node.FastGetSolutionStepValue(WATER_PRESSURE) = + std::pow(node.FastGetSolutionStepValue(WATER_PRESSURE), 0.9725); + } + return true; + } + + void BaseClassFinalizeSolutionStep() override + { + // to avoid calling the BaseClassFinalizeSolutionStep in tested class. + } + + ModelPart& mrModelPart; +}; + +Node::Pointer CreateNodeWithSolutionStepValues(ModelPart& rModelPart, int Id, double CoordinateX, double CoordinateY, double WaterPressure) +{ + constexpr double coordinate_z = 0.0; + auto p_node = rModelPart.CreateNewNode(Id, CoordinateX, CoordinateY, coordinate_z); + + p_node->FastGetSolutionStepValue(WATER_PRESSURE) = WaterPressure; + const array_1d GravitationalAcceleration{0, -9.81, 0}; + p_node->FastGetSolutionStepValue(VOLUME_ACCELERATION) = GravitationalAcceleration; + + return p_node; +} + +Geometry::Pointer CreateQuadrilateral2D4N(ModelPart& rModelPart, + const std::vector& rNodeIds, + double Xmin, + double Xmax, + double WaterPressureLeft, + double WaterPressureRight) +{ + constexpr double y_min = -0.1; + constexpr double y_max = 0.1; + + return make_shared>( + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[0], Xmin, y_min, WaterPressureLeft), + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[1], Xmax, y_min, WaterPressureRight), + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[2], Xmax, y_max, WaterPressureRight), + CreateNodeWithSolutionStepValues(rModelPart, rNodeIds[3], Xmin, y_max, WaterPressureLeft)); +} + +auto SetupPipingStrategy(Model& rModel) +{ + using SparseSpaceType = UblasSpace; + using LocalSpaceType = UblasSpace; + using LinearSolverType = SkylineLUFactorizationSolver; + using ConvergenceCriteriaType = ConvergenceCriteria; + using MixedGenericCriteriaType = MixedGenericCriteria; + using ConvergenceVariableListType = MixedGenericCriteriaType::ConvergenceVariableListType; + + auto& r_model_part = rModel.CreateModelPart("ModelPart", 1); + r_model_part.AddNodalSolutionStepVariable(WATER_PRESSURE); + r_model_part.AddNodalSolutionStepVariable(VOLUME_ACCELERATION); + const auto& r_process_info = r_model_part.GetProcessInfo(); + + auto p_element_props = r_model_part.CreateNewProperties(0); + p_element_props->SetValue(CONSTITUTIVE_LAW, LinearElastic2DInterfaceLaw().Clone()); + + auto p_element = make_intrusive>( + 0, CreateQuadrilateral2D4N(r_model_part, std::vector{13, 14, 15, 16}, 3.0, 4.0, 2000.0, 2000.0), + p_element_props, std::make_unique()); + p_element->Initialize(r_process_info); + r_model_part.AddElement(p_element); + + // Set the pipe element properties + auto p_piping_element_props = r_model_part.CreateNewProperties(1); + p_piping_element_props->SetValue(PIPE_D_70, 2e-4); + p_piping_element_props->SetValue(PIPE_ETA, 0.25); + p_piping_element_props->SetValue(PIPE_THETA, 30); + p_piping_element_props->SetValue(DENSITY_SOLID, 2650); + p_piping_element_props->SetValue(DENSITY_WATER, 1000); + p_piping_element_props->SetValue(PIPE_ELEMENT_LENGTH, 1.0); + p_piping_element_props->SetValue(PIPE_MODIFIED_D, false); + p_piping_element_props->SetValue(PIPE_MODEL_FACTOR, 1); + p_piping_element_props->SetValue(PIPE_START_ELEMENT, 1); + p_piping_element_props->SetValue(CONSTITUTIVE_LAW, LinearElastic2DInterfaceLaw().Clone()); + + // Create the start piping element + auto p_piping_element = make_intrusive>( + 1, CreateQuadrilateral2D4N(r_model_part, std::vector{1, 2, 3, 4}, 0.0, 1.0, 0.0, 500.0), + p_piping_element_props, std::make_unique()); + p_piping_element->Initialize(r_process_info); + r_model_part.AddElement(p_piping_element); + + // Create other piping elements + p_piping_element = make_intrusive>( + 2, CreateQuadrilateral2D4N(r_model_part, std::vector{5, 6, 7, 8}, 2.0, 3.0, 500.0, 1000.0), + p_piping_element_props, std::make_unique()); + p_piping_element->Initialize(r_process_info); + r_model_part.AddElement(p_piping_element); + + p_piping_element = make_intrusive>( + 3, CreateQuadrilateral2D4N(r_model_part, std::vector{9, 10, 11, 12}, 1.0, 2.0, 1000.0, 1500.0), + p_piping_element_props, std::make_unique()); + p_piping_element->Initialize(r_process_info); + r_model_part.AddElement(p_piping_element); + + Parameters p_parameters(R"( + { + "max_piping_iterations": 25 + } )"); + + const auto p_builder_and_solver = + Kratos::make_shared>(nullptr); + const ConvergenceVariableListType convergence_settings{}; + const ConvergenceCriteriaType::Pointer p_criteria = + std::make_shared(convergence_settings); + + return std::make_unique>( + r_model_part, nullptr, p_criteria, p_builder_and_solver, p_parameters); +} +} // namespace Kratos + +namespace Kratos::Testing +{ +KRATOS_TEST_CASE_IN_SUITE(CheckGetPipingElements, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + // Arrange + Model model; + auto p_solving_strategy = SetupPipingStrategy(model); + + // Act + const auto piping_elements = p_solving_strategy->GetPipingElements(); + + // Assert + const auto elements = model.GetModelPart("ModelPart").Elements(); + constexpr int number_of_elements = 4; + KRATOS_EXPECT_EQ(elements.size(), number_of_elements); + + constexpr int number_of_piping_elements = 3; + KRATOS_EXPECT_EQ(piping_elements.size(), number_of_piping_elements); + KRATOS_EXPECT_EQ(piping_elements[0]->GetId(), 1); + KRATOS_EXPECT_EQ(piping_elements[1]->GetId(), 3); + KRATOS_EXPECT_EQ(piping_elements[2]->GetId(), 2); +} + +KRATOS_TEST_CASE_IN_SUITE(CheckFinalizeSolutionStep, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + // Arrange + Model model; + auto p_solving_strategy = SetupPipingStrategy(model); + + // Act + p_solving_strategy->FinalizeSolutionStep(); + + // Assert + auto elements = model.GetModelPart("ModelPart").Elements(); + + constexpr int active = 1; + KRATOS_EXPECT_EQ(elements[1].GetValue(PIPE_ACTIVE), active); + + constexpr double expected_active_pipe_height = 0.0183333; + KRATOS_EXPECT_NEAR(elements[1].GetValue(PIPE_HEIGHT), expected_active_pipe_height, Defaults::relative_tolerance); + + constexpr int inactive = 0; + KRATOS_EXPECT_EQ(elements[2].GetValue(PIPE_ACTIVE), inactive); + KRATOS_EXPECT_EQ(elements[3].GetValue(PIPE_ACTIVE), inactive); + + constexpr double expected_inactive_pipe_height = 0.0; + KRATOS_EXPECT_NEAR(elements[2].GetValue(PIPE_HEIGHT), expected_inactive_pipe_height, Defaults::relative_tolerance); + KRATOS_EXPECT_NEAR(elements[3].GetValue(PIPE_HEIGHT), expected_inactive_pipe_height, Defaults::relative_tolerance); +} +} // namespace Kratos::Testing \ No newline at end of file