From 102d1fb9d97d5e2ee634f4581469b57aea239d80 Mon Sep 17 00:00:00 2001 From: Sebastian Ares de Parga Regalado Date: Wed, 11 Mar 2026 12:13:00 +0100 Subject: [PATCH 1/4] Working with cpp (utility) and python (tensor adaptor). Cpp runs x3 faster. --- .../add_custom_utilities_to_python.cpp | 6 + .../rotating_frame_utility.cpp | 86 ++ .../custom_utilities/rotating_frame_utility.h | 110 +++ .../python_scripts/rotating_frame_process.py | 249 ++++++ .../rotating_frame_process_test.mdpa | 214 +++++ .../rotating_frame_process_test_results.json | 805 ++++++++++++++++++ .../tests/test_rotating_frame_process.py | 102 +++ 7 files changed, 1572 insertions(+) create mode 100644 applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp create mode 100644 applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h create mode 100644 applications/MeshMovingApplication/python_scripts/rotating_frame_process.py create mode 100644 applications/MeshMovingApplication/tests/test_mdpa_files/rotating_frame_process_test.mdpa create mode 100644 applications/MeshMovingApplication/tests/test_results/rotating_frame_process_test_results.json create mode 100644 applications/MeshMovingApplication/tests/test_rotating_frame_process.py diff --git a/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp index a2fbecb49415..9a3b78d8cd79 100644 --- a/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp @@ -27,6 +27,7 @@ #include "custom_utilities/move_mesh_utilities.h" #include "custom_utilities/affine_transform.h" #include "custom_utilities/parametric_affine_transform.h" +#include "custom_utilities/rotating_frame_utility.h" namespace Kratos { namespace Python { @@ -60,6 +61,11 @@ void AddCustomUtilitiesToPython(pybind11::module& m) { .def(py::init()) .def("Apply", &ParametricAffineTransform::Apply) ; + + py::class_(m, "RotatingFrameUtility") + .def_static("AssignRotationalVelocity", &RotatingFrameUtility::AssignRotationalVelocity) + .def_static("ApplyRotationAndMeshDisplacement", &RotatingFrameUtility::ApplyRotationAndMeshDisplacement) + ; void (*CalculateMeshVelocitiesBDF1)(ModelPart&, const TimeDiscretization::BDF1&) = &MeshVelocityCalculation::CalculateMeshVelocities; void (*CalculateMeshVelocitiesBDF2)(ModelPart&, const TimeDiscretization::BDF2&) = &MeshVelocityCalculation::CalculateMeshVelocities; diff --git a/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp b/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp new file mode 100644 index 000000000000..209e75c4af9e --- /dev/null +++ b/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp @@ -0,0 +1,86 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main author: Sebastian Ares de Parga Regalado +// + +// System includes + +// External includes + +// Project includes +#include "custom_utilities/rotating_frame_utility.h" +#include "utilities/parallel_utilities.h" + +namespace Kratos +{ + +void RotatingFrameUtility::AssignRotationalVelocity( + ModelPart& rModelPart, + const array_1d& rAxisOfRotation, + const double Omega, + const array_1d& rCenterOfRotation) +{ + KRATOS_TRY + + const array_1d angular_velocity_vector = Omega * rAxisOfRotation; + + block_for_each(rModelPart.Nodes(), [&](Node& rNode) { + const auto& r_point = rNode.Coordinates(); + const array_1d position_vector = r_point - rCenterOfRotation; + const array_1d velocity_vector = + MathUtils::CrossProduct(angular_velocity_vector, position_vector); + + noalias(rNode.FastGetSolutionStepValue(VELOCITY)) = velocity_vector; + }); + + KRATOS_CATCH(""); +} + +void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( + ModelPart& rModelPart, + const array_1d& rAxisOfRotation, + const double Theta, + const array_1d& rCenterOfRotation) +{ + KRATOS_TRY + + const double a = std::cos(Theta / 2.0); + const double b = -rAxisOfRotation[0] * std::sin(Theta / 2.0); + const double c = -rAxisOfRotation[1] * std::sin(Theta / 2.0); + const double d = -rAxisOfRotation[2] * std::sin(Theta / 2.0); + + BoundedMatrix 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; + + block_for_each(rModelPart.Nodes(), [&](Node& rNode) { + const auto& r_initial_position = rNode.GetInitialPosition().Coordinates(); + const array_1d centered_point = r_initial_position - rCenterOfRotation; + + array_1d rotated_point = prod(centered_point, rot_matrix); + rotated_point += rCenterOfRotation; + + noalias(rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT)) = + rotated_point - rNode.GetInitialPosition(); + + noalias(rNode.Coordinates()) = rotated_point; + }); + + KRATOS_CATCH(""); +} + +} // namespace Kratos \ No newline at end of file diff --git a/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h b/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h new file mode 100644 index 000000000000..1bdb6bb9d539 --- /dev/null +++ b/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h @@ -0,0 +1,110 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main author: Sebastian Ares de Parga Regalado +// + +#pragma once + +// System includes +#include + +// Project includes +#include "includes/define.h" +#include "includes/mesh_moving_variables.h" +#include "utilities/math_utils.h" +#include "utilities/variable_utils.h" + +namespace Kratos +{ + +///@addtogroup MeshMovingApplication +///@{ + +/** + * @class RotatingFrameUtility + * @ingroup MeshMovingApplication + * @brief Utility providing rigid-body rotational kinematics for model parts. + * + * @details + * This utility provides helper functions to: + * + * - Compute and assign rigid-body rotational velocity fields + * - Apply a rigid rotation to a model part based on reference coordinates + * - Compute the corresponding mesh displacement induced by the rotation + * + * The rotation is defined by: + * - an axis of rotation + * - a rotation angle or angular velocity + * - a center of rotation + * + * The utility operates purely at the kinematic level and **does not enforce + * solver policies such as DOF fixity**. Such enforcement should be handled + * at the process level if required. + * + * The implementation is designed to operate efficiently on model part nodes + * and is compatible with Kratos parallel utilities. + */ + +class KRATOS_API(MESH_MOVING_APPLICATION) RotatingFrameUtility +{ +public: + + /** + * @brief Assign rigid-body rotational velocity to nodes. + * + * Computes the velocity induced by rigid-body rotation: + * + * v = ω × r + * + * where: + * - ω is the angular velocity vector + * - r is the position relative to the center of rotation + * + * The computed velocity is assigned to the VELOCITY variable of the nodes + * in the provided model part. + * + * @param rModelPart Model part whose nodes receive the rotational velocity + * @param rAxisOfRotation Axis of rotation + * @param Omega Angular velocity magnitude + * @param rCenterOfRotation Center of rotation + */ + static void AssignRotationalVelocity( + ModelPart& rModelPart, + const array_1d& rAxisOfRotation, + const double Omega, + const array_1d& rCenterOfRotation); + + /** + * @brief Apply rigid rotation and compute mesh displacement. + * + * Applies a rigid-body rotation to the nodes of the provided model part + * using their reference coordinates. The resulting mesh displacement is + * computed as: + * + * MESH_DISPLACEMENT = rotated_position − initial_position + * + * The node coordinates are then updated consistently. + * + * @param rModelPart Model part whose nodes will be rotated + * @param rAxisOfRotation Axis of rotation + * @param Theta Rotation angle + * @param rCenterOfRotation Center of rotation + */ + static void ApplyRotationAndMeshDisplacement( + ModelPart& rModelPart, + const array_1d& rAxisOfRotation, + const double Theta, + const array_1d& rCenterOfRotation); + +}; + +///@} + +} // namespace Kratos \ No newline at end of file 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..ddd42ec4b3de --- /dev/null +++ b/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py @@ -0,0 +1,249 @@ +import KratosMultiphysics as KM +import KratosMultiphysics.MeshMovingApplication as KratosMeshMoving +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 may be applied either directly or after a linear ramp-up + defined by an acceleration time. + """ + + def __init__(self, model, settings): + KM.Process.__init__(self) + + default_settings = KM.Parameters("""{ + "rotating_frame_model_part_name": "", + "rotating_object_model_part_name": "", + "center_of_rotation": [0.0, 0.0, 0.0], + "axis_of_rotation": [1.0, 0.0, 0.0], + "target_angular_velocity_radians": 0.0, + "acceleration_time": 0.0, + "implementation": "cpp", + "echo_level": 0, + "fix_mesh_displacement": false, + "fix_velocity": false + }""") + + 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.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.target_angular_velocity_radians = settings["target_angular_velocity_radians"].GetDouble() + + self.acceleration_time = settings["acceleration_time"].GetDouble() + if self.acceleration_time < 0.0: + raise Exception("The 'acceleration_time' parameter must be non-negative.") + + self.implementation = settings["implementation"].GetString().strip().lower() + if self.implementation in ("c++", "cxx"): + self.implementation = "cpp" + if self.implementation not in ("cpp", "python"): + raise Exception( + f"Unsupported 'implementation': '{self.implementation}'. Use 'cpp' or 'python'." + ) + + 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._python_cache_initialized = False + + def ExecuteInitializeSolutionStep(self): + current_time = self.rotating_frame_model_part.ProcessInfo[KM.TIME] + + theta, omega = self._GetThetaAndOmega(current_time) + + if self.implementation == "cpp": + t0 = time.perf_counter() + KratosMeshMoving.RotatingFrameUtility.ApplyRotationAndMeshDisplacement( + self.rotating_frame_model_part, + self.axis_of_rotation, + theta, + self.center_of_rotation + ) + rot_time = time.perf_counter() - t0 + + t0 = time.perf_counter() + KratosMeshMoving.RotatingFrameUtility.AssignRotationalVelocity( + self.rotating_object_model_part, + self.axis_of_rotation, + omega, + self.center_of_rotation + ) + vel_time = time.perf_counter() - t0 + else: + if not self._python_cache_initialized: + self._InitializePythonBackendCaches() + + t0 = time.perf_counter() + self._ApplyRotationAndMeshDisplacementPython(theta) + rot_time = time.perf_counter() - t0 + + t0 = time.perf_counter() + self._AssignRotationalVelocityPython(omega) + vel_time = time.perf_counter() - t0 + + if self.echo_level > 0: + KM.Logger.PrintInfo( + "RotatingFrameProcess", + f"[implementation={self.implementation}] " + f"rotation={rot_time*1.0e3:.3f} ms, " + f"velocity={vel_time*1.0e3:.3f} ms" + ) + + if self.fix_mesh_displacement: + self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT) + + if self.fix_velocity: + self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY) + + def _GetThetaAndOmega(self, time): + if np.isclose(self.acceleration_time, 0.0): + omega = self.target_angular_velocity_radians + theta = omega * time + return theta, omega + + alpha = self.target_angular_velocity_radians / self.acceleration_time + + if time <= self.acceleration_time: + omega = alpha * time + theta = 0.5 * alpha * time**2 + else: + omega = self.target_angular_velocity_radians + theta = 0.5 * alpha * self.acceleration_time**2 + omega * (time - self.acceleration_time) + + return theta, omega + + 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): + 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, True, model_part.Nodes) + + def _ApplyRotationAndMeshDisplacementPython(self, theta): + sin_half_theta = np.sin(theta / 2.0) + a = np.cos(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 + + current_positions = self._rf_current_pos_ta.data + np.subtract(self._rf_initial_positions, self.center_of_rotation, out=current_positions) + current_positions[:] = current_positions @ rot_matrix + current_positions += self.center_of_rotation + + np.subtract(current_positions, self._rf_initial_positions, out=self._rf_mesh_disp_ta.data) + + 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() + np.subtract( + self._ro_current_pos_ta.data, + self.center_of_rotation, + out=self._ro_relative_positions + ) + + # 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..1c86e466c999 --- /dev/null +++ b/applications/MeshMovingApplication/tests/test_rotating_frame_process.py @@ -0,0 +1,102 @@ +# 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): + + def test_rotating_frame_process(self): + model, model_part = GenerateModel() + + self.print_reference_values = False + + parameters = 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], + "target_angular_velocity_radians": -10.0, + "acceleration_time": 0.0, + "fix_mesh_displacement": true, + "fix_velocity": true + }""") + + 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() + + +if __name__ == "__main__": + KratosUnittest.main() \ No newline at end of file From c9262615238a37bc2c286fa6f156c1e136feb241 Mon Sep 17 00:00:00 2001 From: Sebastian Ares de Parga Regalado Date: Wed, 11 Mar 2026 12:48:44 +0100 Subject: [PATCH 2/4] All code done in python. --- .../add_custom_utilities_to_python.cpp | 7 +- .../rotating_frame_utility.cpp | 86 -------------- .../custom_utilities/rotating_frame_utility.h | 110 ------------------ .../python_scripts/rotating_frame_process.py | 64 +++------- 4 files changed, 16 insertions(+), 251 deletions(-) delete mode 100644 applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp delete mode 100644 applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h diff --git a/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp index 9a3b78d8cd79..2bb064cb1815 100644 --- a/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp @@ -27,7 +27,6 @@ #include "custom_utilities/move_mesh_utilities.h" #include "custom_utilities/affine_transform.h" #include "custom_utilities/parametric_affine_transform.h" -#include "custom_utilities/rotating_frame_utility.h" namespace Kratos { namespace Python { @@ -61,11 +60,7 @@ void AddCustomUtilitiesToPython(pybind11::module& m) { .def(py::init()) .def("Apply", &ParametricAffineTransform::Apply) ; - - py::class_(m, "RotatingFrameUtility") - .def_static("AssignRotationalVelocity", &RotatingFrameUtility::AssignRotationalVelocity) - .def_static("ApplyRotationAndMeshDisplacement", &RotatingFrameUtility::ApplyRotationAndMeshDisplacement) - ; + void (*CalculateMeshVelocitiesBDF1)(ModelPart&, const TimeDiscretization::BDF1&) = &MeshVelocityCalculation::CalculateMeshVelocities; void (*CalculateMeshVelocitiesBDF2)(ModelPart&, const TimeDiscretization::BDF2&) = &MeshVelocityCalculation::CalculateMeshVelocities; diff --git a/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp b/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp deleted file mode 100644 index 209e75c4af9e..000000000000 --- a/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.cpp +++ /dev/null @@ -1,86 +0,0 @@ -// | / | -// ' / __| _` | __| _ \ __| -// . \ | ( | | ( |\__ ` -// _|\_\_| \__,_|\__|\___/ ____/ -// Multi-Physics -// -// License: BSD License -// Kratos default license: kratos/license.txt -// -// Main author: Sebastian Ares de Parga Regalado -// - -// System includes - -// External includes - -// Project includes -#include "custom_utilities/rotating_frame_utility.h" -#include "utilities/parallel_utilities.h" - -namespace Kratos -{ - -void RotatingFrameUtility::AssignRotationalVelocity( - ModelPart& rModelPart, - const array_1d& rAxisOfRotation, - const double Omega, - const array_1d& rCenterOfRotation) -{ - KRATOS_TRY - - const array_1d angular_velocity_vector = Omega * rAxisOfRotation; - - block_for_each(rModelPart.Nodes(), [&](Node& rNode) { - const auto& r_point = rNode.Coordinates(); - const array_1d position_vector = r_point - rCenterOfRotation; - const array_1d velocity_vector = - MathUtils::CrossProduct(angular_velocity_vector, position_vector); - - noalias(rNode.FastGetSolutionStepValue(VELOCITY)) = velocity_vector; - }); - - KRATOS_CATCH(""); -} - -void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( - ModelPart& rModelPart, - const array_1d& rAxisOfRotation, - const double Theta, - const array_1d& rCenterOfRotation) -{ - KRATOS_TRY - - const double a = std::cos(Theta / 2.0); - const double b = -rAxisOfRotation[0] * std::sin(Theta / 2.0); - const double c = -rAxisOfRotation[1] * std::sin(Theta / 2.0); - const double d = -rAxisOfRotation[2] * std::sin(Theta / 2.0); - - BoundedMatrix 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; - - block_for_each(rModelPart.Nodes(), [&](Node& rNode) { - const auto& r_initial_position = rNode.GetInitialPosition().Coordinates(); - const array_1d centered_point = r_initial_position - rCenterOfRotation; - - array_1d rotated_point = prod(centered_point, rot_matrix); - rotated_point += rCenterOfRotation; - - noalias(rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT)) = - rotated_point - rNode.GetInitialPosition(); - - noalias(rNode.Coordinates()) = rotated_point; - }); - - KRATOS_CATCH(""); -} - -} // namespace Kratos \ No newline at end of file diff --git a/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h b/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h deleted file mode 100644 index 1bdb6bb9d539..000000000000 --- a/applications/MeshMovingApplication/custom_utilities/rotating_frame_utility.h +++ /dev/null @@ -1,110 +0,0 @@ -// | / | -// ' / __| _` | __| _ \ __| -// . \ | ( | | ( |\__ ` -// _|\_\_| \__,_|\__|\___/ ____/ -// Multi-Physics -// -// License: BSD License -// Kratos default license: kratos/license.txt -// -// Main author: Sebastian Ares de Parga Regalado -// - -#pragma once - -// System includes -#include - -// Project includes -#include "includes/define.h" -#include "includes/mesh_moving_variables.h" -#include "utilities/math_utils.h" -#include "utilities/variable_utils.h" - -namespace Kratos -{ - -///@addtogroup MeshMovingApplication -///@{ - -/** - * @class RotatingFrameUtility - * @ingroup MeshMovingApplication - * @brief Utility providing rigid-body rotational kinematics for model parts. - * - * @details - * This utility provides helper functions to: - * - * - Compute and assign rigid-body rotational velocity fields - * - Apply a rigid rotation to a model part based on reference coordinates - * - Compute the corresponding mesh displacement induced by the rotation - * - * The rotation is defined by: - * - an axis of rotation - * - a rotation angle or angular velocity - * - a center of rotation - * - * The utility operates purely at the kinematic level and **does not enforce - * solver policies such as DOF fixity**. Such enforcement should be handled - * at the process level if required. - * - * The implementation is designed to operate efficiently on model part nodes - * and is compatible with Kratos parallel utilities. - */ - -class KRATOS_API(MESH_MOVING_APPLICATION) RotatingFrameUtility -{ -public: - - /** - * @brief Assign rigid-body rotational velocity to nodes. - * - * Computes the velocity induced by rigid-body rotation: - * - * v = ω × r - * - * where: - * - ω is the angular velocity vector - * - r is the position relative to the center of rotation - * - * The computed velocity is assigned to the VELOCITY variable of the nodes - * in the provided model part. - * - * @param rModelPart Model part whose nodes receive the rotational velocity - * @param rAxisOfRotation Axis of rotation - * @param Omega Angular velocity magnitude - * @param rCenterOfRotation Center of rotation - */ - static void AssignRotationalVelocity( - ModelPart& rModelPart, - const array_1d& rAxisOfRotation, - const double Omega, - const array_1d& rCenterOfRotation); - - /** - * @brief Apply rigid rotation and compute mesh displacement. - * - * Applies a rigid-body rotation to the nodes of the provided model part - * using their reference coordinates. The resulting mesh displacement is - * computed as: - * - * MESH_DISPLACEMENT = rotated_position − initial_position - * - * The node coordinates are then updated consistently. - * - * @param rModelPart Model part whose nodes will be rotated - * @param rAxisOfRotation Axis of rotation - * @param Theta Rotation angle - * @param rCenterOfRotation Center of rotation - */ - static void ApplyRotationAndMeshDisplacement( - ModelPart& rModelPart, - const array_1d& rAxisOfRotation, - const double Theta, - const array_1d& rCenterOfRotation); - -}; - -///@} - -} // namespace Kratos \ No newline at end of file diff --git a/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py b/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py index ddd42ec4b3de..4847f4e69f52 100644 --- a/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py +++ b/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py @@ -1,5 +1,4 @@ import KratosMultiphysics as KM -import KratosMultiphysics.MeshMovingApplication as KratosMeshMoving import numpy as np import time @@ -31,7 +30,6 @@ def __init__(self, model, settings): "axis_of_rotation": [1.0, 0.0, 0.0], "target_angular_velocity_radians": 0.0, "acceleration_time": 0.0, - "implementation": "cpp", "echo_level": 0, "fix_mesh_displacement": false, "fix_velocity": false @@ -76,14 +74,6 @@ def __init__(self, model, settings): if self.acceleration_time < 0.0: raise Exception("The 'acceleration_time' parameter must be non-negative.") - self.implementation = settings["implementation"].GetString().strip().lower() - if self.implementation in ("c++", "cxx"): - self.implementation = "cpp" - if self.implementation not in ("cpp", "python"): - raise Exception( - f"Unsupported 'implementation': '{self.implementation}'. Use 'cpp' or 'python'." - ) - self.echo_level = settings["echo_level"].GetInt() if self.echo_level < 0: raise Exception("The 'echo_level' parameter must be >= 0.") @@ -97,42 +87,20 @@ def ExecuteInitializeSolutionStep(self): theta, omega = self._GetThetaAndOmega(current_time) - if self.implementation == "cpp": - t0 = time.perf_counter() - KratosMeshMoving.RotatingFrameUtility.ApplyRotationAndMeshDisplacement( - self.rotating_frame_model_part, - self.axis_of_rotation, - theta, - self.center_of_rotation - ) - rot_time = time.perf_counter() - t0 - - t0 = time.perf_counter() - KratosMeshMoving.RotatingFrameUtility.AssignRotationalVelocity( - self.rotating_object_model_part, - self.axis_of_rotation, - omega, - self.center_of_rotation - ) - vel_time = time.perf_counter() - t0 - else: - if not self._python_cache_initialized: - self._InitializePythonBackendCaches() + if not self._python_cache_initialized: + self._InitializePythonBackendCaches() - t0 = time.perf_counter() - self._ApplyRotationAndMeshDisplacementPython(theta) - rot_time = time.perf_counter() - t0 + t0 = time.perf_counter() + self._ApplyRotationAndMeshDisplacementPython(theta) + rot_time = time.perf_counter() - t0 - t0 = time.perf_counter() - self._AssignRotationalVelocityPython(omega) - vel_time = time.perf_counter() - t0 + t0 = time.perf_counter() + self._AssignRotationalVelocityPython(omega) + vel_time = time.perf_counter() - t0 if self.echo_level > 0: KM.Logger.PrintInfo( - "RotatingFrameProcess", - f"[implementation={self.implementation}] " - f"rotation={rot_time*1.0e3:.3f} ms, " - f"velocity={vel_time*1.0e3:.3f} ms" + "RotatingFrameProcess", f" rotation={rot_time*1.0e3:.3f} ms, velocity={vel_time*1.0e3:.3f} ms" ) if self.fix_mesh_displacement: @@ -219,12 +187,14 @@ def _ApplyRotationAndMeshDisplacementPython(self, theta): rot_matrix[2, 1] = 2.0 * (c * d + a * b) rot_matrix[2, 2] = a * a + d * d - b * b - c * c + # Rotate reference coordinates around the center: + # x = (x0 - c) @ R + c current_positions = self._rf_current_pos_ta.data - np.subtract(self._rf_initial_positions, self.center_of_rotation, out=current_positions) - current_positions[:] = current_positions @ rot_matrix + current_positions[:] = (self._rf_initial_positions - self.center_of_rotation) @ rot_matrix current_positions += self.center_of_rotation - np.subtract(current_positions, self._rf_initial_positions, out=self._rf_mesh_disp_ta.data) + # 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() @@ -233,11 +203,7 @@ def _AssignRotationalVelocityPython(self, omega): angular_velocity_vector = omega * self.axis_of_rotation self._ro_current_pos_ta.CollectData() - np.subtract( - self._ro_current_pos_ta.data, - self.center_of_rotation, - out=self._ro_relative_positions - ) + 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 From 1cdb78bceb46b20dba47a94336002fd0bbe5fdee Mon Sep 17 00:00:00 2001 From: Sebastian Ares de Parga Regalado Date: Wed, 11 Mar 2026 12:51:12 +0100 Subject: [PATCH 3/4] Detail --- .../custom_python/add_custom_utilities_to_python.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp index 2bb064cb1815..a2fbecb49415 100644 --- a/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/MeshMovingApplication/custom_python/add_custom_utilities_to_python.cpp @@ -61,7 +61,6 @@ void AddCustomUtilitiesToPython(pybind11::module& m) { .def("Apply", &ParametricAffineTransform::Apply) ; - void (*CalculateMeshVelocitiesBDF1)(ModelPart&, const TimeDiscretization::BDF1&) = &MeshVelocityCalculation::CalculateMeshVelocities; void (*CalculateMeshVelocitiesBDF2)(ModelPart&, const TimeDiscretization::BDF2&) = &MeshVelocityCalculation::CalculateMeshVelocities; void (*CalculateMeshVelocitiesNewmark)(ModelPart&, const TimeDiscretization::Newmark&) = &MeshVelocityCalculation::CalculateMeshVelocities; From 2a4ce30e041e8cfa4730b5b1a5905e22b835414b Mon Sep 17 00:00:00 2001 From: Sebastian Ares de Parga Regalado Date: Thu, 19 Mar 2026 14:10:56 +0100 Subject: [PATCH 4/4] Refactor to function-based, interval-aware angular velocity MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Commit body Replace legacy target_angular_velocity_radians + acceleration_time input with a single angular_velocity_radians (constant or f(x,y,z,t) string). Add standard interval support to align with Kratos process conventions. Change rotation update to incremental step integration (delta_theta) over the overlap: [t-Δt, t] ∩ [interval_begin, interval_end] Apply incremental rotation on current coordinates and recompute MESH_DISPLACEMENT from initial coordinates. Keep rotational VELOCITY assignment active only when process interval is active. Manage fixity activation/deactivation consistently when entering/leaving interval. Enable equivalent behavior between: one process with piecewise angular velocity function, and two process instances split by interval. Update/add tests for: constant numeric vs constant string input equivalence, no-op outside interval, split-interval vs single piecewise-function equivalence. Migration note Update project parameters from: target_angular_velocity_radians, acceleration_time To: angular_velocity_radians (+ optional interval). --- .../python_scripts/rotating_frame_process.py | 128 +++++++++++------ .../tests/test_rotating_frame_process.py | 129 ++++++++++++++++-- 2 files changed, 209 insertions(+), 48 deletions(-) diff --git a/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py b/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py index 4847f4e69f52..d1beb7e060ee 100644 --- a/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py +++ b/applications/MeshMovingApplication/python_scripts/rotating_frame_process.py @@ -16,8 +16,10 @@ class RotatingFrameProcess(KM.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 may be applied either directly or after a linear ramp-up - defined by an acceleration time. + 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): @@ -26,15 +28,22 @@ def __init__(self, model, settings): 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], - "target_angular_velocity_radians": 0.0, - "acceleration_time": 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 @@ -50,6 +59,8 @@ def __init__(self, model, settings): 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) @@ -68,11 +79,10 @@ def __init__(self, model, settings): ) self.axis_of_rotation /= axis_norm - self.target_angular_velocity_radians = settings["target_angular_velocity_radians"].GetDouble() - - self.acceleration_time = settings["acceleration_time"].GetDouble() - if self.acceleration_time < 0.0: - raise Exception("The 'acceleration_time' parameter must be non-negative.") + 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: @@ -80,22 +90,34 @@ def __init__(self, model, settings): 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] - - theta, omega = self._GetThetaAndOmega(current_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() - self._ApplyRotationAndMeshDisplacementPython(theta) + if not np.isclose(delta_theta, 0.0): + self._ApplyRotationAndMeshDisplacementPython(delta_theta) rot_time = time.perf_counter() - t0 t0 = time.perf_counter() - self._AssignRotationalVelocityPython(omega) + if is_active: + self._AssignRotationalVelocityPython(omega) vel_time = time.perf_counter() - t0 if self.echo_level > 0: @@ -103,28 +125,56 @@ def ExecuteInitializeSolutionStep(self): "RotatingFrameProcess", f" rotation={rot_time*1.0e3:.3f} ms, velocity={vel_time*1.0e3:.3f} ms" ) - if self.fix_mesh_displacement: - self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT) + 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) + if self.fix_velocity: + self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY, True) - def _GetThetaAndOmega(self, time): - if np.isclose(self.acceleration_time, 0.0): - omega = self.target_angular_velocity_radians - theta = omega * time - return theta, omega + 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 - alpha = self.target_angular_velocity_radians / self.acceleration_time - - if time <= self.acceleration_time: - omega = alpha * time - theta = 0.5 * alpha * time**2 + @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: - omega = self.target_angular_velocity_radians - theta = 0.5 * alpha * self.acceleration_time**2 + omega * (time - self.acceleration_time) - - return theta, omega + 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( @@ -156,7 +206,7 @@ def _InitializePythonBackendCaches(self): self._python_cache_initialized = True @staticmethod - def _FixVectorVariable(model_part, variable): + 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), @@ -167,11 +217,11 @@ def _FixVectorVariable(model_part, variable): vu = KM.VariableUtils() for comp in components_map[variable]: - vu.ApplyFixity(comp, True, model_part.Nodes) + vu.ApplyFixity(comp, is_fixed, model_part.Nodes) - def _ApplyRotationAndMeshDisplacementPython(self, theta): - sin_half_theta = np.sin(theta / 2.0) - a = np.cos(theta / 2.0) + 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 @@ -187,10 +237,10 @@ def _ApplyRotationAndMeshDisplacementPython(self, theta): rot_matrix[2, 1] = 2.0 * (c * d + a * b) rot_matrix[2, 2] = a * a + d * d - b * b - c * c - # Rotate reference coordinates around the center: - # x = (x0 - c) @ R + 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[:] = (self._rf_initial_positions - self.center_of_rotation) @ rot_matrix + 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 diff --git a/applications/MeshMovingApplication/tests/test_rotating_frame_process.py b/applications/MeshMovingApplication/tests/test_rotating_frame_process.py index 1c86e466c999..852da88837d5 100644 --- a/applications/MeshMovingApplication/tests/test_rotating_frame_process.py +++ b/applications/MeshMovingApplication/tests/test_rotating_frame_process.py @@ -43,23 +43,26 @@ def GenerateModel(): @KratosUnittest.skipUnless(mesh_moving_is_available, "MeshMovingApplication is not available") class TestRotatingFrameProcess(KratosUnittest.TestCase): - def test_rotating_frame_process(self): - model, model_part = GenerateModel() - - self.print_reference_values = False - - parameters = KratosMultiphysics.Parameters(""" + @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], - "target_angular_velocity_radians": -10.0, - "acceleration_time": 0.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) @@ -97,6 +100,114 @@ def test_rotating_frame_process(self): 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() \ No newline at end of file + KratosUnittest.main()