rotating_frame_process.py

import KratosMultiphysics as KM
import numpy as np
import time


def Factory(settings, model):
    if not isinstance(settings, KM.Parameters):
        raise Exception("expected input shall be a Parameters object, encapsulating a json string")
    return RotatingFrameProcess(model, settings["Parameters"])


class RotatingFrameProcess(KM.Process):
    """Process to apply rigid-body rotational kinematics to model parts.

    This process:
      - rotates a model part and assigns the corresponding MESH_DISPLACEMENT
      - assigns the corresponding rigid-body rotational VELOCITY to another model part

    The angular velocity is provided as the standard scalar input
    'angular_velocity_radians = f(x,y,z,t)' (or a constant number).
    The rotation is advanced by integrating the angular-velocity increment
    at each solution step.
    """

    def __init__(self, model, settings):
        KM.Process.__init__(self)

        default_settings = KM.Parameters("""{
            "rotating_frame_model_part_name": "",
            "rotating_object_model_part_name": "",
            "interval": [0.0, 1e30],
            "center_of_rotation": [0.0, 0.0, 0.0],
            "axis_of_rotation": [1.0, 0.0, 0.0],
            "angular_velocity_radians": 0.0,
            "echo_level": 0,
            "fix_mesh_displacement": false,
            "fix_velocity": false
        }""")

        # Allow a string function expression for 'angular_velocity_radians'.
        if settings.Has("angular_velocity_radians") and settings["angular_velocity_radians"].IsString():
            default_settings["angular_velocity_radians"].SetString("0.0")

        # Assign this here since it may normalize the "interval" before validation
        self.interval = KM.IntervalUtility(settings)

        settings.ValidateAndAssignDefaults(default_settings)

        self.model = model
        self.settings = settings

        if not settings["rotating_frame_model_part_name"].GetString():
            raise Exception("'rotating_frame_model_part_name' not provided. Please specify the model part to be rotated.")
        self.rotating_frame_model_part_name = settings["rotating_frame_model_part_name"].GetString()

        if not settings["rotating_object_model_part_name"].GetString():
            raise Exception("'rotating_object_model_part_name' not provided. Please specify the model part to which rotational velocity will be assigned.")
        self.rotating_object_model_part_name = settings["rotating_object_model_part_name"].GetString()

        self.rotating_frame_model_part = self.model.GetModelPart(self.rotating_frame_model_part_name)
        self.rotating_object_model_part = self.model.GetModelPart(self.rotating_object_model_part_name)
        self._interval_begin = settings["interval"][0].GetDouble()
        self._interval_end = settings["interval"][1].GetDouble()

        self.center_of_rotation = np.array(settings["center_of_rotation"].GetVector(), dtype=float)

        self.axis_of_rotation = np.array(settings["axis_of_rotation"].GetVector(), dtype=float)
        if self.axis_of_rotation.size == 0:
            raise Exception("The 'axis_of_rotation' vector is empty.")

        axis_norm = np.linalg.norm(self.axis_of_rotation)
        if np.isclose(axis_norm, 0.0):
            raise Exception("The 'axis_of_rotation' vector must have non-zero norm.")

        if not np.isclose(axis_norm, 1.0, rtol=1e-6):
            KM.Logger.PrintWarning(
                "RotatingFrameProcess",
                "The 'axis_of_rotation' vector is not a unit vector. It will be normalized."
            )
            self.axis_of_rotation /= axis_norm

        self.angular_velocity_function = self._CreateScalarFunction(
            settings["angular_velocity_radians"],
            "angular_velocity_radians"
        )

        self.echo_level = settings["echo_level"].GetInt()
        if self.echo_level < 0:
            raise Exception("The 'echo_level' parameter must be >= 0.")

        self.fix_mesh_displacement = settings["fix_mesh_displacement"].GetBool()
        self.fix_velocity = settings["fix_velocity"].GetBool()
        self._fixity_applied = False
        self._python_cache_initialized = False

    def ExecuteInitializeSolutionStep(self):
        current_time = self.rotating_frame_model_part.ProcessInfo[KM.TIME]
        delta_theta, omega, is_active = self._GetDeltaThetaAndOmegaFromFunction(current_time)

        if np.isclose(delta_theta, 0.0) and not is_active:
            # If we leave the active interval, release any fixity that this process applied.
            if self._fixity_applied:
                if self.fix_mesh_displacement:
                    self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT, False)
                if self.fix_velocity:
                    self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY, False)
                self._fixity_applied = False
            return

        if not self._python_cache_initialized:
            self._InitializePythonBackendCaches()

        t0 = time.perf_counter()
        if not np.isclose(delta_theta, 0.0):
            self._ApplyRotationAndMeshDisplacementPython(delta_theta)
        rot_time = time.perf_counter() - t0

        t0 = time.perf_counter()
        if is_active:
            self._AssignRotationalVelocityPython(omega)
        vel_time = time.perf_counter() - t0

        if self.echo_level > 0:
            KM.Logger.PrintInfo(
                "RotatingFrameProcess", f" rotation={rot_time*1.0e3:.3f} ms, velocity={vel_time*1.0e3:.3f} ms"
            )

        if is_active:
            if self.fix_mesh_displacement:
                self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT, True)

            if self.fix_velocity:
                self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY, True)

            self._fixity_applied = self.fix_mesh_displacement or self.fix_velocity
        elif self._fixity_applied:
            if self.fix_mesh_displacement:
                self._FixVectorVariable(self.rotating_frame_model_part, KM.MESH_DISPLACEMENT, False)
            if self.fix_velocity:
                self._FixVectorVariable(self.rotating_object_model_part, KM.VELOCITY, False)
            self._fixity_applied = False

    @staticmethod
    def _CreateScalarFunction(value, value_name):
        if value.IsNumber():
            function_string = str(value.GetDouble())
        elif value.IsString():
            function_string = value.GetString()
            if function_string == "":
                raise Exception(f"The '{value_name}' function string is empty.")
        else:
            raise Exception(
                f"The '{value_name}' parameter must be provided as a number or as a function string."
            )
        return KM.GenericFunctionUtility(function_string)

    def _EvaluateAngularVelocityFunction(self, time):
        x = self.center_of_rotation[0]
        y = self.center_of_rotation[1]
        z = self.center_of_rotation[2]
        return self.angular_velocity_function.CallFunction(x, y, z, time, 0.0, 0.0, 0.0)

    def _GetDeltaThetaAndOmegaFromFunction(self, time):
        omega = self._EvaluateAngularVelocityFunction(time)
        delta_time = self.rotating_frame_model_part.ProcessInfo[KM.DELTA_TIME]
        delta_theta = 0.0
        if delta_time > 0.0:
            step_begin = time - delta_time
            overlap_begin = max(step_begin, self._interval_begin)
            overlap_end = min(time, self._interval_end)

            if overlap_end > overlap_begin:
                omega_begin = self._EvaluateAngularVelocityFunction(overlap_begin)
                omega_end = self._EvaluateAngularVelocityFunction(overlap_end)
                delta_theta = 0.5 * (omega_begin + omega_end) * (overlap_end - overlap_begin)

        return delta_theta, omega, self.interval.IsInInterval(time)

    def _InitializePythonBackendCaches(self):
        self._rf_initial_pos_ta = KM.TensorAdaptors.NodePositionTensorAdaptor(
            self.rotating_frame_model_part.Nodes, KM.Configuration.Initial
        )
        self._rf_current_pos_ta = KM.TensorAdaptors.NodePositionTensorAdaptor(
            self.rotating_frame_model_part.Nodes, KM.Configuration.Current
        )
        self._rf_mesh_disp_ta = KM.TensorAdaptors.HistoricalVariableTensorAdaptor(
            self.rotating_frame_model_part.Nodes, KM.MESH_DISPLACEMENT
        )
        self._ro_current_pos_ta = KM.TensorAdaptors.NodePositionTensorAdaptor(
            self.rotating_object_model_part.Nodes, KM.Configuration.Current
        )
        self._ro_velocity_ta = KM.TensorAdaptors.HistoricalVariableTensorAdaptor(
            self.rotating_object_model_part.Nodes, KM.VELOCITY
        )

        # Initial frame coordinates are constant for a static mesh.
        self._rf_initial_pos_ta.CollectData()
        self._rf_initial_positions = self._rf_initial_pos_ta.data.copy()

        self._rf_rot_matrix = np.empty((3, 3), dtype=float)
        self._ro_relative_positions = np.empty(
            (self.rotating_object_model_part.NumberOfNodes(), 3),
            dtype=float
        )

        self._python_cache_initialized = True

    @staticmethod
    def _FixVectorVariable(model_part, variable, is_fixed):
        components_map = {
            KM.VELOCITY: (KM.VELOCITY_X, KM.VELOCITY_Y, KM.VELOCITY_Z),
            KM.MESH_DISPLACEMENT: (KM.MESH_DISPLACEMENT_X, KM.MESH_DISPLACEMENT_Y, KM.MESH_DISPLACEMENT_Z),
        }

        if variable not in components_map:
            raise Exception(f"Fixing variable '{variable.Name()}' is not supported by this process.")

        vu = KM.VariableUtils()
        for comp in components_map[variable]:
            vu.ApplyFixity(comp, is_fixed, model_part.Nodes)

    def _ApplyRotationAndMeshDisplacementPython(self, delta_theta):
        sin_half_theta = np.sin(delta_theta / 2.0)
        a = np.cos(delta_theta / 2.0)
        b = -self.axis_of_rotation[0] * sin_half_theta
        c = -self.axis_of_rotation[1] * sin_half_theta
        d = -self.axis_of_rotation[2] * sin_half_theta

        rot_matrix = self._rf_rot_matrix
        rot_matrix[0, 0] = a * a + b * b - c * c - d * d
        rot_matrix[0, 1] = 2.0 * (b * c - a * d)
        rot_matrix[0, 2] = 2.0 * (b * d + a * c)
        rot_matrix[1, 0] = 2.0 * (b * c + a * d)
        rot_matrix[1, 1] = a * a + c * c - b * b - d * d
        rot_matrix[1, 2] = 2.0 * (c * d - a * b)
        rot_matrix[2, 0] = 2.0 * (b * d - a * c)
        rot_matrix[2, 1] = 2.0 * (c * d + a * b)
        rot_matrix[2, 2] = a * a + d * d - b * b - c * c

        # Rotate current coordinates around the center with the incremental angle.
        self._rf_current_pos_ta.CollectData()
        current_positions = self._rf_current_pos_ta.data
        current_positions[:] = (current_positions - self.center_of_rotation) @ rot_matrix
        current_positions += self.center_of_rotation

        # Mesh displacement from reference configuration: u_mesh = x - x0
        self._rf_mesh_disp_ta.data[:] = current_positions - self._rf_initial_positions

        self._rf_mesh_disp_ta.StoreData()
        self._rf_current_pos_ta.StoreData()

    def _AssignRotationalVelocityPython(self, omega):
        angular_velocity_vector = omega * self.axis_of_rotation

        self._ro_current_pos_ta.CollectData()
        self._ro_relative_positions[:] = self._ro_current_pos_ta.data - self.center_of_rotation

        # v = omega x r, matching C++ MathUtils<double>::CrossProduct(omega, r)
        v = self._ro_velocity_ta.data
        r = self._ro_relative_positions
        v[:, 0] = angular_velocity_vector[1] * r[:, 2] - angular_velocity_vector[2] * r[:, 1]
        v[:, 1] = angular_velocity_vector[2] * r[:, 0] - angular_velocity_vector[0] * r[:, 2]
        v[:, 2] = angular_velocity_vector[0] * r[:, 1] - angular_velocity_vector[1] * r[:, 0]

        self._ro_velocity_ta.StoreData()