diff --git a/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py b/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py new file mode 100644 index 000000000000..d1beb7e060ee --- /dev/null +++ b/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py @@ -0,0 +1,265 @@ +import KratosMultiphysics as KM +import numpy as np +import time + + +def Factory(settings, model): + if not isinstance(settings, KM.Parameters): + raise Exception("expected input shall be a Parameters object, encapsulating a json string") + return RotatingFrameProcess(model, settings["Parameters"]) + + +class RotatingFrameProcess(KM.Process): + """Process to apply rigid-body rotational kinematics to model parts. + + This process: + - rotates a model part and assigns the corresponding MESH_DISPLACEMENT + - assigns the corresponding rigid-body rotational VELOCITY to another model part + + The angular velocity is provided as the standard scalar input + 'angular_velocity_radians = f(x,y,z,t)' (or a constant number). + The rotation is advanced by integrating the angular-velocity increment + at each solution step. + """ + + def __init__(self, model, settings): + KM.Process.__init__(self) + + default_settings = KM.Parameters("""{ + "rotating_frame_model_part_name": "", + "rotating_object_model_part_name": "", + "interval": [0.0, 1e30], + "center_of_rotation": [0.0, 0.0, 0.0], + "axis_of_rotation": [1.0, 0.0, 0.0], + "angular_velocity_radians": 0.0, + "echo_level": 0, + "fix_mesh_displacement": false, + "fix_velocity": false + }""") + + # Allow a string function expression for 'angular_velocity_radians'. + if settings.Has("angular_velocity_radians") and settings["angular_velocity_radians"].IsString(): + default_settings["angular_velocity_radians"].SetString("0.0") + + # Assign this here since it may normalize the "interval" before validation + self.interval = KM.IntervalUtility(settings) + + settings.ValidateAndAssignDefaults(default_settings) + + self.model = model + self.settings = settings + + if not settings["rotating_frame_model_part_name"].GetString(): + raise Exception("'rotating_frame_model_part_name' not provided. Please specify the model part to be rotated.") + self.rotating_frame_model_part_name = settings["rotating_frame_model_part_name"].GetString() + + if not settings["rotating_object_model_part_name"].GetString(): + raise Exception("'rotating_object_model_part_name' not provided. Please specify the model part to which rotational velocity will be assigned.") + self.rotating_object_model_part_name = settings["rotating_object_model_part_name"].GetString() + + self.rotating_frame_model_part = self.model.GetModelPart(self.rotating_frame_model_part_name) + self.rotating_object_model_part = self.model.GetModelPart(self.rotating_object_model_part_name) + self._interval_begin = settings["interval"][0].GetDouble() + self._interval_end = settings["interval"][1].GetDouble() + + self.center_of_rotation = np.array(settings["center_of_rotation"].GetVector(), dtype=float) + + self.axis_of_rotation = np.array(settings["axis_of_rotation"].GetVector(), dtype=float) + if self.axis_of_rotation.size == 0: + raise Exception("The 'axis_of_rotation' vector is empty.") + + axis_norm = np.linalg.norm(self.axis_of_rotation) + if np.isclose(axis_norm, 0.0): + raise Exception("The 'axis_of_rotation' vector must have non-zero norm.") + + if not np.isclose(axis_norm, 1.0, rtol=1e-6): + KM.Logger.PrintWarning( + "RotatingFrameProcess", + "The 'axis_of_rotation' vector is not a unit vector. It will be normalized." + ) + self.axis_of_rotation /= axis_norm + + self.angular_velocity_function = self._CreateScalarFunction( + settings["angular_velocity_radians"], + "angular_velocity_radians" + ) + + self.echo_level = settings["echo_level"].GetInt() + if self.echo_level < 0: + raise Exception("The 'echo_level' parameter must be >= 0.") + + self.fix_mesh_displacement = settings["fix_mesh_displacement"].GetBool() + self.fix_velocity = settings["fix_velocity"].GetBool() + self._fixity_applied = False + self._python_cache_initialized = False + + def ExecuteInitializeSolutionStep(self): + current_time = self.rotating_frame_model_part.ProcessInfo[KM.TIME] + delta_theta, omega, is_active = self._GetDeltaThetaAndOmegaFromFunction(current_time) + + if np.isclose(delta_theta, 0.0) and not is_active: + # If we leave the active interval, release any fixity that this process applied. + if self._fixity_applied: + if self.fix_mesh_displacement: + self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT, False) + if self.fix_velocity: + self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY, False) + self._fixity_applied = False + return + + if not self._python_cache_initialized: + self._InitializePythonBackendCaches() + + t0 = time.perf_counter() + if not np.isclose(delta_theta, 0.0): + self._ApplyRotationAndMeshDisplacementPython(delta_theta) + rot_time = time.perf_counter() - t0 + + t0 = time.perf_counter() + if is_active: + self._AssignRotationalVelocityPython(omega) + vel_time = time.perf_counter() - t0 + + if self.echo_level > 0: + KM.Logger.PrintInfo( + "RotatingFrameProcess", f" rotation={rot_time*1.0e3:.3f} ms, velocity={vel_time*1.0e3:.3f} ms" + ) + + if is_active: + if self.fix_mesh_displacement: + self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT, True) + + if self.fix_velocity: + self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY, True) + + self._fixity_applied = self.fix_mesh_displacement or self.fix_velocity + elif self._fixity_applied: + if self.fix_mesh_displacement: + self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT, False) + if self.fix_velocity: + self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY, False) + self._fixity_applied = False + + @staticmethod + def _CreateScalarFunction(value, value_name): + if value.IsNumber(): + function_string = str(value.GetDouble()) + elif value.IsString(): + function_string = value.GetString() + if function_string == "": + raise Exception(f"The '{value_name}' function string is empty.") + else: + raise Exception( + f"The '{value_name}' parameter must be provided as a number or as a function string." + ) + return KM.GenericFunctionUtility(function_string) + + def _EvaluateAngularVelocityFunction(self, time): + x = self.center_of_rotation[0] + y = self.center_of_rotation[1] + z = self.center_of_rotation[2] + return self.angular_velocity_function.CallFunction(x, y, z, time, 0.0, 0.0, 0.0) + + def _GetDeltaThetaAndOmegaFromFunction(self, time): + omega = self._EvaluateAngularVelocityFunction(time) + delta_time = self.rotating_frame_model_part.ProcessInfo[KM.DELTA_TIME] + delta_theta = 0.0 + if delta_time > 0.0: + step_begin = time - delta_time + overlap_begin = max(step_begin, self._interval_begin) + overlap_end = min(time, self._interval_end) + + if overlap_end > overlap_begin: + omega_begin = self._EvaluateAngularVelocityFunction(overlap_begin) + omega_end = self._EvaluateAngularVelocityFunction(overlap_end) + delta_theta = 0.5 * (omega_begin + omega_end) * (overlap_end - overlap_begin) + + return delta_theta, omega, self.interval.IsInInterval(time) + + def _InitializePythonBackendCaches(self): + self._rf_initial_pos_ta = KM.TensorAdaptors.NodePositionTensorAdaptor( + self.rotating_frame_model_part.Nodes, KM.Configuration.Initial + ) + self._rf_current_pos_ta = KM.TensorAdaptors.NodePositionTensorAdaptor( + self.rotating_frame_model_part.Nodes, KM.Configuration.Current + ) + self._rf_mesh_disp_ta = KM.TensorAdaptors.HistoricalVariableTensorAdaptor( + self.rotating_frame_model_part.Nodes, KM.MESH_DISPLACEMENT + ) + self._ro_current_pos_ta = KM.TensorAdaptors.NodePositionTensorAdaptor( + self.rotating_object_model_part.Nodes, KM.Configuration.Current + ) + self._ro_velocity_ta = KM.TensorAdaptors.HistoricalVariableTensorAdaptor( + self.rotating_object_model_part.Nodes, KM.VELOCITY + ) + + # Initial frame coordinates are constant for a static mesh. + self._rf_initial_pos_ta.CollectData() + self._rf_initial_positions = self._rf_initial_pos_ta.data.copy() + + self._rf_rot_matrix = np.empty((3, 3), dtype=float) + self._ro_relative_positions = np.empty( + (self.rotating_object_model_part.NumberOfNodes(), 3), + dtype=float + ) + + self._python_cache_initialized = True + + @staticmethod + def _FixVectorVariable(model_part, variable, is_fixed): + components_map = { + KM.VELOCITY: (KM.VELOCITY_X, KM.VELOCITY_Y, KM.VELOCITY_Z), + KM.MESH_DISPLACEMENT: (KM.MESH_DISPLACEMENT_X, KM.MESH_DISPLACEMENT_Y, KM.MESH_DISPLACEMENT_Z), + } + + if variable not in components_map: + raise Exception(f"Fixing variable '{variable.Name()}' is not supported by this process.") + + vu = KM.VariableUtils() + for comp in components_map[variable]: + vu.ApplyFixity(comp, is_fixed, model_part.Nodes) + + def _ApplyRotationAndMeshDisplacementPython(self, delta_theta): + sin_half_theta = np.sin(delta_theta / 2.0) + a = np.cos(delta_theta / 2.0) + b = -self.axis_of_rotation[0] * sin_half_theta + c = -self.axis_of_rotation[1] * sin_half_theta + d = -self.axis_of_rotation[2] * sin_half_theta + + rot_matrix = self._rf_rot_matrix + rot_matrix[0, 0] = a * a + b * b - c * c - d * d + rot_matrix[0, 1] = 2.0 * (b * c - a * d) + rot_matrix[0, 2] = 2.0 * (b * d + a * c) + rot_matrix[1, 0] = 2.0 * (b * c + a * d) + rot_matrix[1, 1] = a * a + c * c - b * b - d * d + rot_matrix[1, 2] = 2.0 * (c * d - a * b) + rot_matrix[2, 0] = 2.0 * (b * d - a * c) + rot_matrix[2, 1] = 2.0 * (c * d + a * b) + rot_matrix[2, 2] = a * a + d * d - b * b - c * c + + # Rotate current coordinates around the center with the incremental angle. + self._rf_current_pos_ta.CollectData() + current_positions = self._rf_current_pos_ta.data + current_positions[:] = (current_positions - self.center_of_rotation) @ rot_matrix + current_positions += self.center_of_rotation + + # Mesh displacement from reference configuration: u_mesh = x - x0 + self._rf_mesh_disp_ta.data[:] = current_positions - self._rf_initial_positions + + self._rf_mesh_disp_ta.StoreData() + self._rf_current_pos_ta.StoreData() + + def _AssignRotationalVelocityPython(self, omega): + angular_velocity_vector = omega * self.axis_of_rotation + + self._ro_current_pos_ta.CollectData() + self._ro_relative_positions[:] = self._ro_current_pos_ta.data - self.center_of_rotation + + # v = omega x r, matching C++ MathUtils::CrossProduct(omega, r) + v = self._ro_velocity_ta.data + r = self._ro_relative_positions + v[:, 0] = angular_velocity_vector[1] * r[:, 2] - angular_velocity_vector[2] * r[:, 1] + v[:, 1] = angular_velocity_vector[2] * r[:, 0] - angular_velocity_vector[0] * r[:, 2] + v[:, 2] = angular_velocity_vector[0] * r[:, 1] - angular_velocity_vector[1] * r[:, 0] + + self._ro_velocity_ta.StoreData() diff --git a/applications/MeshMovingApplication/tests/test_mdpa_files/rotating_frame_process_test.mdpa b/applications/MeshMovingApplication/tests/test_mdpa_files/rotating_frame_process_test.mdpa new file mode 100644 index 000000000000..1b10409a66c8 --- /dev/null +++ b/applications/MeshMovingApplication/tests/test_mdpa_files/rotating_frame_process_test.mdpa @@ -0,0 +1,214 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes + 1 -0.0000000023 -2.0000000000 0.0000000000 + 2 0.2046309403 -1.5665745249 0.0000000000 + 3 0.5176380872 -1.9318516534 0.0000000000 + 4 -0.5176380916 -1.9318516522 0.0000000000 + 5 -0.4626014305 -1.4460953968 0.0000000000 + 6 0.0000000000 -1.0000000000 0.0000000000 + 7 0.9999999970 -1.7320508093 0.0000000000 + 8 -1.0000000010 -1.7320508070 0.0000000000 + 9 -1.0000000000 -1.0000000000 0.0000000000 + 10 1.0000000000 -1.0000000000 0.0000000000 + 11 -1.4142135628 -1.4142135620 0.0000000000 + 12 1.4142135596 -1.4142135651 0.0000000000 + 13 -1.7320508076 -1.0000000000 0.0000000000 + 14 1.7320508076 -1.0000000000 0.0000000000 + 15 -1.0000000000 0.0000000000 0.0000000000 + 16 1.0000000000 0.0000000000 0.0000000000 + 17 1.5282662467 -0.3654381624 0.0000000000 + 18 -1.5320281206 -0.3586999861 0.0000000000 + 19 -1.9318516536 -0.5176380864 0.0000000000 + 20 1.9318516524 -0.5176380908 0.0000000000 + 21 1.5320281348 0.3586999917 0.0000000000 + 22 -1.5282662625 0.3654381666 0.0000000000 + 23 2.0000000000 -0.0000000012 0.0000000000 + 24 -2.0000000000 0.0000000035 0.0000000000 + 25 0.0000000000 1.0000000000 0.0000000000 + 26 -1.0000000000 1.0000000000 0.0000000000 + 27 1.0000000000 1.0000000000 0.0000000000 + 28 1.9318516530 0.5176380888 0.0000000000 + 29 -1.9318516518 0.5176380932 0.0000000000 + 30 1.7320508087 0.9999999980 0.0000000000 + 31 -1.7320508064 1.0000000020 0.0000000000 + 32 0.4626014313 1.4460953964 0.0000000000 + 33 -0.2046309407 1.5665745245 0.0000000000 + 34 1.4142135646 1.4142135601 0.0000000000 + 35 -1.4142135613 1.4142135634 0.0000000000 + 36 1.0000000030 1.7320508058 0.0000000000 + 37 -0.9999999990 1.7320508081 0.0000000000 + 38 -0.5176380896 1.9318516527 0.0000000000 + 39 0.5176380940 1.9318516516 0.0000000000 + 40 0.0000000000 2.0000000000 0.0000000000 +End Nodes + + +Begin Elements Element2D3N// GUI group identifier: RotatingFrame + 1 0 40 38 33 + 2 0 13 11 9 + 3 0 14 20 17 + 4 0 31 29 22 + 5 0 1 3 2 + 6 0 30 34 27 + 7 0 37 35 26 + 8 0 8 4 5 + 9 0 23 28 21 + 10 0 24 19 18 + 11 0 7 12 10 + 12 0 36 39 32 + 13 0 38 37 33 + 14 0 11 8 9 + 15 0 20 23 17 + 16 0 29 24 22 + 17 0 3 7 2 + 18 0 34 36 27 + 19 0 35 31 26 + 20 0 4 1 2 + 21 0 28 30 21 + 22 0 19 13 18 + 23 0 12 14 10 + 24 0 39 40 33 + 25 0 5 4 2 + 26 0 32 39 33 + 27 0 31 22 26 + 28 0 8 5 9 + 29 0 24 18 22 + 30 0 14 17 10 + 31 0 23 21 17 + 32 0 36 32 27 + 33 0 26 22 15 + 34 0 9 5 6 + 35 0 10 17 16 + 36 0 27 32 25 + 37 0 15 22 18 + 38 0 15 18 9 + 39 0 6 5 2 + 40 0 6 2 10 + 41 0 16 17 21 + 42 0 16 21 27 + 43 0 25 32 33 + 44 0 25 33 26 + 45 0 13 9 18 + 46 0 37 26 33 + 47 0 30 27 21 + 48 0 7 10 2 +End Elements + +Begin SubModelPart GENERIC_RotatingObject // Group RotatingObject // Subtree GENERIC + Begin SubModelPartNodes + 6 + 9 + 10 + 15 + 16 + 25 + 26 + 27 + End SubModelPartNodes + Begin SubModelPartElements + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart +Begin SubModelPart GENERIC_RotatingFrame // Group RotatingFrame // Subtree GENERIC + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart \ No newline at end of file diff --git a/applications/MeshMovingApplication/tests/test_results/rotating_frame_process_test_results.json b/applications/MeshMovingApplication/tests/test_results/rotating_frame_process_test_results.json new file mode 100644 index 000000000000..c6b77fbf3b1a --- /dev/null +++ b/applications/MeshMovingApplication/tests/test_results/rotating_frame_process_test_results.json @@ -0,0 +1,805 @@ +{ + "TIME": [ + 1.0 + ], + "NODE_1": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 1.088042226008604 + ], + "MESH_DISPLACEMENT_Y": [ + 3.678143056901656 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_2": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 0.4759186770532115 + ], + "MESH_DISPLACEMENT_Y": [ + 2.992366158384402 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_3": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 0.09899461402101972 + ], + "MESH_DISPLACEMENT_Y": [ + 3.8344194214044034 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_4": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 2.002941558530337 + ], + "MESH_DISPLACEMENT_Y": [ + 3.27120732232944 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_5": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 1.6374635443617294 + ], + "MESH_DISPLACEMENT_Y": [ + 2.407807928463774 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_6": { + "VELOCITY_X": [ + 8.390715290764524 + ], + "VELOCITY_Y": [ + -5.440211108893697 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 0.5440211108893698 + ], + "MESH_DISPLACEMENT_Y": [ + 1.8390715290764525 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_7": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -0.8967993181670197 + ], + "MESH_DISPLACEMENT_Y": [ + 3.729386439554764 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_8": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 2.7813437350564936 + ], + "MESH_DISPLACEMENT_Y": [ + 2.6413442146342025 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_9": { + "VELOCITY_X": [ + 2.9505041818708264 + ], + "VELOCITY_Y": [ + -13.830926399658221 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 2.383092639965822 + ], + "MESH_DISPLACEMENT_Y": [ + 1.2950504181870826 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_10": { + "VELOCITY_X": [ + 13.830926399658221 + ], + "VELOCITY_Y": [ + 2.9505041818708264 + ], + "VELOCITY_Z": [ + -0.0 + ], + "MESH_DISPLACEMENT_X": [ + -1.2950504181870826 + ], + "MESH_DISPLACEMENT_Y": [ + 2.383092639965822 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_11": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 3.3702019324133063 + ], + "MESH_DISPLACEMENT_Y": [ + 1.8314778644387268 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_12": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -1.8314778587737066 + ], + "MESH_DISPLACEMENT_Y": [ + 3.37020193533752 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_13": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 3.7293864380604065 + ], + "MESH_DISPLACEMENT_Y": [ + 0.8967993246090703 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_14": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -2.6413442162816665 + ], + "MESH_DISPLACEMENT_Y": [ + 2.7813437335438347 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_15": { + "VELOCITY_X": [ + -5.440211108893697 + ], + "VELOCITY_Y": [ + -8.390715290764524 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 1.8390715290764525 + ], + "MESH_DISPLACEMENT_Y": [ + -0.5440211108893698 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_16": { + "VELOCITY_X": [ + 5.440211108893697 + ], + "VELOCITY_Y": [ + 8.390715290764524 + ], + "VELOCITY_Z": [ + -0.0 + ], + "MESH_DISPLACEMENT_X": [ + -1.8390715290764525 + ], + "MESH_DISPLACEMENT_Y": [ + 0.5440211108893698 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_17": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -2.611784868084282 + ], + "MESH_DISPLACEMENT_Y": [ + 1.5034760213723186 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_18": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 3.012649663254089 + ], + "MESH_DISPLACEMENT_Y": [ + -0.17378070816593616 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_19": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 3.8344194213370004 + ], + "MESH_DISPLACEMENT_Y": [ + -0.09899461560108136 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_20": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -3.2712073231324705 + ], + "MESH_DISPLACEMENT_Y": [ + 2.002941557167884 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_21": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -3.012649692415423 + ], + "MESH_DISPLACEMENT_Y": [ + 0.17378070559223535 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_22": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 2.611784894856724 + ], + "MESH_DISPLACEMENT_Y": [ + -1.5034760376919525 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_23": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -3.6781430575000797 + ], + "MESH_DISPLACEMENT_Y": [ + 1.0880422239856256 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_24": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 3.678143056248831 + ], + "MESH_DISPLACEMENT_Y": [ + -1.08804222821549 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_25": { + "VELOCITY_X": [ + -8.390715290764524 + ], + "VELOCITY_Y": [ + 5.440211108893697 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -0.5440211108893698 + ], + "MESH_DISPLACEMENT_Y": [ + -1.8390715290764525 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_26": { + "VELOCITY_X": [ + -13.830926399658221 + ], + "VELOCITY_Y": [ + -2.9505041818708264 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 1.2950504181870826 + ], + "MESH_DISPLACEMENT_Y": [ + -2.383092639965822 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_27": { + "VELOCITY_X": [ + -2.9505041818708264 + ], + "VELOCITY_Y": [ + 13.830926399658221 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -2.383092639965822 + ], + "MESH_DISPLACEMENT_Y": [ + -1.2950504181870826 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_28": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -3.8344194215392085 + ], + "MESH_DISPLACEMENT_Y": [ + 0.09899461086089689 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_29": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 3.271207320723377 + ], + "MESH_DISPLACEMENT_Y": [ + -2.002941561255243 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_30": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -3.7293864389953426 + ], + "MESH_DISPLACEMENT_Y": [ + -0.8967993203325041 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_31": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 2.6413442129867386 + ], + "MESH_DISPLACEMENT_Y": [ + -2.781343736569152 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_32": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -1.637463545615378 + ], + "MESH_DISPLACEMENT_Y": [ + -2.4078079272929283 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_33": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -0.4759186760999744 + ], + "MESH_DISPLACEMENT_Y": [ + -2.9923661578663823 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_34": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -3.3702019346899945 + ], + "MESH_DISPLACEMENT_Y": [ + -1.8314778599652533 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_35": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 1.8314778628249642 + ], + "MESH_DISPLACEMENT_Y": [ + -3.3702019331359345 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_36": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -2.7813437380818113 + ], + "MESH_DISPLACEMENT_Y": [ + -2.641344211339274 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_37": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + 0.8967993224979884 + ], + "MESH_DISPLACEMENT_Y": [ + -3.7293864384359203 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_38": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -0.09899460922643322 + ], + "MESH_DISPLACEMENT_Y": [ + -3.834419421422704 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_39": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -2.0029415626176963 + ], + "MESH_DISPLACEMENT_Y": [ + -3.2712073199203457 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_40": { + "VELOCITY_X": [ + 0.0 + ], + "VELOCITY_Y": [ + 0.0 + ], + "VELOCITY_Z": [ + 0.0 + ], + "MESH_DISPLACEMENT_X": [ + -1.0880422217787395 + ], + "MESH_DISPLACEMENT_Y": [ + -3.678143058152905 + ], + "MESH_DISPLACEMENT_Z": [ + 0.0 + ] + } +} \ No newline at end of file diff --git a/applications/MeshMovingApplication/tests/test_rotating_frame_process.py b/applications/MeshMovingApplication/tests/test_rotating_frame_process.py new file mode 100644 index 000000000000..852da88837d5 --- /dev/null +++ b/applications/MeshMovingApplication/tests/test_rotating_frame_process.py @@ -0,0 +1,213 @@ +# Kratos Imports +import os + +import KratosMultiphysics +import KratosMultiphysics.KratosUnittest as KratosUnittest +import KratosMultiphysics.kratos_utilities as kratos_utilities +from KratosMultiphysics.json_output_process import JsonOutputProcess +from KratosMultiphysics.from_json_check_result_process import FromJsonCheckResultProcess + + +mesh_moving_is_available = kratos_utilities.CheckIfApplicationsAvailable("MeshMovingApplication") + + +def GetFilePath(file_name): + return os.path.join(os.path.dirname(os.path.realpath(__file__)), file_name) + + +def GenerateModel(): + model = KratosMultiphysics.Model() + model_part = model.CreateModelPart("Main") + + model_part.AddNodalSolutionStepVariable(KratosMultiphysics.MESH_DISPLACEMENT) + model_part.AddNodalSolutionStepVariable(KratosMultiphysics.VELOCITY) + + model_part_io = KratosMultiphysics.ModelPartIO( + GetFilePath("test_mdpa_files/rotating_frame_process_test") + ) + model_part_io.ReadModelPart(model_part) + + KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.MESH_DISPLACEMENT_X, model_part) + KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.MESH_DISPLACEMENT_Y, model_part) + KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.MESH_DISPLACEMENT_Z, model_part) + KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.VELOCITY_X, model_part) + KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.VELOCITY_Y, model_part) + KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.VELOCITY_Z, model_part) + + model_part.ProcessInfo[KratosMultiphysics.TIME] = 1.0 + model_part.ProcessInfo[KratosMultiphysics.DELTA_TIME] = 1.0 + + return model, model_part + + +@KratosUnittest.skipUnless(mesh_moving_is_available, "MeshMovingApplication is not available") +class TestRotatingFrameProcess(KratosUnittest.TestCase): + + @staticmethod + def _GetBaseParameters(): + return KratosMultiphysics.Parameters(""" + { + "rotating_frame_model_part_name": "Main.GENERIC_RotatingFrame", + "rotating_object_model_part_name": "Main.GENERIC_RotatingObject", + "center_of_rotation": [0.0, 0.0, 0.0], + "axis_of_rotation": [0.0, 0.0, 1.0], + "fix_mesh_displacement": true, + "fix_velocity": true + }""") + + def test_rotating_frame_process(self): + model, model_part = GenerateModel() + + self.print_reference_values = False + + parameters = self._GetBaseParameters() + parameters.AddEmptyValue("angular_velocity_radians").SetDouble(-10.0) + + from KratosMultiphysics.MeshMovingApplication.rotating_frame_process import RotatingFrameProcess + + process = RotatingFrameProcess(model, parameters) + process.ExecuteInitializeSolutionStep() + + if self.print_reference_values: + json_output_settings = KratosMultiphysics.Parameters(r'''{ + "output_variables": ["VELOCITY", "MESH_DISPLACEMENT"], + "output_file_name": "", + "model_part_name": "Main", + "time_frequency": 0.0 + }''') + json_output_settings["output_file_name"].SetString( + GetFilePath("test_results/rotating_frame_process_test_results.json") + ) + + json_output_process = JsonOutputProcess(model, json_output_settings) + json_output_process.ExecuteInitialize() + json_output_process.ExecuteBeforeSolutionLoop() + json_output_process.ExecuteFinalizeSolutionStep() + else: + json_check_parameters = KratosMultiphysics.Parameters(r'''{ + "check_variables": ["VELOCITY", "MESH_DISPLACEMENT"], + "input_file_name": "", + "model_part_name": "Main", + "time_frequency": 0.0, + "tolerance": 1e-6 + }''') + json_check_parameters["input_file_name"].SetString( + GetFilePath("test_results/rotating_frame_process_test_results.json") + ) + + json_check_process = FromJsonCheckResultProcess(model, json_check_parameters) + json_check_process.ExecuteInitialize() + json_check_process.ExecuteBeforeSolutionLoop() + json_check_process.ExecuteFinalizeSolutionStep() + + def test_rotating_frame_process_with_angular_velocity_string(self): + numeric_model, numeric_model_part = GenerateModel() + function_model, function_model_part = GenerateModel() + + from KratosMultiphysics.MeshMovingApplication.rotating_frame_process import RotatingFrameProcess + + numeric_parameters = self._GetBaseParameters() + numeric_parameters.AddEmptyValue("angular_velocity_radians").SetDouble(-10.0) + + function_parameters = self._GetBaseParameters() + function_parameters.AddEmptyValue("angular_velocity_radians").SetString("-10.0") + + numeric_process = RotatingFrameProcess(numeric_model, numeric_parameters) + function_process = RotatingFrameProcess(function_model, function_parameters) + + numeric_process.ExecuteInitializeSolutionStep() + function_process.ExecuteInitializeSolutionStep() + + for numeric_node in numeric_model_part.Nodes: + function_node = function_model_part.GetNode(numeric_node.Id) + self.assertVectorAlmostEqual( + numeric_node.GetSolutionStepValue(KratosMultiphysics.MESH_DISPLACEMENT), + function_node.GetSolutionStepValue(KratosMultiphysics.MESH_DISPLACEMENT), + 8 + ) + self.assertVectorAlmostEqual( + numeric_node.GetSolutionStepValue(KratosMultiphysics.VELOCITY), + function_node.GetSolutionStepValue(KratosMultiphysics.VELOCITY), + 8 + ) + + def test_rotating_frame_process_outside_interval_is_no_op(self): + model, model_part = GenerateModel() + model_part.ProcessInfo[KratosMultiphysics.TIME] = 1.0 + + from KratosMultiphysics.MeshMovingApplication.rotating_frame_process import RotatingFrameProcess + + parameters = self._GetBaseParameters() + parameters.AddValue("interval", KratosMultiphysics.Parameters('[2.0, "End"]')) + parameters.AddEmptyValue("angular_velocity_radians").SetDouble(-10.0) + + process = RotatingFrameProcess(model, parameters) + process.ExecuteInitializeSolutionStep() + + for node in model_part.Nodes: + self.assertVectorAlmostEqual( + node.GetSolutionStepValue(KratosMultiphysics.MESH_DISPLACEMENT), + [0.0, 0.0, 0.0], + 12 + ) + self.assertVectorAlmostEqual( + node.GetSolutionStepValue(KratosMultiphysics.VELOCITY), + [0.0, 0.0, 0.0], + 12 + ) + + def test_rotating_frame_process_split_intervals_match_piecewise_function(self): + piecewise_model, piecewise_model_part = GenerateModel() + split_model, split_model_part = GenerateModel() + + from KratosMultiphysics.MeshMovingApplication.rotating_frame_process import RotatingFrameProcess + + piecewise_parameters = self._GetBaseParameters() + piecewise_parameters["fix_mesh_displacement"].SetBool(False) + piecewise_parameters["fix_velocity"].SetBool(False) + piecewise_parameters.AddEmptyValue("angular_velocity_radians").SetString("-10.0*t if t<1.0 else -10.0") + + split_parameters_1 = self._GetBaseParameters() + split_parameters_1["fix_mesh_displacement"].SetBool(False) + split_parameters_1["fix_velocity"].SetBool(False) + split_parameters_1.AddValue("interval", KratosMultiphysics.Parameters('[0.0, 1.0]')) + split_parameters_1.AddEmptyValue("angular_velocity_radians").SetString("-10.0*t") + + split_parameters_2 = self._GetBaseParameters() + split_parameters_2["fix_mesh_displacement"].SetBool(False) + split_parameters_2["fix_velocity"].SetBool(False) + split_parameters_2.AddValue("interval", KratosMultiphysics.Parameters('[1.0, "End"]')) + split_parameters_2.AddEmptyValue("angular_velocity_radians").SetDouble(-10.0) + + piecewise_process = RotatingFrameProcess(piecewise_model, piecewise_parameters) + split_process_1 = RotatingFrameProcess(split_model, split_parameters_1) + split_process_2 = RotatingFrameProcess(split_model, split_parameters_2) + + # Step at t=1.0 (dt=1.0) + piecewise_process.ExecuteInitializeSolutionStep() + split_process_1.ExecuteInitializeSolutionStep() + split_process_2.ExecuteInitializeSolutionStep() + + # Step at t=2.0 (dt=1.0) + piecewise_model_part.ProcessInfo[KratosMultiphysics.TIME] = 2.0 + split_model_part.ProcessInfo[KratosMultiphysics.TIME] = 2.0 + piecewise_process.ExecuteInitializeSolutionStep() + split_process_1.ExecuteInitializeSolutionStep() + split_process_2.ExecuteInitializeSolutionStep() + + for piecewise_node in piecewise_model_part.Nodes: + split_node = split_model_part.GetNode(piecewise_node.Id) + self.assertVectorAlmostEqual( + piecewise_node.GetSolutionStepValue(KratosMultiphysics.MESH_DISPLACEMENT), + split_node.GetSolutionStepValue(KratosMultiphysics.MESH_DISPLACEMENT), + 8 + ) + self.assertVectorAlmostEqual( + piecewise_node.GetSolutionStepValue(KratosMultiphysics.VELOCITY), + split_node.GetSolutionStepValue(KratosMultiphysics.VELOCITY), + 8 + ) + + +if __name__ == "__main__": + KratosUnittest.main()