diff --git a/applications/DEMApplication/DEM_application.cpp b/applications/DEMApplication/DEM_application.cpp index 6cb1ef5dafba..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) @@ -153,6 +154,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) @@ -612,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) @@ -660,6 +663,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..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) @@ -125,6 +126,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..ce7366a02521 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[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); + } } if (r_process_info[ENERGY_CALCULATION_OPTION]){ 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 94f647c3c0e6..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() @@ -530,32 +544,98 @@ 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: + 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.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): + 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): + + # 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_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 + + 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 UpdateSearchStartegyAndCPlusPlusStrategy(self): + 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) + 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): delta_time = self.spheres_model_part.ProcessInfo.GetValue(DELTA_TIME) - move_velocity = self.DEM_parameters["BoundingBoxMoveVelocity"].GetDouble() + # 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, @@ -571,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 @@ -742,7 +829,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 +868,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 @@ -1107,6 +1291,47 @@ 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 + + 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 @@ -1643,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 52476b1f7d73..e35bae680caa 100644 --- a/applications/DEMApplication/python_scripts/dem_default_input_parameters.py +++ b/applications/DEMApplication/python_scripts/dem_default_input_parameters.py @@ -20,7 +20,17 @@ 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, + "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 d9ccadeebd4e..3d8599bcbfe5 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.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) + # PRINTING VARIABLES self.spheres_model_part.ProcessInfo.SetValue(PRINT_EXPORT_ID, self.print_export_id) 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, 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()