From 599231acaf945d165950d51992238a6a3834f27c Mon Sep 17 00:00:00 2001 From: SADPR Date: Wed, 17 May 2023 18:48:07 +0200 Subject: [PATCH 01/20] First commit. --- .../FluidDynamicsApplication/CMakeLists.txt | 3 +- .../add_custom_utilities_to_python.cpp | 7 + .../rotating_frame_utility.cpp | 185 ++++++++++ .../custom_utilities/rotating_frame_utility.h | 109 ++++++ .../python_scripts/rotating_frame_process.py | 348 ++++++++++++++++++ 5 files changed, 651 insertions(+), 1 deletion(-) create mode 100644 applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp create mode 100644 applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h create mode 100644 applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py diff --git a/applications/FluidDynamicsApplication/CMakeLists.txt b/applications/FluidDynamicsApplication/CMakeLists.txt index 462d47a49b86..6fbf0e0d6fe6 100644 --- a/applications/FluidDynamicsApplication/CMakeLists.txt +++ b/applications/FluidDynamicsApplication/CMakeLists.txt @@ -12,6 +12,7 @@ set( KRATOS_FLUID_DYNAMICS_APPLICATION_CORE_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/fluid_dynamics_application_variables.cpp # utilities (compiled first because they are used within the element cpps) + ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/rotating_frame_utility.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/estimate_dt_utilities.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/fluid_auxiliary_utilities.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/fluid_calculation_utilities.cpp @@ -176,4 +177,4 @@ endif((${USE_MPI} MATCHES ON) AND (${TRILINOS_FOUND})) # Define custom targets set(KRATOS_KERNEL "${KRATOS_KERNEL};KratosFluidDynamicsCore" PARENT_SCOPE) -set(KRATOS_PYTHON_INTERFACE "${KRATOS_PYTHON_INTERFACE};KratosFluidDynamicsApplication" PARENT_SCOPE) \ No newline at end of file +set(KRATOS_PYTHON_INTERFACE "${KRATOS_PYTHON_INTERFACE};KratosFluidDynamicsApplication" PARENT_SCOPE) diff --git a/applications/FluidDynamicsApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/FluidDynamicsApplication/custom_python/add_custom_utilities_to_python.cpp index cffb73f30bea..3ec8aa8a190b 100644 --- a/applications/FluidDynamicsApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/FluidDynamicsApplication/custom_python/add_custom_utilities_to_python.cpp @@ -36,6 +36,7 @@ #include "custom_utilities/compressible_element_rotation_utility.h" #include "custom_utilities/acceleration_limitation_utilities.h" #include "custom_utilities/fluid_test_utilities.h" +#include "custom_utilities/rotating_frame_utility.h" #include "utilities/split_tetrahedra.h" @@ -199,6 +200,12 @@ void AddCustomUtilitiesToPython(pybind11::module& m) .def_static("RandomFillNonHistoricalVariable", py::overload_cast>&, const IndexType, const double, const double>(&FluidTestUtilities::RandomFillNonHistoricalVariable>)) .def_static("RandomFillNonHistoricalVariable", py::overload_cast>&, const std::string&, const IndexType, const double, const double>(&FluidTestUtilities::RandomFillNonHistoricalVariable>)) ; + + // Rotating Frame Utility + py::class_(m, "RotatingFrameUtility") + .def_static("ApplyVelocityToRotatingObject", &RotatingFrameUtility::ApplyVelocityToRotatingObject) + .def_static("ApplyRotationAndMeshDisplacement", &RotatingFrameUtility::ApplyRotationAndMeshDisplacement) + ; } } // namespace Python. diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp new file mode 100644 index 000000000000..6e82d1151b92 --- /dev/null +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -0,0 +1,185 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main author: Sebastian Ares de Parga Regalado +// + +// System includes +#include +#include + +// External includes + +// Application includes +#include "rotating_frame_utility.h" + +namespace Kratos +{ + +void RotatingFrameUtility::ApplyVelocityToRotatingObject( + ModelPart& rRotatingObjectModelPart, + const double Omega, + const array_1d& rAxisOfRotation) +{ + KRATOS_TRY; + + // Get the array of nodes in the rotating object model part + auto& nodes_array = rRotatingObjectModelPart.Nodes(); + + // Declare thread-local variables outside the block_for_each + thread_local double radius; + thread_local double velocity_magnitude; + thread_local array_1d velocity_direction; + thread_local array_1d velocity_vector; + thread_local array_1d velocity; + + // Apply the velocity calculations to each node in parallel + block_for_each(nodes_array, [&](Node& rNode) { + radius = std::sqrt(rNode.X() * rNode.X() + rNode.Y() * rNode.Y() + rNode.Z() * rNode.Z()); + + // Apply component by component the corresponding velocity. + velocity_magnitude = Omega * radius; + velocity_direction = MathUtils::CrossProduct(rAxisOfRotation, rNode.Coordinates()); + velocity_vector = velocity_direction / norm_2(velocity_direction); + velocity = velocity_magnitude * velocity_vector; + + // Set the solution step values for the velocity components using FastGetSolutionStepValue + rNode.FastGetSolutionStepValue(VELOCITY_X) = velocity[0]; + rNode.FastGetSolutionStepValue(VELOCITY_Y) = velocity[1]; + rNode.FastGetSolutionStepValue(VELOCITY_Z) = velocity[2]; + + // Fix the velocity components to ensure they remain constant during the simulation + rNode.Fix(VELOCITY_X); + rNode.Fix(VELOCITY_Y); + rNode.Fix(VELOCITY_Z); + }); + + KRATOS_CATCH(""); +} + +void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( + ModelPart& rRotatingFrameModelPart, + const array_1d& rAxis, + const double& rTheta, + const array_1d& rCenter) +{ + KRATOS_TRY; + + // Quaternion constants + const double a = std::cos(rTheta / 2); + const double b = -rAxis[0] * std::sin(rTheta / 2); + const double c = -rAxis[1] * std::sin(rTheta / 2); + const double d = -rAxis[2] * std::sin(rTheta / 2); + + // Creating a quaternion rotation matrix + BoundedMatrix rot_matrix; + rot_matrix(0, 0) = a*a+b*b-c*c-d*d; + rot_matrix(0, 1) = 2*(b*c-a*d); + rot_matrix(0, 2) = 2*(b*d+a*c); + rot_matrix(1, 0) = 2*(b*c+a*d); + rot_matrix(1, 1) = a*a+c*c-b*b-d*d; + rot_matrix(1, 2) = 2*(c*d-a*b); + rot_matrix(2, 0) = 2*(b*d-a*c); + rot_matrix(2, 1) = 2*(c*d+a*b); + rot_matrix(2, 2) = a*a+d*d-b*b-c*c; + + // Apply the rotation and mesh displacement calculations to each node in parallel + block_for_each(rRotatingFrameModelPart.Nodes(), [&](Node& rNode) { + // Getting the initial coordinates of the node + auto& point = rNode.GetInitialPosition().Coordinates(); + // point[0] = node.X0(); + // point[0] = node.X0(); + // point[0] = node.X0(); + + // Shifting the rotation center to the origin + point -= rCenter; + + // Applying the rotation + array_1d rotated_point = prod(rot_matrix, point); + + // Shifting the point back and updating the node coordinates + rotated_point += rCenter; + rNode.FastGetSolutionStepValue(DISPLACEMENT_X) = rotated_point[0] - rNode.X0(); + rNode.FastGetSolutionStepValue(DISPLACEMENT_Y) = rotated_point[1] - rNode.Y0(); + rNode.FastGetSolutionStepValue(DISPLACEMENT_Z) = rotated_point[2] - rNode.Z0(); + rNode.Fix(DISPLACEMENT_X); + rNode.Fix(DISPLACEMENT_Y); + rNode.Fix(DISPLACEMENT_Z); + rNode.X() = rotated_point[0]; + rNode.Y() = rotated_point[1]; + rNode.Z() = rotated_point[2]; + }); + + KRATOS_CATCH(""); +} + + +// void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( +// ModelPart& rRotatingFrameModelPart, +// const Matrix& rTransformationMatrix) +// { +// KRATOS_TRY; + +// // Get the array of nodes in the rotating frame model part +// auto& nodes_array = rRotatingFrameModelPart.Nodes(); + +// // Declare thread-local variables outside the block_for_each +// thread_local array_1d homogeneous_coordinates; +// thread_local array_1d rotated_coordinates; +// thread_local array_1d rotated_coor; +// thread_local double mesh_displacement_x; +// thread_local double mesh_displacement_y; +// thread_local double mesh_displacement_z; + +// // Apply the rotation and mesh displacement calculations to each node in parallel +// block_for_each(nodes_array, [&](Node& rNode) { + +// // Convert point to homogeneous coordinates +// // homogeneous_coordinates.resize(4); +// homogeneous_coordinates[0] = rNode.X0(); +// homogeneous_coordinates[1] = rNode.Y0(); +// homogeneous_coordinates[2] = rNode.Z0(); +// homogeneous_coordinates[3] = 1.0; + +// // Apply transformation to homogeneous coordinates +// // rotated_coordinates.resize(4); +// noalias(rotated_coordinates) = prod(rTransformationMatrix, homogeneous_coordinates); + +// // Convert back to 3D coordinates +// // rotated_coor.resize(3); +// rotated_coor[0] = rotated_coordinates[0] / rotated_coordinates[3]; +// rotated_coor[1] = rotated_coordinates[1] / rotated_coordinates[3]; +// rotated_coor[2] = rotated_coordinates[2] / rotated_coordinates[3]; + +// // Calculate mesh displacement +// mesh_displacement_x = rotated_coor[0] - rNode.X0(); +// mesh_displacement_y = rotated_coor[1] - rNode.Y0(); +// mesh_displacement_z = rotated_coor[2] - rNode.Z0(); + +// // Set the solution step values for mesh displacement using FastGetSolutionStepValue +// rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_X) = mesh_displacement_x; +// rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Y) = mesh_displacement_y; +// rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Z) = mesh_displacement_z; + +// // Fix the mesh displacement components to ensure they remain constant during the simulation +// rNode.Fix(MESH_DISPLACEMENT_X); +// rNode.Fix(MESH_DISPLACEMENT_Y); +// rNode.Fix(MESH_DISPLACEMENT_Z); + +// // Update the coordinates of the node using FastGetSolutionStepValue +// rNode.X() = rotated_coor[0]; +// rNode.Y() = rotated_coor[1]; +// rNode.Z() = rotated_coor[2]; +// }); + +// KRATOS_CATCH(""); +// } + +} // namespace Kratos. + diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h new file mode 100644 index 000000000000..62d9b6011f7b --- /dev/null +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h @@ -0,0 +1,109 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main author: Sebastian Ares de Parga Regalado +// + +#if !defined(KRATOS_ROTATING_FRAME_UTILITY_H) +#define KRATOS_ROTATING_FRAME_UTILITY_H + +// System includes +#include +#include + +// External includes +#include "includes/mesh_moving_variables.h" + +// Include kratos definitions +#include "includes/define.h" + +// Application includes +#include "utilities/variable_utils.h" + +// Project includes +#include "utilities/builtin_timer.h" +#include "utilities/math_utils.h" + +namespace Kratos { +///@addtogroup FluidDynamicsApplication +///@{ + +///@name Kratos classes +///@{ + +/** + * @class RotateFrameUtility + * @ingroup Kratos Core + * @brief Utility for rotating a frame with a rotating object. + * @details A utility that rotates a frame (model part) with a rotating object (solid) + * based on the specified angular velocity and axis of rotation. It enables the + * transformation of coordinates and displacements between the rotating and + * non-rotating frames, ensuring consistent updates of the frame's entities + * during simulation. Designed for use within the Kratos framework for various + * applications such as sliding mesh problems and rotating machinery simulations. + * @note Thread-safe and supports parallel execution using Kratos block for each. + * @author Sebastian Ares de Parga Regalado +*/ + +class KRATOS_API(FLUID_DYNAMICS_APPLICATION) RotatingFrameUtility +{ + public: + ///@name Type Definitions + ///@{ + + ///@} + ///@name Static Operations + ///@{ + + /** + * @brief Apply velocity to rotating object. + * + * This function applies velocity to a rotating object model part. + * The velocity is computed based on the angular velocity and the + * axis of rotation. The function calculates the radius of each node, + * and then applies the corresponding velocity component-wise. The + * velocity components are fixed to remain constant during the simulation, + * and the solution step values for the velocity components are set + * accordingly. + * + * @param rRotatingObjectModelPart The rotating object model part. + * @param Omega The angular velocity. + * @param rAxisOfRotation The axis of rotation. + */ + static void ApplyVelocityToRotatingObject( + ModelPart& rRotatingObjectModelPart, + const double Omega, + const array_1d& rAxisOfRotation); + + /** + * @brief Apply rotation and mesh displacement to a rotating frame. + * This function applies rotation and mesh displacement to a rotating frame model part based on the provided rotation axis, angle and rotation center. + * It constructs a rotation matrix based on the rotation axis and angle. The rotation is applied to each node, which is then displaced + * by the difference between the rotated coordinates and the initial coordinates. + * The solution step values for the mesh displacement components are set accordingly, and the mesh displacement components + * are fixed to remain constant during the simulation. Finally, the node coordinates are updated with the rotated + * coordinates. + * @param rRotatingFrameModelPart The rotating frame model part to which the rotation and mesh displacement will be applied. + * @param rAxis The rotation axis vector. + * @param rTheta The rotation angle. + * @param rCenter The center of rotation. + */ + static void ApplyRotationAndMeshDisplacement( + ModelPart& rRotatingFrameModelPart, + const array_1d& rAxis, + const double& rTheta, + const array_1d& rCenter); + +}; // Class RotatingFrameUtility + +///@} + +} // namespace Kratos. + +#endif // KRATOS_ROTATING_FRAME_UTILITY_H \ No newline at end of file diff --git a/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py new file mode 100644 index 000000000000..507f651e25c2 --- /dev/null +++ b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py @@ -0,0 +1,348 @@ +import KratosMultiphysics as KM +import KratosMultiphysics.FluidDynamicsApplication as KratosCFD +import numpy as np + +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"]) + +## All the processes python should be derived from "Process" +class RotatingFrameProcess(KM.Process): + """A Kratos process that rotates a given RotatingFrameModelPart + around a specified axis of rotation. This process uses the Kratos + RotateNodesUtility to rotate the nodes of the RotatingFrameModelPart + and assigns the corresponding MeshDisplacement to the rotated nodes. + The process also assigns the corresponding rotational mesh velocity to + the nodes in a RotatingObjectModelPart. The process takes as input the + model part to be rotated, the axis of rotation vector, the final angular + velocity in radians per second, and the accelerating time in seconds. + + Public member variables: + model -- the container of the different model parts. + settings -- Kratos parameters containing process settings. + """ + + def __init__(self, model, settings): + """ The default constructor of the class + + Keyword arguments: + self -- It signifies an instance of a class. + Model -- the container of the different model parts. + settings -- Kratos parameters containing process settings. + """ + KM.Process.__init__(self) # calling the baseclass constructor + + 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], + "angular_velocity_radians": 0.0, + "acceleration_time": 0.0 + }""") + + # Add missing settings that the user did not provide but that + # are necessary for this process + settings.ValidateAndAssignDefaults(default_settings) + # Alternative: + # settings.RecursivelyValidateAndAssignDefaults(default_settings) + + #Assign settings + self.model = model + + # Get the rotating frame model part name + if not settings["rotating_frame_model_part_name"].GetString(): + raise Exception("\'rotating_frame_model_part_name\' not provided. Please specify the model part to rotate the frame.") + self.rotating_frame_model_part_name = settings["rotating_frame_model_part_name"].GetString() + + # Get the rotating object model part name + if not settings["rotating_object_model_part_name"].GetString(): + raise Exception("\'rotating_object_model_part_name\' not provided. Please specify the slave model part.") + self.rotating_object_model_part_name = settings["rotating_object_model_part_name"].GetString() + + ## Assign rotating frame and object model parts + 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) + + # Get the center of rotation + if settings.Has("center_of_rotation"): + self.center_of_rotation = np.array(settings["center_of_rotation"].GetVector()) + else: + raise Exception("The center_of_rotation parameter is missing from the settings.") + + # Get the axis of rotation + if settings.Has("axis_of_rotation"): + self.axis_of_rotation = np.array(settings["axis_of_rotation"].GetVector()) + 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 not np.isclose(np.linalg.norm(axis_norm), 1.0, rtol=1e-6): + KM.Logger.PrintWarning("RotatingFrameProcess", "The axis_of_rotation vector is not a unit vector... normalizing the vector.") + self.axis_of_rotation = self.axis_of_rotation / axis_norm + else: + raise Exception("The axis_of_rotation parameter is missing from the settings.") + + # Get the angular velocity in radians + if settings.Has("angular_velocity_radians"): + self.angular_velocity_radians = settings["angular_velocity_radians"].GetDouble() + else: + raise Exception("The angular_velocity_radians parameter is missing from the settings.") + + # Get the acceleration time + if settings.Has("acceleration_time"): + self.acceleration_time = settings["acceleration_time"].GetDouble() + if self.acceleration_time < 0.0: + raise Exception("The acceleration_time parameter must be non-negative.") + else: + raise Exception("The acceleration_time parameter is missing from the settings.") + + def ExecuteInitialize(self): + """ This method is executed at the begining to initialize the process + + Keyword arguments: + self -- It signifies an instance of a class. + """ + pass + + def Check(self): + """ This method verifies that the input is correct + + Keyword arguments: + self -- It signifies an instance of a class. + """ + pass + + def ExecuteBeforeSolutionLoop(self): + """ This method is executed just before the solution-loop + + Keyword arguments: + self -- It signifies an instance of a class. + """ + pass + + def ExecuteInitializeSolutionStep(self): + """ This method is executed in order to initialize the current step + + Keyword arguments: + self -- It signifies an instance of a class. + """ + # Get Current time + self.time = self.rotating_frame_model_part.ProcessInfo[KM.TIME] + + # Calculate the angle (theta) for the specified time + self.get_theta_and_omega() #Potentially python + + # Calculate the rotation matrix. + # transformation_matrix = self.get_rotation_matrix(self.theta) #Potentially C++ + + # Apply Rotation and Mesh Displacement to nodes. + # self.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, transformation_matrix) #Potentially C++ + # KratosCFD.RotatingFrameUtility.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, transformation_matrix) + self.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, self.axis_of_rotation, self.theta, self.center_of_rotation) + + # Apply Velocity to object. + self.set_node_velocities(self.rotating_object_model_part, self.axis_of_rotation, self.omega, self.center_of_rotation) + # self.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.omega, self.center_of_rotation, self.axis_of_rotation) + # KratosCFD.RotatingFrameUtility.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.omega, self.axis_of_rotation) + + def ApplyVelocityToRotatingObject(self, rotating_object_model_part, omega, center_of_rotation, axis_of_rotation): + for node in rotating_object_model_part.Nodes: + relative_coordinates = [node.X, node.Y, node.Z]-center_of_rotation + radius = np.linalg.norm(relative_coordinates) + #Apply component by component the corresponding velocity. + velocity_magnitude = omega * radius + velocity_direction = np.cross(axis_of_rotation, relative_coordinates) + velocity_vector = velocity_direction / np.linalg.norm(velocity_direction) + velocity = velocity_magnitude * velocity_vector + node.SetSolutionStepValue(KM.VELOCITY_X, velocity[0]) + node.SetSolutionStepValue(KM.VELOCITY_Y, velocity[1]) + node.SetSolutionStepValue(KM.VELOCITY_Z, velocity[2]) + node.Fix(KM.VELOCITY_X) + node.Fix(KM.VELOCITY_Y) + node.Fix(KM.VELOCITY_Z) + + def calculate_magnitude_difference(self, point, center_of_rotation, axis_of_rotation): + # Step 1: Calculate the vector from the target point to the given point + vector_difference = point - center_of_rotation + + # Step 2: Calculate the dot product between the vector difference and the target axis + dot_product = np.dot(vector_difference, axis_of_rotation) + + # Step 3: Calculate the magnitude of the target axis + magnitude_axis_of_rotation = np.linalg.norm(axis_of_rotation) + + # Step 4: Divide the dot product by the magnitude of the target axis to obtain the projection + projection = dot_product / magnitude_axis_of_rotation + + # Step 5: Calculate the magnitude of the projection to obtain the desired magnitude of difference + magnitude_difference = np.abs(projection) + + return magnitude_difference + + # def ApplyVelocityToRotatingObject(self, rotating_object_model_part, omega, center_of_rotation, axis_of_rotation): + # for node in rotating_object_model_part.Nodes: + # # Calculate the magnitude difference + # point = np.array([node.X, node.Y, node.Z]) + # magnitude_difference = self.calculate_magnitude_difference(point, center_of_rotation, axis_of_rotation) + + # # Apply component by component the corresponding velocity. + # velocity_magnitude = omega * magnitude_difference + # velocity_direction = np.cross(axis_of_rotation, point) + # velocity_vector = velocity_direction / np.linalg.norm(velocity_direction) + # velocity = velocity_magnitude * velocity_vector + + # # Set the solution step values for velocity components + # node.SetSolutionStepValue(KM.VELOCITY_X, velocity[0]) + # node.SetSolutionStepValue(KM.VELOCITY_Y, velocity[1]) + # node.SetSolutionStepValue(KM.VELOCITY_Z, velocity[2]) + + # # Fix the velocity components to ensure they remain constant during the simulation + # node.Fix(KM.VELOCITY_X) + # node.Fix(KM.VELOCITY_Y) + # node.Fix(KM.VELOCITY_Z) + + def ApplyRotationAndMeshDisplacement(self, model_part, axis, theta, center): + # Normalizing the rotation axis + axis = axis / np.sqrt(np.dot(axis, axis)) + a = np.cos(theta / 2) + b, c, d = -axis * np.sin(theta / 2) + + # Creating a quaternion rotation matrix + rot_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], + [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], + [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) + + for node in model_part.Nodes: + # Getting the initial coordinates of the node + point = np.array([node.X0, node.Y0, node.Z0]) + + # Shifting the rotation center to the origin + point -= center + + # Applying the rotation + rotated_point = np.dot(point, rot_matrix) + + # Shifting the point back and updating the node coordinates + rotated_point += center + node.X, node.Y, node.Z = rotated_point + node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_X, rotated_point[0] - node.X0) + node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Y, rotated_point[1] - node.Y0) + node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Z, rotated_point[2] - node.Z0) + node.Fix(KM.MESH_DISPLACEMENT_X) + node.Fix(KM.MESH_DISPLACEMENT_Y) + node.Fix(KM.MESH_DISPLACEMENT_Z) + + def set_node_velocities(self, model_part, axis, omega, center): + # Normalizing the rotation axis and creating the angular velocity vector + axis = axis / np.sqrt(np.dot(axis, axis)) + angular_velocity_vector = omega * axis + + for node in model_part.Nodes: + # Getting the current coordinates of the node + point = np.array([node.X, node.Y, node.Z]) + + # Calculating the position vector (relative to the rotation center) + position_vector = point - center + + # Computing the velocity due to rotation (v = omega cross r) + velocity_vector = np.cross(angular_velocity_vector, position_vector) + + # Setting the node's velocity + node.SetSolutionStepValue(KM.VELOCITY_X, velocity_vector[0]) + node.SetSolutionStepValue(KM.VELOCITY_Y, velocity_vector[1]) + node.SetSolutionStepValue(KM.VELOCITY_Z, velocity_vector[2]) + + # Fix the velocity components to ensure they remain constant during the simulation + node.Fix(KM.VELOCITY_X) + node.Fix(KM.VELOCITY_Y) + node.Fix(KM.VELOCITY_Z) + + # def ApplyRotationAndMeshDisplacement(self, rotating_frame_model_part, transformation_matrix): + # for node in rotating_frame_model_part.Nodes: + # coor = np.array([node.X0, node.Y0, node.Z0, 1.0]) # Convert point to homogeneous coordinates + # rotated_coor = transformation_matrix @ coor + # rotated_coor = rotated_coor[:3] / rotated_coor[3] # Convert back to 3D coordinates + # node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_X, rotated_coor[0] - node.X0) + # node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Y, rotated_coor[1] - node.Y0) + # node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Z, rotated_coor[2] - node.Z0) + # node.Fix(KM.MESH_DISPLACEMENT_X) + # node.Fix(KM.MESH_DISPLACEMENT_Y) + # node.Fix(KM.MESH_DISPLACEMENT_Z) + # node.X = rotated_coor[0] + # node.Y = rotated_coor[1] + # node.Z = rotated_coor[2] + + def get_theta_and_omega(self): + # Calculate the angular acceleration (alpha) + if np.isclose(self.acceleration_time, 0.0): + alpha = 0.0 + else: + alpha = self.angular_velocity_radians / self.acceleration_time + + # Calculate the angle (theta) based on current time + if self.time <= self.acceleration_time: + self.omega = alpha * self.time + self.theta = 0.5 * alpha * self.time**2 + else: + self.omega = self.angular_velocity_radians + self.theta = 0.5 * alpha * self.acceleration_time**2 + self.omega * (self.time - self.acceleration_time) + + def get_rotation_matrix(self, theta): + # Create a rotation matrix for the specified axis of rotation + axis = self.axis_of_rotation + c = np.cos(theta) + s = np.sin(theta) + C = 1 - c + x, y, z = axis + rotation_matrix = np.array([[x*x*C + c, x*y*C - z*s, x*z*C + y*s], + [y*x*C + z*s, y*y*C + c, y*z*C - x*s], + [z*x*C - y*s, z*y*C + x*s, z*z*C + c]]) + + # Create translation matrices for moving the center of rotation + translation_to_origin = np.eye(4) + translation_to_origin[:3, 3] = -self.center_of_rotation + translation_back = np.eye(4) + translation_back[:3, 3] = self.center_of_rotation + + # Create a homogeneous rotation matrix + homogeneous_rotation = np.eye(4) + homogeneous_rotation[:3, :3] = rotation_matrix + + # Combine the translation and rotation matrices to get the final transformation matrix + transformation_matrix = np.dot(np.dot(translation_back, homogeneous_rotation), translation_to_origin) + return transformation_matrix + + def ExecuteBeforeOutputStep(self): + """ This method is executed before writing the output (if output + is being written in this step) + + Keyword arguments: + self -- It signifies an instance of a class. + """ + pass + + def ExecuteAfterOutputStep(self): + """ This method is executed after writing the output (if output + is being written in this step) + + Keyword arguments: + self -- It signifies an instance of a class. + """ + pass + + def ExecuteFinalizeSolutionStep(self): + """ This method is executed in order to finalize the current step + + Keyword arguments: + self -- It signifies an instance of a class. + """ + pass + + def ExecuteFinalize(self): + """ This method is executed after the computations, at the end of the solution-loop + + Keyword arguments: + self -- It signifies an instance of a class. + """ + pass \ No newline at end of file From 856dc2860c11664866034c7f067ce7b6ad813b3c Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 18 May 2023 10:20:44 +0200 Subject: [PATCH 02/20] Working in serial (ApplyMeshDispAndRot). Cleaned. --- .../rotating_frame_utility.cpp | 35 ++-- .../custom_utilities/rotating_frame_utility.h | 20 ++- .../python_scripts/rotating_frame_process.py | 161 ++---------------- 3 files changed, 39 insertions(+), 177 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp index 6e82d1151b92..4d38bb27fe53 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -25,7 +25,7 @@ namespace Kratos void RotatingFrameUtility::ApplyVelocityToRotatingObject( ModelPart& rRotatingObjectModelPart, const double Omega, - const array_1d& rAxisOfRotation) + const array_1d& rAxisOfRotationOfRotation) { KRATOS_TRY; @@ -45,7 +45,7 @@ void RotatingFrameUtility::ApplyVelocityToRotatingObject( // Apply component by component the corresponding velocity. velocity_magnitude = Omega * radius; - velocity_direction = MathUtils::CrossProduct(rAxisOfRotation, rNode.Coordinates()); + velocity_direction = MathUtils::CrossProduct(rAxisOfRotationOfRotation, rNode.Coordinates()); velocity_vector = velocity_direction / norm_2(velocity_direction); velocity = velocity_magnitude * velocity_vector; @@ -65,17 +65,17 @@ void RotatingFrameUtility::ApplyVelocityToRotatingObject( void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( ModelPart& rRotatingFrameModelPart, - const array_1d& rAxis, + const array_1d& rAxisOfRotation, const double& rTheta, - const array_1d& rCenter) + const array_1d& rCenterOfRotation) { KRATOS_TRY; // Quaternion constants const double a = std::cos(rTheta / 2); - const double b = -rAxis[0] * std::sin(rTheta / 2); - const double c = -rAxis[1] * std::sin(rTheta / 2); - const double d = -rAxis[2] * std::sin(rTheta / 2); + const double b = -rAxisOfRotation[0] * std::sin(rTheta / 2); + const double c = -rAxisOfRotation[1] * std::sin(rTheta / 2); + const double d = -rAxisOfRotation[2] * std::sin(rTheta / 2); // Creating a quaternion rotation matrix BoundedMatrix rot_matrix; @@ -93,24 +93,21 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( block_for_each(rRotatingFrameModelPart.Nodes(), [&](Node& rNode) { // Getting the initial coordinates of the node auto& point = rNode.GetInitialPosition().Coordinates(); - // point[0] = node.X0(); - // point[0] = node.X0(); - // point[0] = node.X0(); // Shifting the rotation center to the origin - point -= rCenter; + auto centered_point = point - rCenterOfRotation; // Applying the rotation - array_1d rotated_point = prod(rot_matrix, point); + array_1d rotated_point = prod(centered_point, rot_matrix); // Shifting the point back and updating the node coordinates - rotated_point += rCenter; - rNode.FastGetSolutionStepValue(DISPLACEMENT_X) = rotated_point[0] - rNode.X0(); - rNode.FastGetSolutionStepValue(DISPLACEMENT_Y) = rotated_point[1] - rNode.Y0(); - rNode.FastGetSolutionStepValue(DISPLACEMENT_Z) = rotated_point[2] - rNode.Z0(); - rNode.Fix(DISPLACEMENT_X); - rNode.Fix(DISPLACEMENT_Y); - rNode.Fix(DISPLACEMENT_Z); + rotated_point += rCenterOfRotation; + rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_X) = rotated_point[0] - rNode.X0(); + rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Y) = rotated_point[1] - rNode.Y0(); + rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Z) = rotated_point[2] - rNode.Z0(); + rNode.Fix(MESH_DISPLACEMENT_X); + rNode.Fix(MESH_DISPLACEMENT_Y); + rNode.Fix(MESH_DISPLACEMENT_Z); rNode.X() = rotated_point[0]; rNode.Y() = rotated_point[1]; rNode.Z() = rotated_point[2]; diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h index 62d9b6011f7b..6b7fa5426b3e 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h @@ -83,22 +83,24 @@ class KRATOS_API(FLUID_DYNAMICS_APPLICATION) RotatingFrameUtility /** * @brief Apply rotation and mesh displacement to a rotating frame. - * This function applies rotation and mesh displacement to a rotating frame model part based on the provided rotation axis, angle and rotation center. - * It constructs a rotation matrix based on the rotation axis and angle. The rotation is applied to each node, which is then displaced + * This function applies rotation and mesh displacement to a rotating frame model + * part based on the provided rotation axis, angle and rotation center. + * It constructs a rotation matrix based on the rotation axis and angle. + * The rotation is applied to each node, which is then displaced * by the difference between the rotated coordinates and the initial coordinates. - * The solution step values for the mesh displacement components are set accordingly, and the mesh displacement components - * are fixed to remain constant during the simulation. Finally, the node coordinates are updated with the rotated - * coordinates. + * The solution step values for the mesh displacement components are set accordingly, + * and the mesh displacement components are fixed to remain constant during the simulation. + * Finally, the node coordinates are updated with the rotated coordinates. * @param rRotatingFrameModelPart The rotating frame model part to which the rotation and mesh displacement will be applied. - * @param rAxis The rotation axis vector. + * @param rAxisOfRotation The rotation axis vector. * @param rTheta The rotation angle. - * @param rCenter The center of rotation. + * @param rCenterOfRotation The center of rotation. */ static void ApplyRotationAndMeshDisplacement( ModelPart& rRotatingFrameModelPart, - const array_1d& rAxis, + const array_1d& rAxisOfRotation, const double& rTheta, - const array_1d& rCenter); + const array_1d& rCenterOfRotation); }; // Class RotatingFrameUtility diff --git a/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py index 507f651e25c2..19bf261ce81e 100644 --- a/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py +++ b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py @@ -38,7 +38,7 @@ def __init__(self, model, settings): "rotating_object_model_part_name": "", "center_of_rotation": [0.0,0.0,0.0], "axis_of_rotation": [1.0,0.0,0.0], - "angular_velocity_radians": 0.0, + "target_angular_velocity_radians": 0.0, "acceleration_time": 0.0 }""") @@ -84,10 +84,10 @@ def __init__(self, model, settings): raise Exception("The axis_of_rotation parameter is missing from the settings.") # Get the angular velocity in radians - if settings.Has("angular_velocity_radians"): - self.angular_velocity_radians = settings["angular_velocity_radians"].GetDouble() + if settings.Has("target_angular_velocity_radians"): + self.target_angular_velocity_radians = settings["target_angular_velocity_radians"].GetDouble() else: - raise Exception("The angular_velocity_radians parameter is missing from the settings.") + raise Exception("The target_angular_velocity_radians parameter is missing from the settings.") # Get the acceleration time if settings.Has("acceleration_time"): @@ -131,80 +131,17 @@ def ExecuteInitializeSolutionStep(self): self.time = self.rotating_frame_model_part.ProcessInfo[KM.TIME] # Calculate the angle (theta) for the specified time - self.get_theta_and_omega() #Potentially python - - # Calculate the rotation matrix. - # transformation_matrix = self.get_rotation_matrix(self.theta) #Potentially C++ + self.__GetThetaAndOmega() # Apply Rotation and Mesh Displacement to nodes. - # self.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, transformation_matrix) #Potentially C++ - # KratosCFD.RotatingFrameUtility.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, transformation_matrix) - self.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, self.axis_of_rotation, self.theta, self.center_of_rotation) + KratosCFD.RotatingFrameUtility.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, self.axis_of_rotation, self.theta, self.center_of_rotation) # Apply Velocity to object. - self.set_node_velocities(self.rotating_object_model_part, self.axis_of_rotation, self.omega, self.center_of_rotation) - # self.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.omega, self.center_of_rotation, self.axis_of_rotation) + self.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.axis_of_rotation, self.omega, self.center_of_rotation) # KratosCFD.RotatingFrameUtility.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.omega, self.axis_of_rotation) - def ApplyVelocityToRotatingObject(self, rotating_object_model_part, omega, center_of_rotation, axis_of_rotation): - for node in rotating_object_model_part.Nodes: - relative_coordinates = [node.X, node.Y, node.Z]-center_of_rotation - radius = np.linalg.norm(relative_coordinates) - #Apply component by component the corresponding velocity. - velocity_magnitude = omega * radius - velocity_direction = np.cross(axis_of_rotation, relative_coordinates) - velocity_vector = velocity_direction / np.linalg.norm(velocity_direction) - velocity = velocity_magnitude * velocity_vector - node.SetSolutionStepValue(KM.VELOCITY_X, velocity[0]) - node.SetSolutionStepValue(KM.VELOCITY_Y, velocity[1]) - node.SetSolutionStepValue(KM.VELOCITY_Z, velocity[2]) - node.Fix(KM.VELOCITY_X) - node.Fix(KM.VELOCITY_Y) - node.Fix(KM.VELOCITY_Z) - - def calculate_magnitude_difference(self, point, center_of_rotation, axis_of_rotation): - # Step 1: Calculate the vector from the target point to the given point - vector_difference = point - center_of_rotation - - # Step 2: Calculate the dot product between the vector difference and the target axis - dot_product = np.dot(vector_difference, axis_of_rotation) - - # Step 3: Calculate the magnitude of the target axis - magnitude_axis_of_rotation = np.linalg.norm(axis_of_rotation) - - # Step 4: Divide the dot product by the magnitude of the target axis to obtain the projection - projection = dot_product / magnitude_axis_of_rotation - - # Step 5: Calculate the magnitude of the projection to obtain the desired magnitude of difference - magnitude_difference = np.abs(projection) - - return magnitude_difference - - # def ApplyVelocityToRotatingObject(self, rotating_object_model_part, omega, center_of_rotation, axis_of_rotation): - # for node in rotating_object_model_part.Nodes: - # # Calculate the magnitude difference - # point = np.array([node.X, node.Y, node.Z]) - # magnitude_difference = self.calculate_magnitude_difference(point, center_of_rotation, axis_of_rotation) - - # # Apply component by component the corresponding velocity. - # velocity_magnitude = omega * magnitude_difference - # velocity_direction = np.cross(axis_of_rotation, point) - # velocity_vector = velocity_direction / np.linalg.norm(velocity_direction) - # velocity = velocity_magnitude * velocity_vector - - # # Set the solution step values for velocity components - # node.SetSolutionStepValue(KM.VELOCITY_X, velocity[0]) - # node.SetSolutionStepValue(KM.VELOCITY_Y, velocity[1]) - # node.SetSolutionStepValue(KM.VELOCITY_Z, velocity[2]) - - # # Fix the velocity components to ensure they remain constant during the simulation - # node.Fix(KM.VELOCITY_X) - # node.Fix(KM.VELOCITY_Y) - # node.Fix(KM.VELOCITY_Z) - def ApplyRotationAndMeshDisplacement(self, model_part, axis, theta, center): # Normalizing the rotation axis - axis = axis / np.sqrt(np.dot(axis, axis)) a = np.cos(theta / 2) b, c, d = -axis * np.sin(theta / 2) @@ -233,7 +170,7 @@ def ApplyRotationAndMeshDisplacement(self, model_part, axis, theta, center): node.Fix(KM.MESH_DISPLACEMENT_Y) node.Fix(KM.MESH_DISPLACEMENT_Z) - def set_node_velocities(self, model_part, axis, omega, center): + def ApplyVelocityToRotatingObject(self, model_part, axis, omega, center): # Normalizing the rotation axis and creating the angular velocity vector axis = axis / np.sqrt(np.dot(axis, axis)) angular_velocity_vector = omega * axis @@ -258,91 +195,17 @@ def set_node_velocities(self, model_part, axis, omega, center): node.Fix(KM.VELOCITY_Y) node.Fix(KM.VELOCITY_Z) - # def ApplyRotationAndMeshDisplacement(self, rotating_frame_model_part, transformation_matrix): - # for node in rotating_frame_model_part.Nodes: - # coor = np.array([node.X0, node.Y0, node.Z0, 1.0]) # Convert point to homogeneous coordinates - # rotated_coor = transformation_matrix @ coor - # rotated_coor = rotated_coor[:3] / rotated_coor[3] # Convert back to 3D coordinates - # node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_X, rotated_coor[0] - node.X0) - # node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Y, rotated_coor[1] - node.Y0) - # node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Z, rotated_coor[2] - node.Z0) - # node.Fix(KM.MESH_DISPLACEMENT_X) - # node.Fix(KM.MESH_DISPLACEMENT_Y) - # node.Fix(KM.MESH_DISPLACEMENT_Z) - # node.X = rotated_coor[0] - # node.Y = rotated_coor[1] - # node.Z = rotated_coor[2] - - def get_theta_and_omega(self): + def __GetThetaAndOmega(self): # Calculate the angular acceleration (alpha) if np.isclose(self.acceleration_time, 0.0): alpha = 0.0 else: - alpha = self.angular_velocity_radians / self.acceleration_time + alpha = self.target_angular_velocity_radians / self.acceleration_time # Calculate the angle (theta) based on current time if self.time <= self.acceleration_time: self.omega = alpha * self.time self.theta = 0.5 * alpha * self.time**2 else: - self.omega = self.angular_velocity_radians - self.theta = 0.5 * alpha * self.acceleration_time**2 + self.omega * (self.time - self.acceleration_time) - - def get_rotation_matrix(self, theta): - # Create a rotation matrix for the specified axis of rotation - axis = self.axis_of_rotation - c = np.cos(theta) - s = np.sin(theta) - C = 1 - c - x, y, z = axis - rotation_matrix = np.array([[x*x*C + c, x*y*C - z*s, x*z*C + y*s], - [y*x*C + z*s, y*y*C + c, y*z*C - x*s], - [z*x*C - y*s, z*y*C + x*s, z*z*C + c]]) - - # Create translation matrices for moving the center of rotation - translation_to_origin = np.eye(4) - translation_to_origin[:3, 3] = -self.center_of_rotation - translation_back = np.eye(4) - translation_back[:3, 3] = self.center_of_rotation - - # Create a homogeneous rotation matrix - homogeneous_rotation = np.eye(4) - homogeneous_rotation[:3, :3] = rotation_matrix - - # Combine the translation and rotation matrices to get the final transformation matrix - transformation_matrix = np.dot(np.dot(translation_back, homogeneous_rotation), translation_to_origin) - return transformation_matrix - - def ExecuteBeforeOutputStep(self): - """ This method is executed before writing the output (if output - is being written in this step) - - Keyword arguments: - self -- It signifies an instance of a class. - """ - pass - - def ExecuteAfterOutputStep(self): - """ This method is executed after writing the output (if output - is being written in this step) - - Keyword arguments: - self -- It signifies an instance of a class. - """ - pass - - def ExecuteFinalizeSolutionStep(self): - """ This method is executed in order to finalize the current step - - Keyword arguments: - self -- It signifies an instance of a class. - """ - pass - - def ExecuteFinalize(self): - """ This method is executed after the computations, at the end of the solution-loop - - Keyword arguments: - self -- It signifies an instance of a class. - """ - pass \ No newline at end of file + self.omega = self.target_angular_velocity_radians + self.theta = 0.5 * alpha * self.acceleration_time**2 + self.omega * (self.time - self.acceleration_time) \ No newline at end of file From 0e60412c0f7b99ca408bf9f47173d13d6c9f169b Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 18 May 2023 13:12:04 +0200 Subject: [PATCH 03/20] Working, only needs last checks, cleans, and tests. --- .../rotating_frame_utility.cpp | 112 ++++-------------- .../custom_utilities/rotating_frame_utility.h | 32 ++--- .../python_scripts/rotating_frame_process.py | 58 +-------- 3 files changed, 43 insertions(+), 159 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp index 4d38bb27fe53..01b67696c902 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -23,38 +23,36 @@ namespace Kratos { void RotatingFrameUtility::ApplyVelocityToRotatingObject( - ModelPart& rRotatingObjectModelPart, - const double Omega, - const array_1d& rAxisOfRotationOfRotation) + ModelPart& rModelPart, + const array_1d& rAxisOfRotation, + const double& rOmega, + const array_1d& rCenterOfRotation) { KRATOS_TRY; - // Get the array of nodes in the rotating object model part - auto& nodes_array = rRotatingObjectModelPart.Nodes(); + // Creating the angular velocity vector + array_1d angular_velocity_vector = rOmega * rAxisOfRotation; - // Declare thread-local variables outside the block_for_each - thread_local double radius; - thread_local double velocity_magnitude; - thread_local array_1d velocity_direction; + thread_local array_1d position_vector; thread_local array_1d velocity_vector; - thread_local array_1d velocity; // Apply the velocity calculations to each node in parallel - block_for_each(nodes_array, [&](Node& rNode) { - radius = std::sqrt(rNode.X() * rNode.X() + rNode.Y() * rNode.Y() + rNode.Z() * rNode.Z()); + block_for_each(rModelPart.Nodes(), [&](Node& rNode) { + // Getting the current coordinates of the node + auto& r_point = rNode.Coordinates(); + + // Calculating the position vector (relative to the rotation center) + position_vector = r_point - rCenterOfRotation; - // Apply component by component the corresponding velocity. - velocity_magnitude = Omega * radius; - velocity_direction = MathUtils::CrossProduct(rAxisOfRotationOfRotation, rNode.Coordinates()); - velocity_vector = velocity_direction / norm_2(velocity_direction); - velocity = velocity_magnitude * velocity_vector; + // Computing the velocity due to rotation (v = omega cross r) + velocity_vector = MathUtils::CrossProduct(angular_velocity_vector, position_vector); - // Set the solution step values for the velocity components using FastGetSolutionStepValue - rNode.FastGetSolutionStepValue(VELOCITY_X) = velocity[0]; - rNode.FastGetSolutionStepValue(VELOCITY_Y) = velocity[1]; - rNode.FastGetSolutionStepValue(VELOCITY_Z) = velocity[2]; + // Setting the node's velocity + rNode.FastGetSolutionStepValue(VELOCITY_X) = velocity_vector[0]; + rNode.FastGetSolutionStepValue(VELOCITY_Y) = velocity_vector[1]; + rNode.FastGetSolutionStepValue(VELOCITY_Z) = velocity_vector[2]; - // Fix the velocity components to ensure they remain constant during the simulation + // Fix the velocity components rNode.Fix(VELOCITY_X); rNode.Fix(VELOCITY_Y); rNode.Fix(VELOCITY_Z); @@ -88,17 +86,19 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( rot_matrix(2, 0) = 2*(b*d-a*c); rot_matrix(2, 1) = 2*(c*d+a*b); rot_matrix(2, 2) = a*a+d*d-b*b-c*c; + + thread_local array_1d rotated_point; // Apply the rotation and mesh displacement calculations to each node in parallel block_for_each(rRotatingFrameModelPart.Nodes(), [&](Node& rNode) { // Getting the initial coordinates of the node - auto& point = rNode.GetInitialPosition().Coordinates(); + auto& r_point = rNode.GetInitialPosition().Coordinates(); // Shifting the rotation center to the origin - auto centered_point = point - rCenterOfRotation; + auto centered_point = r_point - rCenterOfRotation; // Applying the rotation - array_1d rotated_point = prod(centered_point, rot_matrix); + rotated_point = prod(centered_point, rot_matrix); // Shifting the point back and updating the node coordinates rotated_point += rCenterOfRotation; @@ -116,67 +116,5 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( KRATOS_CATCH(""); } - -// void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( -// ModelPart& rRotatingFrameModelPart, -// const Matrix& rTransformationMatrix) -// { -// KRATOS_TRY; - -// // Get the array of nodes in the rotating frame model part -// auto& nodes_array = rRotatingFrameModelPart.Nodes(); - -// // Declare thread-local variables outside the block_for_each -// thread_local array_1d homogeneous_coordinates; -// thread_local array_1d rotated_coordinates; -// thread_local array_1d rotated_coor; -// thread_local double mesh_displacement_x; -// thread_local double mesh_displacement_y; -// thread_local double mesh_displacement_z; - -// // Apply the rotation and mesh displacement calculations to each node in parallel -// block_for_each(nodes_array, [&](Node& rNode) { - -// // Convert point to homogeneous coordinates -// // homogeneous_coordinates.resize(4); -// homogeneous_coordinates[0] = rNode.X0(); -// homogeneous_coordinates[1] = rNode.Y0(); -// homogeneous_coordinates[2] = rNode.Z0(); -// homogeneous_coordinates[3] = 1.0; - -// // Apply transformation to homogeneous coordinates -// // rotated_coordinates.resize(4); -// noalias(rotated_coordinates) = prod(rTransformationMatrix, homogeneous_coordinates); - -// // Convert back to 3D coordinates -// // rotated_coor.resize(3); -// rotated_coor[0] = rotated_coordinates[0] / rotated_coordinates[3]; -// rotated_coor[1] = rotated_coordinates[1] / rotated_coordinates[3]; -// rotated_coor[2] = rotated_coordinates[2] / rotated_coordinates[3]; - -// // Calculate mesh displacement -// mesh_displacement_x = rotated_coor[0] - rNode.X0(); -// mesh_displacement_y = rotated_coor[1] - rNode.Y0(); -// mesh_displacement_z = rotated_coor[2] - rNode.Z0(); - -// // Set the solution step values for mesh displacement using FastGetSolutionStepValue -// rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_X) = mesh_displacement_x; -// rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Y) = mesh_displacement_y; -// rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Z) = mesh_displacement_z; - -// // Fix the mesh displacement components to ensure they remain constant during the simulation -// rNode.Fix(MESH_DISPLACEMENT_X); -// rNode.Fix(MESH_DISPLACEMENT_Y); -// rNode.Fix(MESH_DISPLACEMENT_Z); - -// // Update the coordinates of the node using FastGetSolutionStepValue -// rNode.X() = rotated_coor[0]; -// rNode.Y() = rotated_coor[1]; -// rNode.Z() = rotated_coor[2]; -// }); - -// KRATOS_CATCH(""); -// } - } // namespace Kratos. diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h index 6b7fa5426b3e..7892b34c32f9 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h @@ -62,24 +62,26 @@ class KRATOS_API(FLUID_DYNAMICS_APPLICATION) RotatingFrameUtility ///@{ /** - * @brief Apply velocity to rotating object. - * - * This function applies velocity to a rotating object model part. - * The velocity is computed based on the angular velocity and the - * axis of rotation. The function calculates the radius of each node, - * and then applies the corresponding velocity component-wise. The - * velocity components are fixed to remain constant during the simulation, - * and the solution step values for the velocity components are set - * accordingly. - * - * @param rRotatingObjectModelPart The rotating object model part. - * @param Omega The angular velocity. + * @brief Applies an angular velocity to nodes in a rotating frame. + * + * This function applies an angular velocity to each node in a provided model part. + * The angular velocity vector is calculated from a provided rotation axis and rotation speed. + * The position vector for each node is determined relative to a provided rotation center. + * The velocity due to rotation is then computed as the cross product of the angular velocity vector and the position vector. + * The calculated velocity is set as the solution step value for each node. + * Finally, the function fixes the velocity components to ensure they remain constant during the simulation. + * + * @param rModelPart The model part in which to apply the velocity. * @param rAxisOfRotation The axis of rotation. + * @param rOmega The rotation speed. + * @param rCenterOfRotation The center of rotation. */ static void ApplyVelocityToRotatingObject( - ModelPart& rRotatingObjectModelPart, - const double Omega, - const array_1d& rAxisOfRotation); + ModelPart& rModelPart, + const array_1d& rAxisOfRotation, + const double& rOmega, + const array_1d& rCenterOfRotation); + /** * @brief Apply rotation and mesh displacement to a rotating frame. diff --git a/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py index 19bf261ce81e..40bdbcebf7ec 100644 --- a/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py +++ b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py @@ -137,63 +137,7 @@ def ExecuteInitializeSolutionStep(self): KratosCFD.RotatingFrameUtility.ApplyRotationAndMeshDisplacement(self.rotating_frame_model_part, self.axis_of_rotation, self.theta, self.center_of_rotation) # Apply Velocity to object. - self.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.axis_of_rotation, self.omega, self.center_of_rotation) - # KratosCFD.RotatingFrameUtility.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.omega, self.axis_of_rotation) - - def ApplyRotationAndMeshDisplacement(self, model_part, axis, theta, center): - # Normalizing the rotation axis - a = np.cos(theta / 2) - b, c, d = -axis * np.sin(theta / 2) - - # Creating a quaternion rotation matrix - rot_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], - [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], - [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) - - for node in model_part.Nodes: - # Getting the initial coordinates of the node - point = np.array([node.X0, node.Y0, node.Z0]) - - # Shifting the rotation center to the origin - point -= center - - # Applying the rotation - rotated_point = np.dot(point, rot_matrix) - - # Shifting the point back and updating the node coordinates - rotated_point += center - node.X, node.Y, node.Z = rotated_point - node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_X, rotated_point[0] - node.X0) - node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Y, rotated_point[1] - node.Y0) - node.SetSolutionStepValue(KM.MESH_DISPLACEMENT_Z, rotated_point[2] - node.Z0) - node.Fix(KM.MESH_DISPLACEMENT_X) - node.Fix(KM.MESH_DISPLACEMENT_Y) - node.Fix(KM.MESH_DISPLACEMENT_Z) - - def ApplyVelocityToRotatingObject(self, model_part, axis, omega, center): - # Normalizing the rotation axis and creating the angular velocity vector - axis = axis / np.sqrt(np.dot(axis, axis)) - angular_velocity_vector = omega * axis - - for node in model_part.Nodes: - # Getting the current coordinates of the node - point = np.array([node.X, node.Y, node.Z]) - - # Calculating the position vector (relative to the rotation center) - position_vector = point - center - - # Computing the velocity due to rotation (v = omega cross r) - velocity_vector = np.cross(angular_velocity_vector, position_vector) - - # Setting the node's velocity - node.SetSolutionStepValue(KM.VELOCITY_X, velocity_vector[0]) - node.SetSolutionStepValue(KM.VELOCITY_Y, velocity_vector[1]) - node.SetSolutionStepValue(KM.VELOCITY_Z, velocity_vector[2]) - - # Fix the velocity components to ensure they remain constant during the simulation - node.Fix(KM.VELOCITY_X) - node.Fix(KM.VELOCITY_Y) - node.Fix(KM.VELOCITY_Z) + KratosCFD.RotatingFrameUtility.ApplyVelocityToRotatingObject(self.rotating_object_model_part, self.axis_of_rotation, self.omega, self.center_of_rotation) def __GetThetaAndOmega(self): # Calculate the angular acceleration (alpha) From ba331e2729286a4cf01b13f8a646689853c65371 Mon Sep 17 00:00:00 2001 From: SADPR Date: Fri, 19 May 2023 10:56:42 +0200 Subject: [PATCH 04/20] Cleaning Python process --- .../python_scripts/rotating_frame_process.py | 24 ------------------- 1 file changed, 24 deletions(-) diff --git a/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py index 40bdbcebf7ec..6cbcdf356d16 100644 --- a/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py +++ b/applications/FluidDynamicsApplication/python_scripts/rotating_frame_process.py @@ -97,30 +97,6 @@ def __init__(self, model, settings): else: raise Exception("The acceleration_time parameter is missing from the settings.") - def ExecuteInitialize(self): - """ This method is executed at the begining to initialize the process - - Keyword arguments: - self -- It signifies an instance of a class. - """ - pass - - def Check(self): - """ This method verifies that the input is correct - - Keyword arguments: - self -- It signifies an instance of a class. - """ - pass - - def ExecuteBeforeSolutionLoop(self): - """ This method is executed just before the solution-loop - - Keyword arguments: - self -- It signifies an instance of a class. - """ - pass - def ExecuteInitializeSolutionStep(self): """ This method is executed in order to initialize the current step From 91164223ea1605f8b0fa3479abfa2e6044d44dbf Mon Sep 17 00:00:00 2001 From: SADPR Date: Fri, 19 May 2023 10:57:16 +0200 Subject: [PATCH 05/20] Adding MDPA. Maybe clean. --- .../RotatingFrameProcessTest.mdpa | 357 ++++++++++++++++++ 1 file changed, 357 insertions(+) create mode 100644 applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa diff --git a/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa b/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa new file mode 100644 index 000000000000..c96915c90128 --- /dev/null +++ b/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa @@ -0,0 +1,357 @@ +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 Conditions LineCondition2D2N// GUI group identifier: _HIDDEN__SKIN_ + 1 0 40 38 + 2 0 38 37 + 3 0 37 35 + 4 0 35 31 + 5 0 31 29 + 6 0 29 24 + 7 0 24 19 + 8 0 19 13 + 9 0 13 11 + 10 0 11 8 + 11 0 8 4 + 12 0 4 1 + 13 0 1 3 + 14 0 3 7 + 15 0 7 12 + 16 0 12 14 + 17 0 14 20 + 18 0 20 23 + 19 0 23 28 + 20 0 28 30 + 21 0 30 34 + 22 0 34 36 + 23 0 36 39 + 24 0 39 40 +End Conditions + +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_Pressure // Group Pressure // Subtree GENERIC + Begin SubModelPartNodes + 9 + 10 + 26 + 27 + End SubModelPartNodes + Begin SubModelPartElements + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart +Begin SubModelPart FluidParts_RotatingFrame // Group RotatingFrame // Subtree FluidParts + 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 +Begin SubModelPart VELOCITY_RotatingFrame // Group RotatingFrame // Subtree VELOCITY + 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 + End SubModelPartElements + Begin SubModelPartConditions + End SubModelPartConditions +End SubModelPart +Begin SubModelPart NoSlip2D_NoSlip // Group NoSlip // Subtree NoSlip2D + Begin SubModelPartNodes + 1 + 3 + 4 + 7 + 8 + 11 + 12 + 13 + 14 + 19 + 20 + 23 + 24 + 28 + 29 + 30 + 31 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + End SubModelPartNodes + Begin SubModelPartElements + End SubModelPartElements + Begin SubModelPartConditions + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + End SubModelPartConditions +End SubModelPart From 765b89dfce20d956403503c706ea373e48bdab71 Mon Sep 17 00:00:00 2001 From: SADPR Date: Fri, 19 May 2023 10:57:34 +0200 Subject: [PATCH 06/20] Results file (JSON). Test. --- .../RotatingFrameProcessTestResults.json | 805 ++++++++++++++++++ 1 file changed, 805 insertions(+) create mode 100644 applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTestResults.json diff --git a/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTestResults.json b/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTestResults.json new file mode 100644 index 000000000000..c6b77fbf3b1a --- /dev/null +++ b/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTestResults.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 From 6d5498876b4f2fa7ff270e7f85091f2ea5887f0e Mon Sep 17 00:00:00 2001 From: SADPR Date: Fri, 19 May 2023 10:57:57 +0200 Subject: [PATCH 07/20] Added test. --- .../tests/test_rotating_frame_process.py | 84 +++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py diff --git a/applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py b/applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py new file mode 100644 index 000000000000..8058171761d6 --- /dev/null +++ b/applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py @@ -0,0 +1,84 @@ +# Kratos Imports +import KratosMultiphysics +import KratosMultiphysics.KratosUnittest as KratosUnittest +import os + +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") + +@KratosUnittest.skipUnless(mesh_moving_is_available,"MeshMovingApplication is not available") +def GetFilePath(fileName): + return os.path.join(os.path.dirname(os.path.realpath(__file__)), fileName) + +def GenerateModel(): + # Set up the model + 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("RotatingFrameProcessTest/RotatingFrameProcessTest")) + 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 + model_part.ProcessInfo[KratosMultiphysics.DELTA_TIME] = 1 + + return model, model_part + + +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.FluidParts_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 + }""") + + # Create and execute process + from KratosMultiphysics.FluidDynamicsApplication.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("RotatingFrameProcessTest/RotatingFrameProcessTestResults.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("RotatingFrameProcessTest/RotatingFrameProcessTestResults.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 1e59954f454a89869c42319130fd28e0ac138646 Mon Sep 17 00:00:00 2001 From: SADPR Date: Fri, 19 May 2023 11:51:11 +0200 Subject: [PATCH 08/20] Cleaned MDPA. --- .../RotatingFrameProcessTest.mdpa | 145 +----------------- 1 file changed, 1 insertion(+), 144 deletions(-) diff --git a/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa b/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa index c96915c90128..e82556723203 100644 --- a/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa +++ b/applications/FluidDynamicsApplication/tests/RotatingFrameProcessTest/RotatingFrameProcessTest.mdpa @@ -100,33 +100,6 @@ Begin Elements Element2D3N// GUI group identifier: RotatingFrame 48 0 7 10 2 End Elements -Begin Conditions LineCondition2D2N// GUI group identifier: _HIDDEN__SKIN_ - 1 0 40 38 - 2 0 38 37 - 3 0 37 35 - 4 0 35 31 - 5 0 31 29 - 6 0 29 24 - 7 0 24 19 - 8 0 19 13 - 9 0 13 11 - 10 0 11 8 - 11 0 8 4 - 12 0 4 1 - 13 0 1 3 - 14 0 3 7 - 15 0 7 12 - 16 0 12 14 - 17 0 14 20 - 18 0 20 23 - 19 0 23 28 - 20 0 28 30 - 21 0 30 34 - 22 0 34 36 - 23 0 36 39 - 24 0 39 40 -End Conditions - Begin SubModelPart GENERIC_RotatingObject // Group RotatingObject // Subtree GENERIC Begin SubModelPartNodes 6 @@ -143,19 +116,7 @@ Begin SubModelPart GENERIC_RotatingObject // Group RotatingObject // Subtree GEN Begin SubModelPartConditions End SubModelPartConditions End SubModelPart -Begin SubModelPart GENERIC_Pressure // Group Pressure // Subtree GENERIC - Begin SubModelPartNodes - 9 - 10 - 26 - 27 - End SubModelPartNodes - Begin SubModelPartElements - End SubModelPartElements - Begin SubModelPartConditions - End SubModelPartConditions -End SubModelPart -Begin SubModelPart FluidParts_RotatingFrame // Group RotatingFrame // Subtree FluidParts +Begin SubModelPart GENERIC_RotatingFrame // Group RotatingFrame // Subtree GENERIC Begin SubModelPartNodes 1 2 @@ -251,107 +212,3 @@ Begin SubModelPart FluidParts_RotatingFrame // Group RotatingFrame // Subtree Fl Begin SubModelPartConditions End SubModelPartConditions End SubModelPart -Begin SubModelPart VELOCITY_RotatingFrame // Group RotatingFrame // Subtree VELOCITY - 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 - End SubModelPartElements - Begin SubModelPartConditions - End SubModelPartConditions -End SubModelPart -Begin SubModelPart NoSlip2D_NoSlip // Group NoSlip // Subtree NoSlip2D - Begin SubModelPartNodes - 1 - 3 - 4 - 7 - 8 - 11 - 12 - 13 - 14 - 19 - 20 - 23 - 24 - 28 - 29 - 30 - 31 - 34 - 35 - 36 - 37 - 38 - 39 - 40 - End SubModelPartNodes - Begin SubModelPartElements - End SubModelPartElements - Begin SubModelPartConditions - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9 - 10 - 11 - 12 - 13 - 14 - 15 - 16 - 17 - 18 - 19 - 20 - 21 - 22 - 23 - 24 - End SubModelPartConditions -End SubModelPart From 1c4ebe5be1c5b2da5f2774ad6f8d81a2a5c8e8a0 Mon Sep 17 00:00:00 2001 From: SADPR Date: Fri, 19 May 2023 11:51:40 +0200 Subject: [PATCH 09/20] Changing name of Model Part --- .../tests/test_rotating_frame_process.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py b/applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py index 8058171761d6..d150dd5bd971 100644 --- a/applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py +++ b/applications/FluidDynamicsApplication/tests/test_rotating_frame_process.py @@ -41,7 +41,7 @@ def test_rotating_frame_process(self): parameters = KratosMultiphysics.Parameters(""" { - "rotating_frame_model_part_name": "Main.FluidParts_RotatingFrame", + "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], From 4bf94c2fc52ccfb501f2d9894ec6c63cb4cd20f9 Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 22 Aug 2023 11:01:06 +0200 Subject: [PATCH 10/20] Assign values using the whole vector --- .../custom_utilities/rotating_frame_utility.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp index 01b67696c902..e4a22b689588 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -48,9 +48,10 @@ void RotatingFrameUtility::ApplyVelocityToRotatingObject( velocity_vector = MathUtils::CrossProduct(angular_velocity_vector, position_vector); // Setting the node's velocity - rNode.FastGetSolutionStepValue(VELOCITY_X) = velocity_vector[0]; - rNode.FastGetSolutionStepValue(VELOCITY_Y) = velocity_vector[1]; - rNode.FastGetSolutionStepValue(VELOCITY_Z) = velocity_vector[2]; + // rNode.FastGetSolutionStepValue(VELOCITY_X) = velocity_vector[0]; + // rNode.FastGetSolutionStepValue(VELOCITY_Y) = velocity_vector[1]; + // rNode.FastGetSolutionStepValue(VELOCITY_Z) = velocity_vector[2]; + rNode.FastGetSolutionStepValue(VELOCITY) = velocity_vector; // Fix the velocity components rNode.Fix(VELOCITY_X); From af0c22381b166a5fb416de8c37ced103be16862c Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 22 Aug 2023 11:14:53 +0200 Subject: [PATCH 11/20] Assign vector --- .../custom_utilities/rotating_frame_utility.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp index e4a22b689588..363c9811fcd0 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -89,6 +89,7 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( rot_matrix(2, 2) = a*a+d*d-b*b-c*c; thread_local array_1d rotated_point; + thread_local array_1d mesh_displacement; // Apply the rotation and mesh displacement calculations to each node in parallel block_for_each(rRotatingFrameModelPart.Nodes(), [&](Node& rNode) { @@ -103,9 +104,14 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( // Shifting the point back and updating the node coordinates rotated_point += rCenterOfRotation; - rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_X) = rotated_point[0] - rNode.X0(); - rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Y) = rotated_point[1] - rNode.Y0(); - rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Z) = rotated_point[2] - rNode.Z0(); + mesh_displacement[0] = rotated_point[0] - rNode.X0(); + mesh_displacement[1] = rotated_point[1] - rNode.Y0(); + mesh_displacement[2] = rotated_point[2] - rNode.Z0(); + + rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT) = mesh_displacement; + // rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_X) = rotated_point[0] - rNode.X0(); + // rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Y) = rotated_point[1] - rNode.Y0(); + // rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Z) = rotated_point[2] - rNode.Z0(); rNode.Fix(MESH_DISPLACEMENT_X); rNode.Fix(MESH_DISPLACEMENT_Y); rNode.Fix(MESH_DISPLACEMENT_Z); From 870a4fcb4b4dd289f0bb4e46b7df988622ca4565 Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 22 Aug 2023 11:19:51 +0200 Subject: [PATCH 12/20] Assign by vector (coordinates). --- .../custom_utilities/rotating_frame_utility.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp index 363c9811fcd0..7aa3b4035567 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -115,9 +115,10 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( rNode.Fix(MESH_DISPLACEMENT_X); rNode.Fix(MESH_DISPLACEMENT_Y); rNode.Fix(MESH_DISPLACEMENT_Z); - rNode.X() = rotated_point[0]; - rNode.Y() = rotated_point[1]; - rNode.Z() = rotated_point[2]; + rNode.Coordinates() = rotated_point; + // rNode.X() = rotated_point[0]; + // rNode.Y() = rotated_point[1]; + // rNode.Z() = rotated_point[2]; }); KRATOS_CATCH(""); From cec81ab5a181ecbe4929f8f582c323a7fbdc6da8 Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 22 Aug 2023 11:23:13 +0200 Subject: [PATCH 13/20] Pragma once --- .../custom_utilities/rotating_frame_utility.h | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h index 7892b34c32f9..e0e44c755f15 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h @@ -10,8 +10,7 @@ // Main author: Sebastian Ares de Parga Regalado // -#if !defined(KRATOS_ROTATING_FRAME_UTILITY_H) -#define KRATOS_ROTATING_FRAME_UTILITY_H +#pragma once // System includes #include @@ -108,6 +107,4 @@ class KRATOS_API(FLUID_DYNAMICS_APPLICATION) RotatingFrameUtility ///@} -} // namespace Kratos. - -#endif // KRATOS_ROTATING_FRAME_UTILITY_H \ No newline at end of file +} // namespace Kratos. \ No newline at end of file From 68207b052d37b9ff6fec8196b1d737ddd14f090e Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 22 Aug 2023 11:38:02 +0200 Subject: [PATCH 14/20] Last cleaning. --- .../custom_utilities/rotating_frame_utility.cpp | 9 --------- 1 file changed, 9 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp index 7aa3b4035567..d18a15cf0a83 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -48,9 +48,6 @@ void RotatingFrameUtility::ApplyVelocityToRotatingObject( velocity_vector = MathUtils::CrossProduct(angular_velocity_vector, position_vector); // Setting the node's velocity - // rNode.FastGetSolutionStepValue(VELOCITY_X) = velocity_vector[0]; - // rNode.FastGetSolutionStepValue(VELOCITY_Y) = velocity_vector[1]; - // rNode.FastGetSolutionStepValue(VELOCITY_Z) = velocity_vector[2]; rNode.FastGetSolutionStepValue(VELOCITY) = velocity_vector; // Fix the velocity components @@ -109,16 +106,10 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( mesh_displacement[2] = rotated_point[2] - rNode.Z0(); rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT) = mesh_displacement; - // rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_X) = rotated_point[0] - rNode.X0(); - // rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Y) = rotated_point[1] - rNode.Y0(); - // rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT_Z) = rotated_point[2] - rNode.Z0(); rNode.Fix(MESH_DISPLACEMENT_X); rNode.Fix(MESH_DISPLACEMENT_Y); rNode.Fix(MESH_DISPLACEMENT_Z); rNode.Coordinates() = rotated_point; - // rNode.X() = rotated_point[0]; - // rNode.Y() = rotated_point[1]; - // rNode.Z() = rotated_point[2]; }); KRATOS_CATCH(""); From 3ccb49ef6a43eab473cb41a4b3ec0d74e0c8480d Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 22 Aug 2023 11:51:01 +0200 Subject: [PATCH 15/20] Riccardo's suggestion (from another discussion). --- .../custom_utilities/rotating_frame_utility.cpp | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp index d18a15cf0a83..f8b94c922330 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.cpp @@ -33,19 +33,16 @@ void RotatingFrameUtility::ApplyVelocityToRotatingObject( // Creating the angular velocity vector array_1d angular_velocity_vector = rOmega * rAxisOfRotation; - thread_local array_1d position_vector; - thread_local array_1d velocity_vector; - // Apply the velocity calculations to each node in parallel block_for_each(rModelPart.Nodes(), [&](Node& rNode) { // Getting the current coordinates of the node auto& r_point = rNode.Coordinates(); // Calculating the position vector (relative to the rotation center) - position_vector = r_point - rCenterOfRotation; + array_1d position_vector = r_point - rCenterOfRotation; // Computing the velocity due to rotation (v = omega cross r) - velocity_vector = MathUtils::CrossProduct(angular_velocity_vector, position_vector); + array_1d velocity_vector = MathUtils::CrossProduct(angular_velocity_vector, position_vector); // Setting the node's velocity rNode.FastGetSolutionStepValue(VELOCITY) = velocity_vector; @@ -84,9 +81,6 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( rot_matrix(2, 0) = 2*(b*d-a*c); rot_matrix(2, 1) = 2*(c*d+a*b); rot_matrix(2, 2) = a*a+d*d-b*b-c*c; - - thread_local array_1d rotated_point; - thread_local array_1d mesh_displacement; // Apply the rotation and mesh displacement calculations to each node in parallel block_for_each(rRotatingFrameModelPart.Nodes(), [&](Node& rNode) { @@ -97,15 +91,12 @@ void RotatingFrameUtility::ApplyRotationAndMeshDisplacement( auto centered_point = r_point - rCenterOfRotation; // Applying the rotation - rotated_point = prod(centered_point, rot_matrix); + array_1d rotated_point = prod(centered_point, rot_matrix); // Shifting the point back and updating the node coordinates rotated_point += rCenterOfRotation; - mesh_displacement[0] = rotated_point[0] - rNode.X0(); - mesh_displacement[1] = rotated_point[1] - rNode.Y0(); - mesh_displacement[2] = rotated_point[2] - rNode.Z0(); - rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT) = mesh_displacement; + rNode.FastGetSolutionStepValue(MESH_DISPLACEMENT) = rotated_point - rNode.GetInitialPosition(); rNode.Fix(MESH_DISPLACEMENT_X); rNode.Fix(MESH_DISPLACEMENT_Y); rNode.Fix(MESH_DISPLACEMENT_Z); From 94abf9a0581f3bda3a06f851af31c85436cb6aed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Thu, 19 Oct 2023 18:37:12 +0200 Subject: [PATCH 16/20] Blocking again --- .../custom_utilities/rotating_frame_utility.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h index e0e44c755f15..556f5ded5dac 100644 --- a/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h +++ b/applications/FluidDynamicsApplication/custom_utilities/rotating_frame_utility.h @@ -29,7 +29,8 @@ #include "utilities/builtin_timer.h" #include "utilities/math_utils.h" -namespace Kratos { +namespace Kratos +{ ///@addtogroup FluidDynamicsApplication ///@{ From 54917778e22f21b619b1c8a932ed2ae47af62082 Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 5 Dec 2023 15:45:41 +0100 Subject: [PATCH 17/20] Erase this, is not related to this.!!!!!! I will continue this PR after. --- ...lave_constraints_to_neighbours_utility.cpp | 68 +++++++++++++++++-- ..._slave_constraints_to_neighbours_utility.h | 9 ++- 2 files changed, 70 insertions(+), 7 deletions(-) diff --git a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp index 0cda733a60bf..dbfd9b00f0c1 100644 --- a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp +++ b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp @@ -20,7 +20,7 @@ // Project includes #include "utilities/assign_master_slave_constraints_to_neighbours_utility.h" -namespace Kratos +namespace Kratos { ///@addtogroup Kratos Core ///@{ @@ -43,8 +43,8 @@ AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstraintsToN NodesContainerType::ContainerType& nodes_model_part = rMasterStructureNodes.GetContainer(); mpBins = Kratos::make_unique(nodes_model_part.begin(), nodes_model_part.end()); mMaxNumberOfNodes = rMasterStructureNodes.size(); - KRATOS_CATCH(""); -} + KRATOS_CATCH(""); +} // Destructor AssignMasterSlaveConstraintsToNeighboursUtility::~AssignMasterSlaveConstraintsToNeighboursUtility() {} @@ -97,7 +97,7 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::GetDofsAndCoordinatesForSl } else { std::stringstream variables_str; variables_str << "["; - for (const auto& variable : rVariableList) + for (const auto& variable : rVariableList) variables_str << variable.get() << ", "; variables_str.seekp(-2, std::ios_base::end); // Remove the last ", " variables_str << "]"; @@ -206,7 +206,7 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstrain // Get Dofs and Coordinates for the slave node and the cloud of nodes GetDofsAndCoordinatesForSlaveNode(p_slave_node, rVariableList, r_slave_dofs, r_slave_coordinates); GetDofsAndCoordinatesForCloudOfNodes(r_cloud_of_nodes, rVariableList, r_cloud_of_dofs, r_cloud_of_nodes_coordinates); - + // Calculate shape functions RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_slave_coordinates, r_n_container); @@ -235,5 +235,61 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstrain KRATOS_INFO("AssignMasterSlaveConstraintsToNeighboursUtility") << "Build and Assign Master-Slave Constraints Time: " << build_and_assign_mscs.ElapsedSeconds() << std::endl; KRATOS_CATCH(""); } -} // namespace Kratos. +void AssignMasterSlaveConstraintsToNeighboursUtility::InterpolateVelocityToNodes( + NodesContainerType pDestinationNodes, + double const Radius, + const std::vector>>& rVariableList, + double const MinNumOfNeighNodes + ) +{ + KRATOS_TRY; + + BuiltinTimer build_and_assign_mscs; + + // Declare thread-local storage (TLS) variables + thread_local std::vector r_cloud_of_dofs; // The DOFs of the cloud of nodes + thread_local std::vector r_destination_dofs; // The DOFs of the destination node + thread_local Matrix r_cloud_of_nodes_coordinates; // The coordinates of the cloud of nodes + thread_local array_1d r_destination_coordinates; // The coordinates of the destination node + thread_local Vector r_n_container; // Container for shape functions + + auto& p_nodes_array = pDestinationNodes.GetContainer(); // Get the container of destination nodes + + // Declare a counter variable outside the loop as std::atomic + std::atomic i(0); + + block_for_each(p_nodes_array, [&](Node::Pointer& p_destination_node) { + + double local_radius = Radius; + + ResultNodesContainerType r_cloud_of_nodes(0); // Container for the cloud of nodes + ResultNodesContainerType r_local_results(mMaxNumberOfNodes); // Container for node search results + + while (r_cloud_of_nodes.size() < MinNumOfNeighNodes) { + r_cloud_of_nodes.clear(); + SearchNodesInRadiusForNode(p_destination_node, local_radius, r_cloud_of_nodes, r_local_results); // Search for nodes in the radius + local_radius += Radius; + } + + // Get Dofs and Coordinates for the slave node and the cloud of nodes + GetDofsAndCoordinatesForSlaveNode(p_destination_node, rVariableList, r_destination_dofs, r_destination_coordinates); + GetDofsAndCoordinatesForCloudOfNodes(r_cloud_of_nodes, rVariableList, r_cloud_of_dofs, r_cloud_of_nodes_coordinates); + + // Calculate shape functions + RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_destination_coordinates, r_n_container); + + // IndexType it = i.fetch_add(1) * rVariableList.size(); // Atomically increment the counter and get the previous value + + for (std::size_t j = 0; j < rVariableList.size(); ++j) { + double interpolated_variable_value = 0.0; + for (std::size_t k = 0; k < r_cloud_of_nodes.size(); ++k){ + auto r_current_node = r_cloud_of_nodes[k]; + interpolated_variable_value += r_current_node->FastGetSolutionStepValue(rVariableList[j].get()) * r_n_container[k]; + } + p_destination_node->FastGetSolutionStepValue(rVariableList[j].get(), interpolated_variable_value); + } + }); + KRATOS_CATCH(""); +} // namespace Kratos. +} \ No newline at end of file diff --git a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h index ec88495e7ed5..e4c66cbd225b 100644 --- a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h +++ b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h @@ -57,7 +57,7 @@ class KRATOS_API(KRATOS_CORE) AssignMasterSlaveConstraintsToNeighboursUtility public: ///@name Type Definitions ///@{ - + //Node Types using NodeType = ModelPart::NodeType; @@ -155,6 +155,13 @@ class KRATOS_API(KRATOS_CORE) AssignMasterSlaveConstraintsToNeighboursUtility double const MinNumOfNeighNodes ); + void InterpolateVelocityToNodes( + NodesContainerType pDestinationNodes, + double const Radius, + const std::vector>>& rVariableList, + double const MinNumOfNeighNodes + ); + ///@} ///@name Input and output ///@{ From 5be37fedc1ab27dadab5387178a6f2d58f2e2d63 Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 5 Dec 2023 20:54:49 +0100 Subject: [PATCH 18/20] Auto stash before merge of "Kratos_RotatingFrameProcess" and "origin/Kratos_RotatingFrameProcess" --- .../add_geometrical_utilities_to_python.cpp | 2 + ..._symmetry_master_slave_constraints_process | 181 +++++++++ ...mmetry_master_slave_constraints_process.py | 182 +++++++++ ...lave_constraints_to_neighbours_utility.cpp | 383 +++++++++++++++++- ..._slave_constraints_to_neighbours_utility.h | 76 +++- 5 files changed, 819 insertions(+), 5 deletions(-) create mode 100644 kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process create mode 100644 kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py diff --git a/kratos/python/add_geometrical_utilities_to_python.cpp b/kratos/python/add_geometrical_utilities_to_python.cpp index f2913ad0c6c7..82972ff5c113 100644 --- a/kratos/python/add_geometrical_utilities_to_python.cpp +++ b/kratos/python/add_geometrical_utilities_to_python.cpp @@ -439,6 +439,8 @@ void AddGeometricalUtilitiesToPython(pybind11::module &m) .def(py::init()) .def("AssignMasterSlaveConstraintsToNodes", [](AssignMasterSlaveConstraintsToNeighboursUtility& rAssignMasterSlaveConstraintsToNeighboursUtility, NodesContainerType pNodes, double const Radius, ModelPart& rComputingModelPart, const std::vector>>& rVariableList, double const MinNumOfNeighNodes){ return rAssignMasterSlaveConstraintsToNeighboursUtility.AssignMasterSlaveConstraintsToNodes(pNodes, Radius, rComputingModelPart, rVariableList, MinNumOfNeighNodes);}) + .def("AssignRotationalMasterSlaveConstraintsToNodes", [](AssignMasterSlaveConstraintsToNeighboursUtility& rAssignMasterSlaveConstraintsToNeighboursUtility, NodesContainerType pNodes, double const Radius, ModelPart& rComputingModelPart, const std::vector>>& rVariableList, double const MinNumOfNeighNodes, const std::unordered_map>& GhostToOriginalMapping){ + return rAssignMasterSlaveConstraintsToNeighboursUtility.AssignRotationalMasterSlaveConstraintsToNodes(pNodes, Radius, rComputingModelPart, rVariableList, MinNumOfNeighNodes, GhostToOriginalMapping);}) ; } diff --git a/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process b/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process new file mode 100644 index 000000000000..720de7016de3 --- /dev/null +++ b/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process @@ -0,0 +1,181 @@ +import KratosMultiphysics as KM +import numpy as np +import math + +def Factory(settings, model): + if not isinstance(settings, KM.Parameters): + raise Exception("expected input shall be a Parameters object, encapsulating a json string") + return AssignRotationalSymmetryMasterSlaveConstraintsProcess(model, settings["Parameters"]) + +## All the processes python should be derived from "Process" +class AssignRotationalSymmetryMasterSlaveConstraintsProcess(KM.Process): + """The process facilitates the discovery of neighboring nodes in a + master model part within a designated radius for each node in the + slave model part. Following this, it establishes a master-slave constraint + and calculates its weights using a spatial function that employs radial basis functions. + + Public member variables: + model -- the container of the different model parts. + settings -- Kratos parameters containing process settings. + """ + + def __init__(self, model, settings): + """ The default constructor of the class + + Keyword arguments: + self -- It signifies an instance of a class. + Model -- the container of the different model parts. + settings -- Kratos parameters containing process settings. + """ + KM.Process.__init__(self) # calling the baseclass constructor + + default_settings = KM.Parameters("""{ + "model_part_name": "", + "slave_model_part_name": "", + "master_model_part_name": "", + "variable_names": [], + "search_radius": 1.0, + "minimum_number_of_neighbouring_nodes": 3, + "reform_constraints_at_each_step": false + }""") + + # Add missing settings that the user did not provide but that + # are necessary for this process + settings.ValidateAndAssignDefaults(default_settings) + + # Get the model part on which the MasterSlaveConstraints are going to be applied + if not settings["model_part_name"].GetString(): + raise Exception("\'model_part_name\' not provided. Please specify the model part to apply to MasterSlaveConstraints to.") + model_part_name = settings["model_part_name"].GetString() #MasterSlaveConstraints are applied to computing model part + self.computing_model_part = model.GetModelPart(model_part_name) + + # Get the slave model part + if not settings["slave_model_part_name"].GetString(): + raise Exception("\'slave_model_part_name\' not provided. Please specify the slave model part.") + slave_model_part_name = settings["slave_model_part_name"].GetString() + self.slave_model_part = model.GetModelPart(slave_model_part_name) + + # Get the master model part + if not settings["master_model_part_name"].GetString(): + raise Exception("\'master_model_part_name\' not provided. Please specify the master model part.") + master_model_part_name = settings["master_model_part_name"].GetString() + self.master_model_part = model.GetModelPart(master_model_part_name) + + # Search radius for the MasterSlaveConstraints + self.search_radius = settings["search_radius"].GetDouble() + + # Minimum number of neighboring nodes retrieved from the search radius + self.minimum_number_of_neighbouring_nodes = settings["minimum_number_of_neighbouring_nodes"].GetInt() + + # Apply MasterSlaveConstraints at each time step (True) or only once (False) + self.reform_constraints_at_each_step = settings["reform_constraints_at_each_step"].GetBool() + + # Retrieve and check if variables exist + variable_names = settings["variable_names"].GetStringArray() + if len(variable_names) == 0: + err_msg = "The variable names need to be specified by the user in the \'variable_names\' string array." + raise Exception(err_msg) + if any(variable_names.count(var_name) > 1 for var_name in variable_names): + err_msg = "There are repeated variables in the \'variable_names\' string array." + raise Exception(err_msg) + variable_names.sort() + + self.variables_list = [] # Initialize the list of variables + + for var_name in variable_names: + # Check if the variable exists in KratosGlobals + if not KM.KratosGlobals.HasVariable(var_name): + err_msg = "\'{}\' variable in \'variable_names\' is not in KratosGlobals. Please check the provided value.".format(var_name) + + var_type = KM.KratosGlobals.GetVariableType(var_name) # Get the type of the variable + + # Check the variable type and add it to the variables_list accordingly + if var_type == "Array": + domain_size = self.computing_model_part.ProcessInfo[KM.DOMAIN_SIZE] # Get the domain size from the ProcessInfo + component_suffixes = ["_X", "_Y", "_Z"] # Suffixes for the components of the array variable + for i in range(domain_size): + var_name_with_suffix = f"{var_name}{component_suffixes[i]}" # Append the component suffix to the variable name + self.variables_list.append(KM.KratosGlobals.GetVariable(var_name_with_suffix)) # Add the variable to the list + elif var_type == "Double": + self.variables_list.append(KM.KratosGlobals.GetVariable(var_name)) # Add the variable to the list + else: + raise Exception("Variable " + var_name + " not compatible") # Raise an exception for incompatible variable types + + # Initialization of ghost nodes related attributes + self.model.CreateModelPart("AuxiliarMasterModelPart") + self.aux_master_model_part = self.model.GetModelPart("AuxiliarMasterModelPart") + self.max_node_id = self.computing_model_part.NumberOfNodes() + + def ExecuteInitialize(self): + """ This method is executed at the begining to initialize the process + + Keyword arguments: + self -- It signifies an instance of a class. + """ + + self.__CreateGhostNodes(self.master_model_part.Nodes) + + # Initialize MasterSlaveConstraints to neighbours utility + self.assign_mscs_utility = KM.AssignMasterSlaveConstraintsToNeighboursUtility(self.master_model_part.Nodes) + + # The user may only need to set up the MasterSlaveConstraints only once + if not self.reform_constraints_at_each_step: + for variable in self.variables_list: + self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes,self.search_radius,self.computing_model_part, variable, self.minimum_number_of_neighbouring_nodes) + + + def ExecuteInitializeSolutionStep(self): + """ This method is executed in order to initialize the current step + + Keyword arguments: + self -- It signifies an instance of a class. + """ + # If the user want the mscs to be updated at each time step, this is usefull for moving meshes. + if self.reform_constraints_at_each_step: + self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes,self.search_radius,self.computing_model_part, self.variables_list, self.minimum_number_of_neighbouring_nodes) + + def ExecuteFinalizeSolutionStep(self): + """ This method is executed in order to finalize the current step + + Keyword arguments: + self -- It signifies an instance of a class. + """ + # If MasterSlaveConstraints are updated every time step, these are to me removed before being re-assigned. + if self.reform_constraints_at_each_step: + self.__RemoveConstraints() + + def __RemoveConstraints(self): + #Remove master-slave constraints + KM.VariableUtils().SetFlag(KM.TO_ERASE, True, self.computing_model_part.MasterSlaveConstraints) + self.computing_model_part.RemoveMasterSlaveConstraintsFromAllLevels(KM.TO_ERASE) + + def __CreateGhostNodes(self, nodes): + """Create ghost nodes for each node.""" + ghost_to_original_mapping = {} + rotation_angles = [0, 60, 120, 180, 240, 300] + + self.max_node_id = self.main_model_part.NumberOfNodes() + + for angle in rotation_angles: + # Create ghost nodes + for original_node in nodes: + ghost_node = self.__RotateNode(original_node, angle) + ghost_to_original_mapping[ghost_node] = original_node + self.aux_master_model_part.AddNode(ghost_node) + + return ghost_to_original_mapping + + def __RotateNode(self, original_node, angle_degrees): + """Rotate a node by a given angle (in degrees) about the origin.""" + angle_radians = math.radians(angle_degrees) + + # Apply 2D rotation matrix + new_x = original_node.X * math.cos(angle_radians) - original_node.Y * math.sin(angle_radians) + new_y = original_node.X * math.sin(angle_radians) + original_node.Y * math.cos(angle_radians) + + # Generate a new node ID based on the maximum node ID + new_node_id = self.max_node_id + 1 + self.max_node_id += 1 + + # Create the new node with the ID and coordinates + return KM.Node(new_node_id, new_x, new_y, 0.0) \ No newline at end of file diff --git a/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py b/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py new file mode 100644 index 000000000000..54f4824552e3 --- /dev/null +++ b/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py @@ -0,0 +1,182 @@ +import KratosMultiphysics as KM +import numpy as np +import math + +def Factory(settings, model): + if not isinstance(settings, KM.Parameters): + raise Exception("expected input shall be a Parameters object, encapsulating a json string") + return AssignRotationalSymmetryMasterSlaveConstraintsProcess(model, settings["Parameters"]) + +## All the processes python should be derived from "Process" +class AssignRotationalSymmetryMasterSlaveConstraintsProcess(KM.Process): + """The process facilitates the discovery of neighboring nodes in a + master model part within a designated radius for each node in the + slave model part. Following this, it establishes a master-slave constraint + and calculates its weights using a spatial function that employs radial basis functions. + + Public member variables: + model -- the container of the different model parts. + settings -- Kratos parameters containing process settings. + """ + + def __init__(self, model, settings): + """ The default constructor of the class + + Keyword arguments: + self -- It signifies an instance of a class. + Model -- the container of the different model parts. + settings -- Kratos parameters containing process settings. + """ + KM.Process.__init__(self) # calling the baseclass constructor + + default_settings = KM.Parameters("""{ + "model_part_name": "", + "slave_model_part_name": "", + "master_model_part_name": "", + "variable_names": [], + "search_radius": 1.0, + "minimum_number_of_neighbouring_nodes": 3 + }""") + + # Add missing settings that the user did not provide but that + # are necessary for this process + settings.ValidateAndAssignDefaults(default_settings) + + # Get the model part on which the MasterSlaveConstraints are going to be applied + if not settings["model_part_name"].GetString(): + raise Exception("\'model_part_name\' not provided. Please specify the model part to apply to MasterSlaveConstraints to.") + model_part_name = settings["model_part_name"].GetString() #MasterSlaveConstraints are applied to computing model part + self.computing_model_part = model.GetModelPart(model_part_name) + + # Get the slave model part + if not settings["slave_model_part_name"].GetString(): + raise Exception("\'slave_model_part_name\' not provided. Please specify the slave model part.") + slave_model_part_name = settings["slave_model_part_name"].GetString() + self.slave_model_part = model.GetModelPart(slave_model_part_name) + + # Get the master model part + if not settings["master_model_part_name"].GetString(): + raise Exception("\'master_model_part_name\' not provided. Please specify the master model part.") + master_model_part_name = settings["master_model_part_name"].GetString() + self.master_model_part = model.GetModelPart(master_model_part_name) + + # Search radius for the MasterSlaveConstraints + self.search_radius = settings["search_radius"].GetDouble() + + # Minimum number of neighboring nodes retrieved from the search radius + self.minimum_number_of_neighbouring_nodes = settings["minimum_number_of_neighbouring_nodes"].GetInt() + + # Retrieve and check if variables exist + variable_names = settings["variable_names"].GetStringArray() + if len(variable_names) == 0: + err_msg = "The variable names need to be specified by the user in the \'variable_names\' string array." + raise Exception(err_msg) + if any(variable_names.count(var_name) > 1 for var_name in variable_names): + err_msg = "There are repeated variables in the \'variable_names\' string array." + raise Exception(err_msg) + variable_names.sort() + + self.variables_list = [] # Initialize the list of variables + + for var_name in variable_names: + # Check if the variable exists in KratosGlobals + if not KM.KratosGlobals.HasVariable(var_name): + err_msg = "\'{}\' variable in \'variable_names\' is not in KratosGlobals. Please check the provided value.".format(var_name) + + var_type = KM.KratosGlobals.GetVariableType(var_name) # Get the type of the variable + + # Check the variable type and add it to the variables_list accordingly + if var_type == "Array": + domain_size = self.computing_model_part.ProcessInfo[KM.DOMAIN_SIZE] # Get the domain size from the ProcessInfo + component_suffixes = ["_X", "_Y", "_Z"] # Suffixes for the components of the array variable + for i in range(domain_size): + var_name_with_suffix = f"{var_name}{component_suffixes[i]}" # Append the component suffix to the variable name + self.variables_list.append(KM.KratosGlobals.GetVariable(var_name_with_suffix)) # Add the variable to the list + elif var_type == "Double": + self.variables_list.append(KM.KratosGlobals.GetVariable(var_name)) # Add the variable to the list + else: + raise Exception("Variable " + var_name + " not compatible") # Raise an exception for incompatible variable types + + # Initialization of ghost nodes related attributes + model.CreateModelPart("AuxiliarMasterModelPart") + self.aux_master_model_part = model.GetModelPart("AuxiliarMasterModelPart") + self.max_node_id = self.computing_model_part.NumberOfNodes() + + def ExecuteInitialize(self): + """ This method is executed at the begining to initialize the process + + Keyword arguments: + self -- It signifies an instance of a class. + """ + + # Create ghost nodes in the master model part + self.__CreateGhostNodes(self.master_model_part.Nodes) + + # Initialize MasterSlaveConstraints to neighbours utility with the master model part nodes (which now includes ghost nodes) + self.assign_mscs_utility = KM.AssignMasterSlaveConstraintsToNeighboursUtility(self.aux_master_model_part.Nodes) + + def ExecuteInitializeSolutionStep(self): + """ This method is executed in order to initialize the current step + + Keyword arguments: + self -- It signifies an instance of a class. + """ + # Assign dynamically the master-slave constraints + # self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes, self.search_radius, self.computing_model_part, self.variables_list, self.minimum_number_of_neighbouring_nodes) + self.assign_mscs_utility.AssignRotationalMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes, self.search_radius, self.computing_model_part, self.variables_list, self.minimum_number_of_neighbouring_nodes, self.ghost_to_original_mapping) + + def ExecuteFinalizeSolutionStep(self): + """ This method is executed in order to finalize the current step + + Keyword arguments: + self -- It signifies an instance of a class. + """ + # Since MasterSlaveConstraints are updated every time step, they need to be removed before being re-assigned. + self.__RemoveConstraints() + + def __RemoveConstraints(self): + #Remove master-slave constraints + KM.VariableUtils().SetFlag(KM.TO_ERASE, True, self.computing_model_part.MasterSlaveConstraints) + self.computing_model_part.RemoveMasterSlaveConstraintsFromAllLevels(KM.TO_ERASE) + + def __CreateGhostNodes(self, nodes): + """Create ghost nodes for each node by rotating them about the origin. + + Keyword arguments: + nodes -- The nodes to be rotated to create the ghost nodes. + """ + self.ghost_to_original_mapping = {} # Initialize the mapping as an attribute of the class + rotation_angles = [0, -60, -120, -180, -240, -300] + self.max_node_id = self.computing_model_part.NumberOfNodes() + + for angle in rotation_angles: + # Create ghost nodes + for original_node in nodes: + ghost_node = self.__RotateNode(original_node, angle) + # Store the mapping using node IDs and the rotation angle + self.ghost_to_original_mapping[ghost_node.Id] = (original_node.Id, np.pi*angle/180) + if original_node.Id==922: + print(angle) + print(ghost_node.Id) + self.aux_master_model_part.AddNode(ghost_node) + + + def __RotateNode(self, original_node, angle_degrees): + """Rotate a node by a given angle (in degrees) about the origin. + + Keyword arguments: + original_node -- The node to be rotated. + angle_degrees -- The angle in degrees by which the node is to be rotated. + """ + angle_radians = math.radians(angle_degrees) + + # Apply 2D rotation matrix + new_x = original_node.X * math.cos(angle_radians) - original_node.Y * math.sin(angle_radians) + new_y = original_node.X * math.sin(angle_radians) + original_node.Y * math.cos(angle_radians) + + # Generate a new node ID based on the maximum node ID + new_node_id = self.max_node_id + 1 + self.max_node_id += 1 + + # Create the new node with the ID and coordinates + return KM.Node(new_node_id, new_x, new_y, 0.0) \ No newline at end of file diff --git a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp index dbfd9b00f0c1..6c2a92ad7f0e 100644 --- a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp +++ b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp @@ -131,7 +131,7 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::GetDofsAndCoordinatesForCl // Loop over the cloud of nodes for (std::size_t i = 0; i < CloudOfNodesArray.size(); i++) { - auto p_slave_node = CloudOfNodesArray[i]; + auto p_slave_node = CloudOfNodesArray[i];//////CHANGE SLAVE NODE NAME // Loop over the variable list for (std::size_t j = 0; j < rVariableList.size(); j++) { @@ -154,6 +154,99 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::GetDofsAndCoordinatesForCl KRATOS_CATCH(""); } +void AssignMasterSlaveConstraintsToNeighboursUtility::GetDofsAndCoordinatesForCloudOfNodesFlattened( + const ResultNodesContainerType& CloudOfNodesArray, + const std::vector>>& rVariableList, + DofPointerVectorType& rCloudOfDofs, + Matrix& rCloudOfNodesCoordinates +) +{ + KRATOS_TRY; + + // Resize rCloudOfDofs to accommodate all DOFs for all variables + rCloudOfDofs.resize(rVariableList.size() * CloudOfNodesArray.size()); + + // Resize rCloudOfNodesCoordinates to match the size of the cloud of nodes + rCloudOfNodesCoordinates.resize(CloudOfNodesArray.size(), 3); + + // Loop over the cloud of nodes + for (std::size_t i = 0; i < CloudOfNodesArray.size(); i++) { + auto p_slave_node = CloudOfNodesArray[i]; + + // Loop over the variable list + for (std::size_t j = 0; j < rVariableList.size(); j++) { + const auto& r_variable = rVariableList[j]; + + // Check if the slave node has the required DOF for the current variable + if (p_slave_node->HasDofFor(r_variable.get())) { + std::size_t counter = i * rVariableList.size() + j; + rCloudOfDofs[counter] = p_slave_node->pGetDof(r_variable.get()); + } else { + KRATOS_ERROR << "The node with ID " << p_slave_node->Id() << " does not have the required DOF for the variable " << r_variable.get().Name() << std::endl; + } + } + + // Assign the coordinates of the slave node to the corresponding row in rCloudOfNodesCoordinates + noalias(row(rCloudOfNodesCoordinates, i)) = p_slave_node->Coordinates(); + } + + KRATOS_CATCH(""); +} + +struct CustomComparator { + static constexpr double TOLERANCE = 1e-6; + + bool operator()(const std::tuple& a, const std::tuple& b) const { + return std::fabs(std::get<0>(a) - std::get<0>(b)) > TOLERANCE || + std::fabs(std::get<1>(a) - std::get<1>(b)) > TOLERANCE || + std::fabs(std::get<2>(a) - std::get<2>(b)) > TOLERANCE; + } +}; + +void AssignMasterSlaveConstraintsToNeighboursUtility::FilterUniqueNodesForRBF( + ResultNodesContainerType& CloudOfNodesArray, + Matrix& rCloudOfNodesCoordinates +) { + KRATOS_TRY; + + // Use a set with a custom comparator to keep track of unique coordinates + std::set, CustomComparator> unique_coords; + + // Temporary vector to store unique nodes + std::vector temp_nodes; + + // Loop over the cloud of nodes + for (auto& p_node : CloudOfNodesArray) { + auto coords = p_node->Coordinates(); + std::tuple current_coords(coords[0], coords[1], coords[2]); + + // Check if the coordinates are unique + if (unique_coords.find(current_coords) == unique_coords.end()) { + unique_coords.insert(current_coords); + temp_nodes.push_back(p_node); + } + } + + // Clear the original CloudOfNodesArray and repopulate with unique nodes + CloudOfNodesArray.clear(); + for (auto& p_node : temp_nodes) { + CloudOfNodesArray.push_back(p_node); + } + + // Resize rCloudOfNodesCoordinates to match the size of unique nodes + rCloudOfNodesCoordinates.resize(temp_nodes.size(), 3); + + // Fill rCloudOfNodesCoordinates with unique coordinates + for (std::size_t i = 0; i < temp_nodes.size(); i++) { + auto coords = temp_nodes[i]->Coordinates(); + for (std::size_t j = 0; j < 3; j++) { + rCloudOfNodesCoordinates(i, j) = coords[j]; + } + } + + KRATOS_CATCH(""); +} + void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstraintsToNodes( NodesContainerType pSlaveNodes, double const Radius, @@ -216,6 +309,12 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstrain constant_vector.resize(r_n_container.size()); IndexType it = i.fetch_add(1) * rVariableList.size(); // Atomically increment the counter and get the previous value + if (p_slave_node->Id() == 929){ + KRATOS_WATCH(r_cloud_of_nodes) + KRATOS_WATCH(r_n_container) + KRATOS_WATCH(shape_matrix) + } + for (std::size_t j = 0; j < rVariableList.size(); ++j) { concurrent_constraints.enqueue(r_clone_constraint.Create(prev_num_mscs + it + j + 1, const_cast(r_cloud_of_dofs[j]), r_slave_dofs[j], shape_matrix, constant_vector)); } @@ -236,6 +335,288 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstrain KRATOS_CATCH(""); } +void AssignMasterSlaveConstraintsToNeighboursUtility::AssignRotationalMasterSlaveConstraintsToNodes( + NodesContainerType pSlaveNodes, + double const Radius, + ModelPart& rComputingModelPart, + const std::vector>>& rVariableList, + double const MinNumOfNeighNodes, + const std::unordered_map>& GhostToOriginalMapping +) +{ + KRATOS_TRY; + + BuiltinTimer build_and_assign_mscs; + + thread_local std::vector r_cloud_of_dofs; + thread_local DofPointerVectorType r_cloud_of_master_dofs; + thread_local std::vector r_slave_dofs; + thread_local Matrix r_cloud_of_nodes_coordinates; + thread_local Matrix r_cloud_of_master_nodes_coordinates; + thread_local array_1d r_slave_coordinates; + thread_local Vector r_n_container; + thread_local Vector constant_vector; + thread_local Matrix shape_matrix_x, shape_matrix_y, shape_matrix_z; + + int prev_num_mscs = rComputingModelPart.NumberOfMasterSlaveConstraints(); + + KRATOS_WARNING_IF("AssignRotationalMasterSlaveConstraintsToNodes", prev_num_mscs > 0) << "Previous Master-Slave Constraints exist in the ModelPart. The new Constraints may interact with the existing ones." << std::endl; + + auto& p_nodes_array = pSlaveNodes.GetContainer(); + + const auto& r_clone_constraint = KratosComponents::Get("LinearMasterSlaveConstraint"); + + using MasterSlaveConstraintQueue = moodycamel::ConcurrentQueue; + MasterSlaveConstraintQueue concurrent_constraints; + + std::atomic i(0); + + block_for_each(p_nodes_array, [&](Node::Pointer& p_slave_node) { + + // Skip nodes with IDs 1174 or 1015 950 1111 + if (p_slave_node->Id() == 1174 || p_slave_node->Id() == 1015) { + return; // Skip the current iteration of the loop + } + + double local_radius = Radius; + ResultNodesContainerType r_cloud_of_nodes(0); + ResultNodesContainerType r_local_results(mMaxNumberOfNodes); + + while (r_cloud_of_nodes.size() < MinNumOfNeighNodes) { + r_cloud_of_nodes.clear(); + SearchNodesInRadiusForNode(p_slave_node, local_radius, r_cloud_of_nodes, r_local_results); + local_radius += Radius; + } + + GetDofsAndCoordinatesForSlaveNode(p_slave_node, rVariableList, r_slave_dofs, r_slave_coordinates); + FilterUniqueNodesForRBF(r_cloud_of_nodes, r_cloud_of_nodes_coordinates); + + ResultNodesContainerType r_cloud_of_master_nodes; + r_cloud_of_master_nodes.reserve(r_cloud_of_nodes.size()); + for (auto& ghost_node : r_cloud_of_nodes) { + auto it = GhostToOriginalMapping.find(ghost_node->Id()); + if (it != GhostToOriginalMapping.end()) { + IndexType master_node_id = it->second.first; // Extracting the original/master node ID from the mapping + Node::Pointer p_master_node = rComputingModelPart.pGetNode(master_node_id); + r_cloud_of_master_nodes.push_back(p_master_node); + } else { + KRATOS_ERROR << "Ghost Node ID " << ghost_node->Id() << " not found in GhostToOriginalMapping!" << std::endl; + } + } + + GetDofsAndCoordinatesForCloudOfNodesFlattened(r_cloud_of_master_nodes, rVariableList, r_cloud_of_master_dofs, r_cloud_of_master_nodes_coordinates); + + RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_slave_coordinates, r_n_container); + + Vector ListOfAngles; + ExtractAnglesFromMapping(r_cloud_of_nodes, GhostToOriginalMapping, ListOfAngles); + + ObtainShapeMatrix(shape_matrix_x, shape_matrix_y, shape_matrix_z, r_n_container, ListOfAngles); + constant_vector.resize(shape_matrix_x.size2()); + if (p_slave_node->Id() == 929){ + KRATOS_WATCH(r_cloud_of_nodes) + KRATOS_WATCH(r_cloud_of_nodes_coordinates) + KRATOS_WATCH(r_slave_coordinates) + KRATOS_WATCH(r_cloud_of_master_nodes) + KRATOS_WATCH(r_cloud_of_master_dofs) + KRATOS_WATCH(r_n_container) + KRATOS_WATCH(ListOfAngles) + KRATOS_WATCH(shape_matrix_x) + KRATOS_WATCH(shape_matrix_y) + KRATOS_WATCH(shape_matrix_z) + } + + IndexType it = i.fetch_add(1) * rVariableList.size(); + + std::vector shape_matrices = {shape_matrix_x, shape_matrix_y, shape_matrix_z}; + + for (std::size_t j = 0; j < rVariableList.size(); ++j) { + if (j >= shape_matrices.size()) { + KRATOS_ERROR << "Unexpected index in rVariableList. Only 3 DoFs (x, y, z) are supported." << std::endl; + } + concurrent_constraints.enqueue(r_clone_constraint.Create(prev_num_mscs + it + j + 1, r_cloud_of_master_dofs, r_slave_dofs[j], shape_matrices[j], constant_vector)); + } + + + }); + + ConstraintContainerType temp_constraints; + ModelPart::MasterSlaveConstraintType::Pointer p_constraint; + + while (concurrent_constraints.try_dequeue(p_constraint)) + temp_constraints.push_back(p_constraint); + + rComputingModelPart.AddMasterSlaveConstraints(temp_constraints.begin(), temp_constraints.end()); + + KRATOS_INFO("AssignRotationalMasterSlaveConstraintsToNodes") << "Build and Assign Rotational Master-Slave Constraints Time: " << build_and_assign_mscs.ElapsedSeconds() << std::endl; + + KRATOS_CATCH(""); +} + +void AssignMasterSlaveConstraintsToNeighboursUtility::ObtainShapeMatrix( + Matrix& shape_matrix_x, + Matrix& shape_matrix_y, + Matrix& shape_matrix_z, + Vector& rN_container, + Vector& ListOfAngles +) +{ + std::vector rotation_x_global; + std::vector rotation_y_global; + std::vector rotation_z_global; + std::vector rotation_x(3); + std::vector rotation_y(3); + std::vector rotation_z(3); + for (int i=0; i < static_cast(rN_container.size()); i++){ + //TODO: Any kind of rotation. x-axis, y-axis, z-axis. Make it a function. + rotation_x[0] = rN_container[i]*cos(ListOfAngles[i]); + rotation_x[1] = -rN_container[i]*sin(ListOfAngles[i]); + rotation_x[2] = rN_container[i]*0.0; + rotation_y[0] = rN_container[i]*sin(ListOfAngles[i]); + rotation_y[1] = rN_container[i]*cos(ListOfAngles[i]); + rotation_y[2] = rN_container[i]*0.0; + rotation_z[0] = rN_container[i]*0.0; + rotation_z[1] = rN_container[i]*0.0; + rotation_z[2] = rN_container[i]*1.0; + // noalias(row(shape_matrix,i)) = aux_vector_x; + rotation_x_global.insert(rotation_x_global.end(),rotation_x.begin(),rotation_x.end()); + rotation_y_global.insert(rotation_y_global.end(),rotation_y.begin(),rotation_y.end()); + rotation_z_global.insert(rotation_z_global.end(),rotation_z.begin(),rotation_z.end()); + } + shape_matrix_x.resize(1, rotation_x_global.size()); + shape_matrix_y.resize(1, rotation_y_global.size()); + shape_matrix_z.resize(1, rotation_z_global.size()); + for (int i=0; i(rotation_x_global.size());i++){ + shape_matrix_x(0,i) = rotation_x_global[i]; + shape_matrix_y(0,i) = rotation_y_global[i]; + shape_matrix_z(0,i) = rotation_z_global[i]; + } +} + +// void AssignMasterSlaveConstraintsToNeighboursUtility::AssignRotationalMasterSlaveConstraintsToNodes( +// NodesContainerType pSlaveNodes, +// double const Radius, +// ModelPart& rComputingModelPart, +// const std::vector>>& rVariableList, +// double const MinNumOfNeighNodes, +// const std::unordered_map>& GhostToOriginalMapping +// ) +// { +// KRATOS_TRY; + +// BuiltinTimer build_and_assign_mscs; + +// thread_local std::vector r_cloud_of_dofs; +// thread_local std::vector r_slave_dofs; +// thread_local Matrix r_cloud_of_nodes_coordinates; +// thread_local array_1d r_slave_coordinates; +// thread_local Vector r_n_container; +// thread_local Vector constant_vector; +// thread_local Matrix shape_matrix; + +// int prev_num_mscs = rComputingModelPart.NumberOfMasterSlaveConstraints(); + +// KRATOS_WARNING_IF("AssignRotationalMasterSlaveConstraintsToNodes", prev_num_mscs > 0) << "Previous Master-Slave Constraints exist in the ModelPart. The new Constraints may interact with the existing ones." << std::endl; + +// auto& p_nodes_array = pSlaveNodes.GetContainer(); + +// const auto& r_clone_constraint = KratosComponents::Get("LinearMasterSlaveConstraint"); + +// using MasterSlaveConstraintQueue = moodycamel::ConcurrentQueue; +// MasterSlaveConstraintQueue concurrent_constraints; + +// std::atomic i(0); + +// block_for_each(p_nodes_array, [&](Node::Pointer& p_slave_node) { + +// double local_radius = Radius; +// ResultNodesContainerType r_cloud_of_nodes(0); +// ResultNodesContainerType r_local_results(mMaxNumberOfNodes); + +// while (r_cloud_of_nodes.size() < MinNumOfNeighNodes) { +// r_cloud_of_nodes.clear(); +// SearchNodesInRadiusForNode(p_slave_node, local_radius, r_cloud_of_nodes, r_local_results); +// local_radius += Radius; +// } + +// GetDofsAndCoordinatesForSlaveNode(p_slave_node, rVariableList, r_slave_dofs, r_slave_coordinates); +// GetDofsAndCoordinatesForCloudOfNodes(r_cloud_of_nodes, rVariableList, r_cloud_of_dofs, r_cloud_of_nodes_coordinates); + +// RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_slave_coordinates, r_n_container); + +// shape_matrix.resize(1, r_n_container.size()); +// noalias(row(shape_matrix, 0)) = r_n_container; +// constant_vector.resize(r_n_container.size()); + +// // Check if the current slave node is a ghost node and adjust for rotation if necessary +// if (GhostToOriginalMapping.find(p_slave_node->Id()) != GhostToOriginalMapping.end()) { +// auto mapping_data = GhostToOriginalMapping.at(p_slave_node->Id()); +// double rotation_angle = mapping_data.second; +// AdjustConstraintForRotation(rotation_angle, shape_matrix); +// } + +// IndexType it = i.fetch_add(1) * rVariableList.size(); + +// for (std::size_t j = 0; j < rVariableList.size(); ++j) { +// concurrent_constraints.enqueue(r_clone_constraint.Create(prev_num_mscs + it + j + 1, const_cast(r_cloud_of_dofs[j]), r_slave_dofs[j], shape_matrix, constant_vector)); +// } + +// }); + +// ConstraintContainerType temp_constraints; +// ModelPart::MasterSlaveConstraintType::Pointer p_constraint; + +// while (concurrent_constraints.try_dequeue(p_constraint)) +// temp_constraints.push_back(p_constraint); + +// rComputingModelPart.AddMasterSlaveConstraints(temp_constraints.begin(), temp_constraints.end()); + +// KRATOS_INFO("AssignRotationalMasterSlaveConstraintsToNodes") << "Build and Assign Rotational Master-Slave Constraints Time: " << build_and_assign_mscs.ElapsedSeconds() << std::endl; + +// KRATOS_CATCH(""); +// } + +void AssignMasterSlaveConstraintsToNeighboursUtility::AdjustConstraintForRotation(double rotation_angle, Matrix& shape_matrix) +{ + KRATOS_TRY; + + double cos_theta = std::cos(rotation_angle * M_PI / 180.0); // Convert degrees to radians + double sin_theta = std::sin(rotation_angle * M_PI / 180.0); + + for (std::size_t i = 0; i < shape_matrix.size1(); ++i) { + double original_w1 = shape_matrix(i, 0); + double original_w2 = shape_matrix(i, 1); + + shape_matrix(i, 0) = original_w1 * cos_theta - original_w2 * sin_theta; + shape_matrix(i, 1) = original_w1 * sin_theta + original_w2 * cos_theta; + } + + KRATOS_CATCH(""); +} + +void AssignMasterSlaveConstraintsToNeighboursUtility::ExtractAnglesFromMapping( + const ResultNodesContainerType& Results, + const std::unordered_map>& GhostToOriginalMapping, + Vector& ListOfAngles +) +{ + ListOfAngles.resize(Results.size(), false); + for(int i = 0; i < static_cast(Results.size()); i++) + { + auto it = GhostToOriginalMapping.find(Results[i]->Id()); + if(it != GhostToOriginalMapping.end()) + { + ListOfAngles[i] = it->second.second; // Extracting the angle from the mapping + } + else + { + KRATOS_ERROR << "Node ID " << Results[i]->Id() << " not found in GhostToOriginalMapping!" << std::endl; + } + } +} + +} // namespace Kratos. + void AssignMasterSlaveConstraintsToNeighboursUtility::InterpolateVelocityToNodes( NodesContainerType pDestinationNodes, double const Radius, diff --git a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h index e4c66cbd225b..754e74aa4ccd 100644 --- a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h +++ b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h @@ -137,6 +137,29 @@ class KRATOS_API(KRATOS_CORE) AssignMasterSlaveConstraintsToNeighboursUtility Matrix& rCloudOfNodesCoordinates ); + /** + * @brief Removes duplicate nodes based on their coordinates. + * + * This function is intended to filter out nodes with duplicate coordinates + * to avoid errors when using Radial Basis Functions (RBF) that can't handle + * repetitive coordinates. The function updates the provided CloudOfNodesArray + * to only contain unique nodes and also returns a matrix of their coordinates. + * + * @param CloudOfNodesArray A reference to the container of nodes. This will be updated to only contain unique nodes. + * @param rCloudOfNodesCoordinates A reference to a matrix that will be filled with the coordinates of the unique nodes. + */ + void FilterUniqueNodesForRBF( + ResultNodesContainerType& CloudOfNodesArray, + Matrix& rCloudOfNodesCoordinates + ); + + void GetDofsAndCoordinatesForCloudOfNodesFlattened( + const ResultNodesContainerType& CloudOfNodesArray, + const std::vector>>& rVariableList, + DofPointerVectorType& rCloudOfDofs, + Matrix& rCloudOfNodesCoordinates + ); + /** * @brief Assign Master-Slave Constraints to a set of Nodes. * @details Assign Matser-Slave Constraints to a set of Nodes given a radius of influence @@ -155,12 +178,57 @@ class KRATOS_API(KRATOS_CORE) AssignMasterSlaveConstraintsToNeighboursUtility double const MinNumOfNeighNodes ); - void InterpolateVelocityToNodes( - NodesContainerType pDestinationNodes, + /** + * @brief Assign Rotational Master-Slave Constraints to a set of Nodes. + * @details This function assigns Master-Slave Constraints to a set of nodes considering rotational effects. + * It takes into account the rotation of the nodes and adjusts the constraints accordingly. The function + * uses MLS or RBF shape functions within a given radius of influence to determine the constraints. + * Additionally, it uses a mapping between ghost nodes and their original counterparts to account for the rotation. + * @param pSlaveNodes Nodes to which MasterSlaveConstraints are assigned. + * @param Radius Search radius for neighboring nodes. + * @param rComputingModelPart Model Part to which MasterSlaveConstraints are applied. + * @param rVariableList List of DOFs to assign the MasterSlaveConstraints. + * @param MinNumOfNeighNodes Minimum number of neighboring nodes required for the constraint. + * @param GhostToOriginalMapping Mapping between ghost nodes and their original counterparts, including the rotation angle. + */ + void AssignRotationalMasterSlaveConstraintsToNodes( + NodesContainerType pSlaveNodes, double const Radius, + ModelPart& rComputingModelPart, const std::vector>>& rVariableList, - double const MinNumOfNeighNodes - ); + double const MinNumOfNeighNodes, + const std::unordered_map>& GhostToOriginalMapping + ); + /** + * @brief Adjusts the Master-Slave Constraint for rotational effects. + * @details This function modifies the shape functions matrix (weights) of the Master-Slave Constraint + * to account for the rotation of the nodes. Given a rotation angle, it applies the necessary + * trigonometric transformations to the shape functions to ensure that the constraints are correctly + * adjusted for the rotation. This is crucial when dealing with models that have rotational dynamics + * or when considering ghost nodes that represent rotated versions of the original nodes. + * @param rotation_angle The angle (in degrees) by which the node has been rotated. + * @param shape_matrix The shape functions matrix that needs adjustment for rotation. + */ + void AdjustConstraintForRotation( + double rotation_angle, + Matrix& shape_matrix + ); + + void ExtractAnglesFromMapping( + const ResultNodesContainerType& Results, + const std::unordered_map>& GhostToOriginalMapping, + Vector& ListOfAngles + ); + + void ObtainShapeMatrix( + Matrix& shape_matrix_x, + Matrix& shape_matrix_y, + Matrix& shape_matrix_z, + Vector& rN_container, + Vector& ListOfAngles + ); + + ///@} ///@name Input and output From 37240f878d42d71ce260e5280fc3bf354902908a Mon Sep 17 00:00:00 2001 From: SADPR Date: Tue, 5 Dec 2023 21:12:08 +0100 Subject: [PATCH 19/20] Revert "Auto stash before merge of "Kratos_RotatingFrameProcess" and "origin/Kratos_RotatingFrameProcess"" This reverts commit 5be37fedc1ab27dadab5387178a6f2d58f2e2d63. --- .../add_geometrical_utilities_to_python.cpp | 2 - ..._symmetry_master_slave_constraints_process | 181 --------- ...mmetry_master_slave_constraints_process.py | 182 --------- ...lave_constraints_to_neighbours_utility.cpp | 383 +----------------- ..._slave_constraints_to_neighbours_utility.h | 76 +--- 5 files changed, 5 insertions(+), 819 deletions(-) delete mode 100644 kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process delete mode 100644 kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py diff --git a/kratos/python/add_geometrical_utilities_to_python.cpp b/kratos/python/add_geometrical_utilities_to_python.cpp index 82972ff5c113..f2913ad0c6c7 100644 --- a/kratos/python/add_geometrical_utilities_to_python.cpp +++ b/kratos/python/add_geometrical_utilities_to_python.cpp @@ -439,8 +439,6 @@ void AddGeometricalUtilitiesToPython(pybind11::module &m) .def(py::init()) .def("AssignMasterSlaveConstraintsToNodes", [](AssignMasterSlaveConstraintsToNeighboursUtility& rAssignMasterSlaveConstraintsToNeighboursUtility, NodesContainerType pNodes, double const Radius, ModelPart& rComputingModelPart, const std::vector>>& rVariableList, double const MinNumOfNeighNodes){ return rAssignMasterSlaveConstraintsToNeighboursUtility.AssignMasterSlaveConstraintsToNodes(pNodes, Radius, rComputingModelPart, rVariableList, MinNumOfNeighNodes);}) - .def("AssignRotationalMasterSlaveConstraintsToNodes", [](AssignMasterSlaveConstraintsToNeighboursUtility& rAssignMasterSlaveConstraintsToNeighboursUtility, NodesContainerType pNodes, double const Radius, ModelPart& rComputingModelPart, const std::vector>>& rVariableList, double const MinNumOfNeighNodes, const std::unordered_map>& GhostToOriginalMapping){ - return rAssignMasterSlaveConstraintsToNeighboursUtility.AssignRotationalMasterSlaveConstraintsToNodes(pNodes, Radius, rComputingModelPart, rVariableList, MinNumOfNeighNodes, GhostToOriginalMapping);}) ; } diff --git a/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process b/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process deleted file mode 100644 index 720de7016de3..000000000000 --- a/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process +++ /dev/null @@ -1,181 +0,0 @@ -import KratosMultiphysics as KM -import numpy as np -import math - -def Factory(settings, model): - if not isinstance(settings, KM.Parameters): - raise Exception("expected input shall be a Parameters object, encapsulating a json string") - return AssignRotationalSymmetryMasterSlaveConstraintsProcess(model, settings["Parameters"]) - -## All the processes python should be derived from "Process" -class AssignRotationalSymmetryMasterSlaveConstraintsProcess(KM.Process): - """The process facilitates the discovery of neighboring nodes in a - master model part within a designated radius for each node in the - slave model part. Following this, it establishes a master-slave constraint - and calculates its weights using a spatial function that employs radial basis functions. - - Public member variables: - model -- the container of the different model parts. - settings -- Kratos parameters containing process settings. - """ - - def __init__(self, model, settings): - """ The default constructor of the class - - Keyword arguments: - self -- It signifies an instance of a class. - Model -- the container of the different model parts. - settings -- Kratos parameters containing process settings. - """ - KM.Process.__init__(self) # calling the baseclass constructor - - default_settings = KM.Parameters("""{ - "model_part_name": "", - "slave_model_part_name": "", - "master_model_part_name": "", - "variable_names": [], - "search_radius": 1.0, - "minimum_number_of_neighbouring_nodes": 3, - "reform_constraints_at_each_step": false - }""") - - # Add missing settings that the user did not provide but that - # are necessary for this process - settings.ValidateAndAssignDefaults(default_settings) - - # Get the model part on which the MasterSlaveConstraints are going to be applied - if not settings["model_part_name"].GetString(): - raise Exception("\'model_part_name\' not provided. Please specify the model part to apply to MasterSlaveConstraints to.") - model_part_name = settings["model_part_name"].GetString() #MasterSlaveConstraints are applied to computing model part - self.computing_model_part = model.GetModelPart(model_part_name) - - # Get the slave model part - if not settings["slave_model_part_name"].GetString(): - raise Exception("\'slave_model_part_name\' not provided. Please specify the slave model part.") - slave_model_part_name = settings["slave_model_part_name"].GetString() - self.slave_model_part = model.GetModelPart(slave_model_part_name) - - # Get the master model part - if not settings["master_model_part_name"].GetString(): - raise Exception("\'master_model_part_name\' not provided. Please specify the master model part.") - master_model_part_name = settings["master_model_part_name"].GetString() - self.master_model_part = model.GetModelPart(master_model_part_name) - - # Search radius for the MasterSlaveConstraints - self.search_radius = settings["search_radius"].GetDouble() - - # Minimum number of neighboring nodes retrieved from the search radius - self.minimum_number_of_neighbouring_nodes = settings["minimum_number_of_neighbouring_nodes"].GetInt() - - # Apply MasterSlaveConstraints at each time step (True) or only once (False) - self.reform_constraints_at_each_step = settings["reform_constraints_at_each_step"].GetBool() - - # Retrieve and check if variables exist - variable_names = settings["variable_names"].GetStringArray() - if len(variable_names) == 0: - err_msg = "The variable names need to be specified by the user in the \'variable_names\' string array." - raise Exception(err_msg) - if any(variable_names.count(var_name) > 1 for var_name in variable_names): - err_msg = "There are repeated variables in the \'variable_names\' string array." - raise Exception(err_msg) - variable_names.sort() - - self.variables_list = [] # Initialize the list of variables - - for var_name in variable_names: - # Check if the variable exists in KratosGlobals - if not KM.KratosGlobals.HasVariable(var_name): - err_msg = "\'{}\' variable in \'variable_names\' is not in KratosGlobals. Please check the provided value.".format(var_name) - - var_type = KM.KratosGlobals.GetVariableType(var_name) # Get the type of the variable - - # Check the variable type and add it to the variables_list accordingly - if var_type == "Array": - domain_size = self.computing_model_part.ProcessInfo[KM.DOMAIN_SIZE] # Get the domain size from the ProcessInfo - component_suffixes = ["_X", "_Y", "_Z"] # Suffixes for the components of the array variable - for i in range(domain_size): - var_name_with_suffix = f"{var_name}{component_suffixes[i]}" # Append the component suffix to the variable name - self.variables_list.append(KM.KratosGlobals.GetVariable(var_name_with_suffix)) # Add the variable to the list - elif var_type == "Double": - self.variables_list.append(KM.KratosGlobals.GetVariable(var_name)) # Add the variable to the list - else: - raise Exception("Variable " + var_name + " not compatible") # Raise an exception for incompatible variable types - - # Initialization of ghost nodes related attributes - self.model.CreateModelPart("AuxiliarMasterModelPart") - self.aux_master_model_part = self.model.GetModelPart("AuxiliarMasterModelPart") - self.max_node_id = self.computing_model_part.NumberOfNodes() - - def ExecuteInitialize(self): - """ This method is executed at the begining to initialize the process - - Keyword arguments: - self -- It signifies an instance of a class. - """ - - self.__CreateGhostNodes(self.master_model_part.Nodes) - - # Initialize MasterSlaveConstraints to neighbours utility - self.assign_mscs_utility = KM.AssignMasterSlaveConstraintsToNeighboursUtility(self.master_model_part.Nodes) - - # The user may only need to set up the MasterSlaveConstraints only once - if not self.reform_constraints_at_each_step: - for variable in self.variables_list: - self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes,self.search_radius,self.computing_model_part, variable, self.minimum_number_of_neighbouring_nodes) - - - def ExecuteInitializeSolutionStep(self): - """ This method is executed in order to initialize the current step - - Keyword arguments: - self -- It signifies an instance of a class. - """ - # If the user want the mscs to be updated at each time step, this is usefull for moving meshes. - if self.reform_constraints_at_each_step: - self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes,self.search_radius,self.computing_model_part, self.variables_list, self.minimum_number_of_neighbouring_nodes) - - def ExecuteFinalizeSolutionStep(self): - """ This method is executed in order to finalize the current step - - Keyword arguments: - self -- It signifies an instance of a class. - """ - # If MasterSlaveConstraints are updated every time step, these are to me removed before being re-assigned. - if self.reform_constraints_at_each_step: - self.__RemoveConstraints() - - def __RemoveConstraints(self): - #Remove master-slave constraints - KM.VariableUtils().SetFlag(KM.TO_ERASE, True, self.computing_model_part.MasterSlaveConstraints) - self.computing_model_part.RemoveMasterSlaveConstraintsFromAllLevels(KM.TO_ERASE) - - def __CreateGhostNodes(self, nodes): - """Create ghost nodes for each node.""" - ghost_to_original_mapping = {} - rotation_angles = [0, 60, 120, 180, 240, 300] - - self.max_node_id = self.main_model_part.NumberOfNodes() - - for angle in rotation_angles: - # Create ghost nodes - for original_node in nodes: - ghost_node = self.__RotateNode(original_node, angle) - ghost_to_original_mapping[ghost_node] = original_node - self.aux_master_model_part.AddNode(ghost_node) - - return ghost_to_original_mapping - - def __RotateNode(self, original_node, angle_degrees): - """Rotate a node by a given angle (in degrees) about the origin.""" - angle_radians = math.radians(angle_degrees) - - # Apply 2D rotation matrix - new_x = original_node.X * math.cos(angle_radians) - original_node.Y * math.sin(angle_radians) - new_y = original_node.X * math.sin(angle_radians) + original_node.Y * math.cos(angle_radians) - - # Generate a new node ID based on the maximum node ID - new_node_id = self.max_node_id + 1 - self.max_node_id += 1 - - # Create the new node with the ID and coordinates - return KM.Node(new_node_id, new_x, new_y, 0.0) \ No newline at end of file diff --git a/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py b/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py deleted file mode 100644 index 54f4824552e3..000000000000 --- a/kratos/python_scripts/assign_rotational_symmetry_master_slave_constraints_process.py +++ /dev/null @@ -1,182 +0,0 @@ -import KratosMultiphysics as KM -import numpy as np -import math - -def Factory(settings, model): - if not isinstance(settings, KM.Parameters): - raise Exception("expected input shall be a Parameters object, encapsulating a json string") - return AssignRotationalSymmetryMasterSlaveConstraintsProcess(model, settings["Parameters"]) - -## All the processes python should be derived from "Process" -class AssignRotationalSymmetryMasterSlaveConstraintsProcess(KM.Process): - """The process facilitates the discovery of neighboring nodes in a - master model part within a designated radius for each node in the - slave model part. Following this, it establishes a master-slave constraint - and calculates its weights using a spatial function that employs radial basis functions. - - Public member variables: - model -- the container of the different model parts. - settings -- Kratos parameters containing process settings. - """ - - def __init__(self, model, settings): - """ The default constructor of the class - - Keyword arguments: - self -- It signifies an instance of a class. - Model -- the container of the different model parts. - settings -- Kratos parameters containing process settings. - """ - KM.Process.__init__(self) # calling the baseclass constructor - - default_settings = KM.Parameters("""{ - "model_part_name": "", - "slave_model_part_name": "", - "master_model_part_name": "", - "variable_names": [], - "search_radius": 1.0, - "minimum_number_of_neighbouring_nodes": 3 - }""") - - # Add missing settings that the user did not provide but that - # are necessary for this process - settings.ValidateAndAssignDefaults(default_settings) - - # Get the model part on which the MasterSlaveConstraints are going to be applied - if not settings["model_part_name"].GetString(): - raise Exception("\'model_part_name\' not provided. Please specify the model part to apply to MasterSlaveConstraints to.") - model_part_name = settings["model_part_name"].GetString() #MasterSlaveConstraints are applied to computing model part - self.computing_model_part = model.GetModelPart(model_part_name) - - # Get the slave model part - if not settings["slave_model_part_name"].GetString(): - raise Exception("\'slave_model_part_name\' not provided. Please specify the slave model part.") - slave_model_part_name = settings["slave_model_part_name"].GetString() - self.slave_model_part = model.GetModelPart(slave_model_part_name) - - # Get the master model part - if not settings["master_model_part_name"].GetString(): - raise Exception("\'master_model_part_name\' not provided. Please specify the master model part.") - master_model_part_name = settings["master_model_part_name"].GetString() - self.master_model_part = model.GetModelPart(master_model_part_name) - - # Search radius for the MasterSlaveConstraints - self.search_radius = settings["search_radius"].GetDouble() - - # Minimum number of neighboring nodes retrieved from the search radius - self.minimum_number_of_neighbouring_nodes = settings["minimum_number_of_neighbouring_nodes"].GetInt() - - # Retrieve and check if variables exist - variable_names = settings["variable_names"].GetStringArray() - if len(variable_names) == 0: - err_msg = "The variable names need to be specified by the user in the \'variable_names\' string array." - raise Exception(err_msg) - if any(variable_names.count(var_name) > 1 for var_name in variable_names): - err_msg = "There are repeated variables in the \'variable_names\' string array." - raise Exception(err_msg) - variable_names.sort() - - self.variables_list = [] # Initialize the list of variables - - for var_name in variable_names: - # Check if the variable exists in KratosGlobals - if not KM.KratosGlobals.HasVariable(var_name): - err_msg = "\'{}\' variable in \'variable_names\' is not in KratosGlobals. Please check the provided value.".format(var_name) - - var_type = KM.KratosGlobals.GetVariableType(var_name) # Get the type of the variable - - # Check the variable type and add it to the variables_list accordingly - if var_type == "Array": - domain_size = self.computing_model_part.ProcessInfo[KM.DOMAIN_SIZE] # Get the domain size from the ProcessInfo - component_suffixes = ["_X", "_Y", "_Z"] # Suffixes for the components of the array variable - for i in range(domain_size): - var_name_with_suffix = f"{var_name}{component_suffixes[i]}" # Append the component suffix to the variable name - self.variables_list.append(KM.KratosGlobals.GetVariable(var_name_with_suffix)) # Add the variable to the list - elif var_type == "Double": - self.variables_list.append(KM.KratosGlobals.GetVariable(var_name)) # Add the variable to the list - else: - raise Exception("Variable " + var_name + " not compatible") # Raise an exception for incompatible variable types - - # Initialization of ghost nodes related attributes - model.CreateModelPart("AuxiliarMasterModelPart") - self.aux_master_model_part = model.GetModelPart("AuxiliarMasterModelPart") - self.max_node_id = self.computing_model_part.NumberOfNodes() - - def ExecuteInitialize(self): - """ This method is executed at the begining to initialize the process - - Keyword arguments: - self -- It signifies an instance of a class. - """ - - # Create ghost nodes in the master model part - self.__CreateGhostNodes(self.master_model_part.Nodes) - - # Initialize MasterSlaveConstraints to neighbours utility with the master model part nodes (which now includes ghost nodes) - self.assign_mscs_utility = KM.AssignMasterSlaveConstraintsToNeighboursUtility(self.aux_master_model_part.Nodes) - - def ExecuteInitializeSolutionStep(self): - """ This method is executed in order to initialize the current step - - Keyword arguments: - self -- It signifies an instance of a class. - """ - # Assign dynamically the master-slave constraints - # self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes, self.search_radius, self.computing_model_part, self.variables_list, self.minimum_number_of_neighbouring_nodes) - self.assign_mscs_utility.AssignRotationalMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes, self.search_radius, self.computing_model_part, self.variables_list, self.minimum_number_of_neighbouring_nodes, self.ghost_to_original_mapping) - - def ExecuteFinalizeSolutionStep(self): - """ This method is executed in order to finalize the current step - - Keyword arguments: - self -- It signifies an instance of a class. - """ - # Since MasterSlaveConstraints are updated every time step, they need to be removed before being re-assigned. - self.__RemoveConstraints() - - def __RemoveConstraints(self): - #Remove master-slave constraints - KM.VariableUtils().SetFlag(KM.TO_ERASE, True, self.computing_model_part.MasterSlaveConstraints) - self.computing_model_part.RemoveMasterSlaveConstraintsFromAllLevels(KM.TO_ERASE) - - def __CreateGhostNodes(self, nodes): - """Create ghost nodes for each node by rotating them about the origin. - - Keyword arguments: - nodes -- The nodes to be rotated to create the ghost nodes. - """ - self.ghost_to_original_mapping = {} # Initialize the mapping as an attribute of the class - rotation_angles = [0, -60, -120, -180, -240, -300] - self.max_node_id = self.computing_model_part.NumberOfNodes() - - for angle in rotation_angles: - # Create ghost nodes - for original_node in nodes: - ghost_node = self.__RotateNode(original_node, angle) - # Store the mapping using node IDs and the rotation angle - self.ghost_to_original_mapping[ghost_node.Id] = (original_node.Id, np.pi*angle/180) - if original_node.Id==922: - print(angle) - print(ghost_node.Id) - self.aux_master_model_part.AddNode(ghost_node) - - - def __RotateNode(self, original_node, angle_degrees): - """Rotate a node by a given angle (in degrees) about the origin. - - Keyword arguments: - original_node -- The node to be rotated. - angle_degrees -- The angle in degrees by which the node is to be rotated. - """ - angle_radians = math.radians(angle_degrees) - - # Apply 2D rotation matrix - new_x = original_node.X * math.cos(angle_radians) - original_node.Y * math.sin(angle_radians) - new_y = original_node.X * math.sin(angle_radians) + original_node.Y * math.cos(angle_radians) - - # Generate a new node ID based on the maximum node ID - new_node_id = self.max_node_id + 1 - self.max_node_id += 1 - - # Create the new node with the ID and coordinates - return KM.Node(new_node_id, new_x, new_y, 0.0) \ No newline at end of file diff --git a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp index 6c2a92ad7f0e..dbfd9b00f0c1 100644 --- a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp +++ b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp @@ -131,7 +131,7 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::GetDofsAndCoordinatesForCl // Loop over the cloud of nodes for (std::size_t i = 0; i < CloudOfNodesArray.size(); i++) { - auto p_slave_node = CloudOfNodesArray[i];//////CHANGE SLAVE NODE NAME + auto p_slave_node = CloudOfNodesArray[i]; // Loop over the variable list for (std::size_t j = 0; j < rVariableList.size(); j++) { @@ -154,99 +154,6 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::GetDofsAndCoordinatesForCl KRATOS_CATCH(""); } -void AssignMasterSlaveConstraintsToNeighboursUtility::GetDofsAndCoordinatesForCloudOfNodesFlattened( - const ResultNodesContainerType& CloudOfNodesArray, - const std::vector>>& rVariableList, - DofPointerVectorType& rCloudOfDofs, - Matrix& rCloudOfNodesCoordinates -) -{ - KRATOS_TRY; - - // Resize rCloudOfDofs to accommodate all DOFs for all variables - rCloudOfDofs.resize(rVariableList.size() * CloudOfNodesArray.size()); - - // Resize rCloudOfNodesCoordinates to match the size of the cloud of nodes - rCloudOfNodesCoordinates.resize(CloudOfNodesArray.size(), 3); - - // Loop over the cloud of nodes - for (std::size_t i = 0; i < CloudOfNodesArray.size(); i++) { - auto p_slave_node = CloudOfNodesArray[i]; - - // Loop over the variable list - for (std::size_t j = 0; j < rVariableList.size(); j++) { - const auto& r_variable = rVariableList[j]; - - // Check if the slave node has the required DOF for the current variable - if (p_slave_node->HasDofFor(r_variable.get())) { - std::size_t counter = i * rVariableList.size() + j; - rCloudOfDofs[counter] = p_slave_node->pGetDof(r_variable.get()); - } else { - KRATOS_ERROR << "The node with ID " << p_slave_node->Id() << " does not have the required DOF for the variable " << r_variable.get().Name() << std::endl; - } - } - - // Assign the coordinates of the slave node to the corresponding row in rCloudOfNodesCoordinates - noalias(row(rCloudOfNodesCoordinates, i)) = p_slave_node->Coordinates(); - } - - KRATOS_CATCH(""); -} - -struct CustomComparator { - static constexpr double TOLERANCE = 1e-6; - - bool operator()(const std::tuple& a, const std::tuple& b) const { - return std::fabs(std::get<0>(a) - std::get<0>(b)) > TOLERANCE || - std::fabs(std::get<1>(a) - std::get<1>(b)) > TOLERANCE || - std::fabs(std::get<2>(a) - std::get<2>(b)) > TOLERANCE; - } -}; - -void AssignMasterSlaveConstraintsToNeighboursUtility::FilterUniqueNodesForRBF( - ResultNodesContainerType& CloudOfNodesArray, - Matrix& rCloudOfNodesCoordinates -) { - KRATOS_TRY; - - // Use a set with a custom comparator to keep track of unique coordinates - std::set, CustomComparator> unique_coords; - - // Temporary vector to store unique nodes - std::vector temp_nodes; - - // Loop over the cloud of nodes - for (auto& p_node : CloudOfNodesArray) { - auto coords = p_node->Coordinates(); - std::tuple current_coords(coords[0], coords[1], coords[2]); - - // Check if the coordinates are unique - if (unique_coords.find(current_coords) == unique_coords.end()) { - unique_coords.insert(current_coords); - temp_nodes.push_back(p_node); - } - } - - // Clear the original CloudOfNodesArray and repopulate with unique nodes - CloudOfNodesArray.clear(); - for (auto& p_node : temp_nodes) { - CloudOfNodesArray.push_back(p_node); - } - - // Resize rCloudOfNodesCoordinates to match the size of unique nodes - rCloudOfNodesCoordinates.resize(temp_nodes.size(), 3); - - // Fill rCloudOfNodesCoordinates with unique coordinates - for (std::size_t i = 0; i < temp_nodes.size(); i++) { - auto coords = temp_nodes[i]->Coordinates(); - for (std::size_t j = 0; j < 3; j++) { - rCloudOfNodesCoordinates(i, j) = coords[j]; - } - } - - KRATOS_CATCH(""); -} - void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstraintsToNodes( NodesContainerType pSlaveNodes, double const Radius, @@ -309,12 +216,6 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstrain constant_vector.resize(r_n_container.size()); IndexType it = i.fetch_add(1) * rVariableList.size(); // Atomically increment the counter and get the previous value - if (p_slave_node->Id() == 929){ - KRATOS_WATCH(r_cloud_of_nodes) - KRATOS_WATCH(r_n_container) - KRATOS_WATCH(shape_matrix) - } - for (std::size_t j = 0; j < rVariableList.size(); ++j) { concurrent_constraints.enqueue(r_clone_constraint.Create(prev_num_mscs + it + j + 1, const_cast(r_cloud_of_dofs[j]), r_slave_dofs[j], shape_matrix, constant_vector)); } @@ -335,288 +236,6 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::AssignMasterSlaveConstrain KRATOS_CATCH(""); } -void AssignMasterSlaveConstraintsToNeighboursUtility::AssignRotationalMasterSlaveConstraintsToNodes( - NodesContainerType pSlaveNodes, - double const Radius, - ModelPart& rComputingModelPart, - const std::vector>>& rVariableList, - double const MinNumOfNeighNodes, - const std::unordered_map>& GhostToOriginalMapping -) -{ - KRATOS_TRY; - - BuiltinTimer build_and_assign_mscs; - - thread_local std::vector r_cloud_of_dofs; - thread_local DofPointerVectorType r_cloud_of_master_dofs; - thread_local std::vector r_slave_dofs; - thread_local Matrix r_cloud_of_nodes_coordinates; - thread_local Matrix r_cloud_of_master_nodes_coordinates; - thread_local array_1d r_slave_coordinates; - thread_local Vector r_n_container; - thread_local Vector constant_vector; - thread_local Matrix shape_matrix_x, shape_matrix_y, shape_matrix_z; - - int prev_num_mscs = rComputingModelPart.NumberOfMasterSlaveConstraints(); - - KRATOS_WARNING_IF("AssignRotationalMasterSlaveConstraintsToNodes", prev_num_mscs > 0) << "Previous Master-Slave Constraints exist in the ModelPart. The new Constraints may interact with the existing ones." << std::endl; - - auto& p_nodes_array = pSlaveNodes.GetContainer(); - - const auto& r_clone_constraint = KratosComponents::Get("LinearMasterSlaveConstraint"); - - using MasterSlaveConstraintQueue = moodycamel::ConcurrentQueue; - MasterSlaveConstraintQueue concurrent_constraints; - - std::atomic i(0); - - block_for_each(p_nodes_array, [&](Node::Pointer& p_slave_node) { - - // Skip nodes with IDs 1174 or 1015 950 1111 - if (p_slave_node->Id() == 1174 || p_slave_node->Id() == 1015) { - return; // Skip the current iteration of the loop - } - - double local_radius = Radius; - ResultNodesContainerType r_cloud_of_nodes(0); - ResultNodesContainerType r_local_results(mMaxNumberOfNodes); - - while (r_cloud_of_nodes.size() < MinNumOfNeighNodes) { - r_cloud_of_nodes.clear(); - SearchNodesInRadiusForNode(p_slave_node, local_radius, r_cloud_of_nodes, r_local_results); - local_radius += Radius; - } - - GetDofsAndCoordinatesForSlaveNode(p_slave_node, rVariableList, r_slave_dofs, r_slave_coordinates); - FilterUniqueNodesForRBF(r_cloud_of_nodes, r_cloud_of_nodes_coordinates); - - ResultNodesContainerType r_cloud_of_master_nodes; - r_cloud_of_master_nodes.reserve(r_cloud_of_nodes.size()); - for (auto& ghost_node : r_cloud_of_nodes) { - auto it = GhostToOriginalMapping.find(ghost_node->Id()); - if (it != GhostToOriginalMapping.end()) { - IndexType master_node_id = it->second.first; // Extracting the original/master node ID from the mapping - Node::Pointer p_master_node = rComputingModelPart.pGetNode(master_node_id); - r_cloud_of_master_nodes.push_back(p_master_node); - } else { - KRATOS_ERROR << "Ghost Node ID " << ghost_node->Id() << " not found in GhostToOriginalMapping!" << std::endl; - } - } - - GetDofsAndCoordinatesForCloudOfNodesFlattened(r_cloud_of_master_nodes, rVariableList, r_cloud_of_master_dofs, r_cloud_of_master_nodes_coordinates); - - RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_slave_coordinates, r_n_container); - - Vector ListOfAngles; - ExtractAnglesFromMapping(r_cloud_of_nodes, GhostToOriginalMapping, ListOfAngles); - - ObtainShapeMatrix(shape_matrix_x, shape_matrix_y, shape_matrix_z, r_n_container, ListOfAngles); - constant_vector.resize(shape_matrix_x.size2()); - if (p_slave_node->Id() == 929){ - KRATOS_WATCH(r_cloud_of_nodes) - KRATOS_WATCH(r_cloud_of_nodes_coordinates) - KRATOS_WATCH(r_slave_coordinates) - KRATOS_WATCH(r_cloud_of_master_nodes) - KRATOS_WATCH(r_cloud_of_master_dofs) - KRATOS_WATCH(r_n_container) - KRATOS_WATCH(ListOfAngles) - KRATOS_WATCH(shape_matrix_x) - KRATOS_WATCH(shape_matrix_y) - KRATOS_WATCH(shape_matrix_z) - } - - IndexType it = i.fetch_add(1) * rVariableList.size(); - - std::vector shape_matrices = {shape_matrix_x, shape_matrix_y, shape_matrix_z}; - - for (std::size_t j = 0; j < rVariableList.size(); ++j) { - if (j >= shape_matrices.size()) { - KRATOS_ERROR << "Unexpected index in rVariableList. Only 3 DoFs (x, y, z) are supported." << std::endl; - } - concurrent_constraints.enqueue(r_clone_constraint.Create(prev_num_mscs + it + j + 1, r_cloud_of_master_dofs, r_slave_dofs[j], shape_matrices[j], constant_vector)); - } - - - }); - - ConstraintContainerType temp_constraints; - ModelPart::MasterSlaveConstraintType::Pointer p_constraint; - - while (concurrent_constraints.try_dequeue(p_constraint)) - temp_constraints.push_back(p_constraint); - - rComputingModelPart.AddMasterSlaveConstraints(temp_constraints.begin(), temp_constraints.end()); - - KRATOS_INFO("AssignRotationalMasterSlaveConstraintsToNodes") << "Build and Assign Rotational Master-Slave Constraints Time: " << build_and_assign_mscs.ElapsedSeconds() << std::endl; - - KRATOS_CATCH(""); -} - -void AssignMasterSlaveConstraintsToNeighboursUtility::ObtainShapeMatrix( - Matrix& shape_matrix_x, - Matrix& shape_matrix_y, - Matrix& shape_matrix_z, - Vector& rN_container, - Vector& ListOfAngles -) -{ - std::vector rotation_x_global; - std::vector rotation_y_global; - std::vector rotation_z_global; - std::vector rotation_x(3); - std::vector rotation_y(3); - std::vector rotation_z(3); - for (int i=0; i < static_cast(rN_container.size()); i++){ - //TODO: Any kind of rotation. x-axis, y-axis, z-axis. Make it a function. - rotation_x[0] = rN_container[i]*cos(ListOfAngles[i]); - rotation_x[1] = -rN_container[i]*sin(ListOfAngles[i]); - rotation_x[2] = rN_container[i]*0.0; - rotation_y[0] = rN_container[i]*sin(ListOfAngles[i]); - rotation_y[1] = rN_container[i]*cos(ListOfAngles[i]); - rotation_y[2] = rN_container[i]*0.0; - rotation_z[0] = rN_container[i]*0.0; - rotation_z[1] = rN_container[i]*0.0; - rotation_z[2] = rN_container[i]*1.0; - // noalias(row(shape_matrix,i)) = aux_vector_x; - rotation_x_global.insert(rotation_x_global.end(),rotation_x.begin(),rotation_x.end()); - rotation_y_global.insert(rotation_y_global.end(),rotation_y.begin(),rotation_y.end()); - rotation_z_global.insert(rotation_z_global.end(),rotation_z.begin(),rotation_z.end()); - } - shape_matrix_x.resize(1, rotation_x_global.size()); - shape_matrix_y.resize(1, rotation_y_global.size()); - shape_matrix_z.resize(1, rotation_z_global.size()); - for (int i=0; i(rotation_x_global.size());i++){ - shape_matrix_x(0,i) = rotation_x_global[i]; - shape_matrix_y(0,i) = rotation_y_global[i]; - shape_matrix_z(0,i) = rotation_z_global[i]; - } -} - -// void AssignMasterSlaveConstraintsToNeighboursUtility::AssignRotationalMasterSlaveConstraintsToNodes( -// NodesContainerType pSlaveNodes, -// double const Radius, -// ModelPart& rComputingModelPart, -// const std::vector>>& rVariableList, -// double const MinNumOfNeighNodes, -// const std::unordered_map>& GhostToOriginalMapping -// ) -// { -// KRATOS_TRY; - -// BuiltinTimer build_and_assign_mscs; - -// thread_local std::vector r_cloud_of_dofs; -// thread_local std::vector r_slave_dofs; -// thread_local Matrix r_cloud_of_nodes_coordinates; -// thread_local array_1d r_slave_coordinates; -// thread_local Vector r_n_container; -// thread_local Vector constant_vector; -// thread_local Matrix shape_matrix; - -// int prev_num_mscs = rComputingModelPart.NumberOfMasterSlaveConstraints(); - -// KRATOS_WARNING_IF("AssignRotationalMasterSlaveConstraintsToNodes", prev_num_mscs > 0) << "Previous Master-Slave Constraints exist in the ModelPart. The new Constraints may interact with the existing ones." << std::endl; - -// auto& p_nodes_array = pSlaveNodes.GetContainer(); - -// const auto& r_clone_constraint = KratosComponents::Get("LinearMasterSlaveConstraint"); - -// using MasterSlaveConstraintQueue = moodycamel::ConcurrentQueue; -// MasterSlaveConstraintQueue concurrent_constraints; - -// std::atomic i(0); - -// block_for_each(p_nodes_array, [&](Node::Pointer& p_slave_node) { - -// double local_radius = Radius; -// ResultNodesContainerType r_cloud_of_nodes(0); -// ResultNodesContainerType r_local_results(mMaxNumberOfNodes); - -// while (r_cloud_of_nodes.size() < MinNumOfNeighNodes) { -// r_cloud_of_nodes.clear(); -// SearchNodesInRadiusForNode(p_slave_node, local_radius, r_cloud_of_nodes, r_local_results); -// local_radius += Radius; -// } - -// GetDofsAndCoordinatesForSlaveNode(p_slave_node, rVariableList, r_slave_dofs, r_slave_coordinates); -// GetDofsAndCoordinatesForCloudOfNodes(r_cloud_of_nodes, rVariableList, r_cloud_of_dofs, r_cloud_of_nodes_coordinates); - -// RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_slave_coordinates, r_n_container); - -// shape_matrix.resize(1, r_n_container.size()); -// noalias(row(shape_matrix, 0)) = r_n_container; -// constant_vector.resize(r_n_container.size()); - -// // Check if the current slave node is a ghost node and adjust for rotation if necessary -// if (GhostToOriginalMapping.find(p_slave_node->Id()) != GhostToOriginalMapping.end()) { -// auto mapping_data = GhostToOriginalMapping.at(p_slave_node->Id()); -// double rotation_angle = mapping_data.second; -// AdjustConstraintForRotation(rotation_angle, shape_matrix); -// } - -// IndexType it = i.fetch_add(1) * rVariableList.size(); - -// for (std::size_t j = 0; j < rVariableList.size(); ++j) { -// concurrent_constraints.enqueue(r_clone_constraint.Create(prev_num_mscs + it + j + 1, const_cast(r_cloud_of_dofs[j]), r_slave_dofs[j], shape_matrix, constant_vector)); -// } - -// }); - -// ConstraintContainerType temp_constraints; -// ModelPart::MasterSlaveConstraintType::Pointer p_constraint; - -// while (concurrent_constraints.try_dequeue(p_constraint)) -// temp_constraints.push_back(p_constraint); - -// rComputingModelPart.AddMasterSlaveConstraints(temp_constraints.begin(), temp_constraints.end()); - -// KRATOS_INFO("AssignRotationalMasterSlaveConstraintsToNodes") << "Build and Assign Rotational Master-Slave Constraints Time: " << build_and_assign_mscs.ElapsedSeconds() << std::endl; - -// KRATOS_CATCH(""); -// } - -void AssignMasterSlaveConstraintsToNeighboursUtility::AdjustConstraintForRotation(double rotation_angle, Matrix& shape_matrix) -{ - KRATOS_TRY; - - double cos_theta = std::cos(rotation_angle * M_PI / 180.0); // Convert degrees to radians - double sin_theta = std::sin(rotation_angle * M_PI / 180.0); - - for (std::size_t i = 0; i < shape_matrix.size1(); ++i) { - double original_w1 = shape_matrix(i, 0); - double original_w2 = shape_matrix(i, 1); - - shape_matrix(i, 0) = original_w1 * cos_theta - original_w2 * sin_theta; - shape_matrix(i, 1) = original_w1 * sin_theta + original_w2 * cos_theta; - } - - KRATOS_CATCH(""); -} - -void AssignMasterSlaveConstraintsToNeighboursUtility::ExtractAnglesFromMapping( - const ResultNodesContainerType& Results, - const std::unordered_map>& GhostToOriginalMapping, - Vector& ListOfAngles -) -{ - ListOfAngles.resize(Results.size(), false); - for(int i = 0; i < static_cast(Results.size()); i++) - { - auto it = GhostToOriginalMapping.find(Results[i]->Id()); - if(it != GhostToOriginalMapping.end()) - { - ListOfAngles[i] = it->second.second; // Extracting the angle from the mapping - } - else - { - KRATOS_ERROR << "Node ID " << Results[i]->Id() << " not found in GhostToOriginalMapping!" << std::endl; - } - } -} - -} // namespace Kratos. - void AssignMasterSlaveConstraintsToNeighboursUtility::InterpolateVelocityToNodes( NodesContainerType pDestinationNodes, double const Radius, diff --git a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h index 754e74aa4ccd..e4c66cbd225b 100644 --- a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h +++ b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.h @@ -137,29 +137,6 @@ class KRATOS_API(KRATOS_CORE) AssignMasterSlaveConstraintsToNeighboursUtility Matrix& rCloudOfNodesCoordinates ); - /** - * @brief Removes duplicate nodes based on their coordinates. - * - * This function is intended to filter out nodes with duplicate coordinates - * to avoid errors when using Radial Basis Functions (RBF) that can't handle - * repetitive coordinates. The function updates the provided CloudOfNodesArray - * to only contain unique nodes and also returns a matrix of their coordinates. - * - * @param CloudOfNodesArray A reference to the container of nodes. This will be updated to only contain unique nodes. - * @param rCloudOfNodesCoordinates A reference to a matrix that will be filled with the coordinates of the unique nodes. - */ - void FilterUniqueNodesForRBF( - ResultNodesContainerType& CloudOfNodesArray, - Matrix& rCloudOfNodesCoordinates - ); - - void GetDofsAndCoordinatesForCloudOfNodesFlattened( - const ResultNodesContainerType& CloudOfNodesArray, - const std::vector>>& rVariableList, - DofPointerVectorType& rCloudOfDofs, - Matrix& rCloudOfNodesCoordinates - ); - /** * @brief Assign Master-Slave Constraints to a set of Nodes. * @details Assign Matser-Slave Constraints to a set of Nodes given a radius of influence @@ -178,57 +155,12 @@ class KRATOS_API(KRATOS_CORE) AssignMasterSlaveConstraintsToNeighboursUtility double const MinNumOfNeighNodes ); - /** - * @brief Assign Rotational Master-Slave Constraints to a set of Nodes. - * @details This function assigns Master-Slave Constraints to a set of nodes considering rotational effects. - * It takes into account the rotation of the nodes and adjusts the constraints accordingly. The function - * uses MLS or RBF shape functions within a given radius of influence to determine the constraints. - * Additionally, it uses a mapping between ghost nodes and their original counterparts to account for the rotation. - * @param pSlaveNodes Nodes to which MasterSlaveConstraints are assigned. - * @param Radius Search radius for neighboring nodes. - * @param rComputingModelPart Model Part to which MasterSlaveConstraints are applied. - * @param rVariableList List of DOFs to assign the MasterSlaveConstraints. - * @param MinNumOfNeighNodes Minimum number of neighboring nodes required for the constraint. - * @param GhostToOriginalMapping Mapping between ghost nodes and their original counterparts, including the rotation angle. - */ - void AssignRotationalMasterSlaveConstraintsToNodes( - NodesContainerType pSlaveNodes, + void InterpolateVelocityToNodes( + NodesContainerType pDestinationNodes, double const Radius, - ModelPart& rComputingModelPart, const std::vector>>& rVariableList, - double const MinNumOfNeighNodes, - const std::unordered_map>& GhostToOriginalMapping - ); - /** - * @brief Adjusts the Master-Slave Constraint for rotational effects. - * @details This function modifies the shape functions matrix (weights) of the Master-Slave Constraint - * to account for the rotation of the nodes. Given a rotation angle, it applies the necessary - * trigonometric transformations to the shape functions to ensure that the constraints are correctly - * adjusted for the rotation. This is crucial when dealing with models that have rotational dynamics - * or when considering ghost nodes that represent rotated versions of the original nodes. - * @param rotation_angle The angle (in degrees) by which the node has been rotated. - * @param shape_matrix The shape functions matrix that needs adjustment for rotation. - */ - void AdjustConstraintForRotation( - double rotation_angle, - Matrix& shape_matrix - ); - - void ExtractAnglesFromMapping( - const ResultNodesContainerType& Results, - const std::unordered_map>& GhostToOriginalMapping, - Vector& ListOfAngles - ); - - void ObtainShapeMatrix( - Matrix& shape_matrix_x, - Matrix& shape_matrix_y, - Matrix& shape_matrix_z, - Vector& rN_container, - Vector& ListOfAngles - ); - - + double const MinNumOfNeighNodes + ); ///@} ///@name Input and output From f53cc51f4db453e0eade02e3bc8a5a021c668153 Mon Sep 17 00:00:00 2001 From: SADPR Date: Wed, 20 Dec 2023 12:24:46 +0100 Subject: [PATCH 20/20] Doing tests on interpolations. --- .../add_geometrical_utilities_to_python.cpp | 2 + ...slave_constraints_to_neighbours_process.py | 3 +- ...lave_constraints_to_neighbours_utility.cpp | 58 ++++++++++++++++--- 3 files changed, 54 insertions(+), 9 deletions(-) diff --git a/kratos/python/add_geometrical_utilities_to_python.cpp b/kratos/python/add_geometrical_utilities_to_python.cpp index f2913ad0c6c7..b574692f6644 100644 --- a/kratos/python/add_geometrical_utilities_to_python.cpp +++ b/kratos/python/add_geometrical_utilities_to_python.cpp @@ -439,6 +439,8 @@ void AddGeometricalUtilitiesToPython(pybind11::module &m) .def(py::init()) .def("AssignMasterSlaveConstraintsToNodes", [](AssignMasterSlaveConstraintsToNeighboursUtility& rAssignMasterSlaveConstraintsToNeighboursUtility, NodesContainerType pNodes, double const Radius, ModelPart& rComputingModelPart, const std::vector>>& rVariableList, double const MinNumOfNeighNodes){ return rAssignMasterSlaveConstraintsToNeighboursUtility.AssignMasterSlaveConstraintsToNodes(pNodes, Radius, rComputingModelPart, rVariableList, MinNumOfNeighNodes);}) + .def("InterpolateVelocityToNodes", [](AssignMasterSlaveConstraintsToNeighboursUtility& rAssignMasterSlaveConstraintsToNeighboursUtility, NodesContainerType pNodes, double const Radius, const std::vector>>& rVariableList, double const MinNumOfNeighNodes){ + return rAssignMasterSlaveConstraintsToNeighboursUtility.InterpolateVelocityToNodes(pNodes, Radius, rVariableList, MinNumOfNeighNodes);}) ; } diff --git a/kratos/python_scripts/assign_master_slave_constraints_to_neighbours_process.py b/kratos/python_scripts/assign_master_slave_constraints_to_neighbours_process.py index 91b5b2188884..f937a53d4c89 100644 --- a/kratos/python_scripts/assign_master_slave_constraints_to_neighbours_process.py +++ b/kratos/python_scripts/assign_master_slave_constraints_to_neighbours_process.py @@ -112,8 +112,7 @@ def ExecuteInitialize(self): # The user may only need to set up the MasterSlaveConstraints only once if not self.reform_constraints_at_each_step: - for variable in self.variables_list: - self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes,self.search_radius,self.computing_model_part, variable, self.minimum_number_of_neighbouring_nodes) + self.assign_mscs_utility.AssignMasterSlaveConstraintsToNodes(self.slave_model_part.Nodes,self.search_radius,self.computing_model_part, self.variables_list, self.minimum_number_of_neighbouring_nodes) def ExecuteInitializeSolutionStep(self): diff --git a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp index dbfd9b00f0c1..1339c2db6892 100644 --- a/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp +++ b/kratos/utilities/assign_master_slave_constraints_to_neighbours_utility.cpp @@ -272,22 +272,66 @@ void AssignMasterSlaveConstraintsToNeighboursUtility::InterpolateVelocityToNodes local_radius += Radius; } - // Get Dofs and Coordinates for the slave node and the cloud of nodes - GetDofsAndCoordinatesForSlaveNode(p_destination_node, rVariableList, r_destination_dofs, r_destination_coordinates); - GetDofsAndCoordinatesForCloudOfNodes(r_cloud_of_nodes, rVariableList, r_cloud_of_dofs, r_cloud_of_nodes_coordinates); + r_destination_coordinates = p_destination_node->Coordinates(); - // Calculate shape functions - RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_destination_coordinates, r_n_container); + // Resize rCloudOfNodesCoordinates to match the size of the cloud of nodes + r_cloud_of_nodes_coordinates.resize(r_cloud_of_nodes.size(), 3); - // IndexType it = i.fetch_add(1) * rVariableList.size(); // Atomically increment the counter and get the previous value + for (std::size_t i = 0; i < r_cloud_of_nodes.size(); i++) { + auto p_cloud_node = r_cloud_of_nodes[i]; + // Assign the coordinates of the slave node to the corresponding row in rCloudOfNodesCoordinates + noalias(row(r_cloud_of_nodes_coordinates, i)) = p_cloud_node->Coordinates(); + } + + // New part: Calculate weights using inverse distance + r_n_container.resize(r_cloud_of_nodes.size()); + double weight_sum = 0.0; + bool is_exact_overlap = false; + std::size_t overlap_index = 0; + + // Constants for modified IDW + const double epsilon = 1e-5; + const double p = 2.0; + + // First, check for exact overlap + for (std::size_t i = 0; i < r_cloud_of_nodes.size(); i++) { + auto dist = norm_2(row(r_cloud_of_nodes_coordinates, i) - r_destination_coordinates); + if (dist < 1e-12) { // A small threshold to consider as zero + is_exact_overlap = true; + overlap_index = i; + break; + } + } + + if (is_exact_overlap) { + // Set the weight of the overlapping node to 1 and others to 0 + r_n_container = ZeroVector(r_cloud_of_nodes.size()); + r_n_container[overlap_index] = 1.0; + } else { + // Calculate weights using modified IDW + r_n_container.resize(r_cloud_of_nodes.size()); + for (std::size_t i = 0; i < r_cloud_of_nodes.size(); i++) { + auto dist = norm_2(row(r_cloud_of_nodes_coordinates, i) - r_destination_coordinates); + r_n_container[i] = 1.0 / pow(dist + epsilon, p); // Modified IDW + weight_sum += r_n_container[i]; + } + // Normalize weights + for (auto& weight : r_n_container) { + weight /= weight_sum; + } + } + + // Calculate shape functions + // RBFShapeFunctionsUtility::CalculateShapeFunctions(r_cloud_of_nodes_coordinates, r_destination_coordinates, r_n_container); + for (std::size_t j = 0; j < rVariableList.size(); ++j) { double interpolated_variable_value = 0.0; for (std::size_t k = 0; k < r_cloud_of_nodes.size(); ++k){ auto r_current_node = r_cloud_of_nodes[k]; interpolated_variable_value += r_current_node->FastGetSolutionStepValue(rVariableList[j].get()) * r_n_container[k]; } - p_destination_node->FastGetSolutionStepValue(rVariableList[j].get(), interpolated_variable_value); + p_destination_node->FastGetSolutionStepValue(rVariableList[j].get()) = interpolated_variable_value; } }); KRATOS_CATCH("");