diff --git a/docs/continuum_mechanics/functions/doc_rotation.rst b/docs/continuum_mechanics/functions/doc_rotation.rst new file mode 100644 index 000000000..fc699f4fe --- /dev/null +++ b/docs/continuum_mechanics/functions/doc_rotation.rst @@ -0,0 +1,421 @@ +The Rotation Library +==================== + +The rotation library provides comprehensive tools for 3D rotations in continuum mechanics. +The ``Rotation`` class inherits from ``scipy.spatial.transform.Rotation``, so all scipy +features (batch operations, ``mean()``, ``Slerp``, ``RotationSpline``, etc.) are available +directly, while simcoon adds continuum-mechanics methods for Voigt-notation tensors. + +.. contents:: Contents + :local: + :depth: 2 + +Architecture +------------ + +The rotation system is organized in three layers: + +**C++ core** (``simcoon::Rotation`` in ``include/simcoon/Simulation/Maths/rotation.hpp``): + The full rotation class with quaternion internals, all factory methods, conversions, + algebra (composition, inverse, SLERP), and the Voigt mechanics operations. This is + what C++ code uses directly. + +**pybind11 bridge** (``_CppRotation``): + A minimal binding that exposes only what the Python layer needs from C++: a single + factory method (``from_quat``) and the nine mechanics operations (``apply_stress``, + ``apply_stiffness``, etc.). Everything else — creation, conversion, composition, + interpolation — is handled by scipy on the Python side. + +**Python class** (``simcoon.Rotation``): + Inherits from ``scipy.spatial.transform.Rotation``. Users get the full scipy API + for free (factory methods, conversions, batch operations, ``mean()``, ``Slerp``, + ``RotationSpline``). The mechanics methods delegate to C++ by converting the + quaternion on each call: + +.. code-block:: text + + r.apply_stiffness(L) + | + v + Python: self._to_cpp() --> _CppRotation.from_quat(self.as_quat()) + | ^ scipy provides the quaternion + v + C++: simcoon::Rotation::from_quat(q) --> .apply_stiffness(L) + | + v + Python: returns 6x6 numpy array + +The quaternion transfer (4 doubles) is negligible. This design avoids duplicating +scipy's well-tested rotation algebra while keeping the performance-critical Voigt +operations in C++. + +The Rotation Class +------------------ + +The ``Rotation`` class extends ``scipy.spatial.transform.Rotation`` with methods for +rotating stress, strain, stiffness, and compliance tensors. It uses unit quaternions +internally for numerical stability and efficient composition. + +Creating Rotations +~~~~~~~~~~~~~~~~~~ + +All scipy factory methods are available and return a ``simcoon.Rotation`` instance: + +.. code-block:: python + + import simcoon as smc + import numpy as np + + # Identity rotation (no rotation) + r = smc.Rotation.identity() + + # From Euler angles (scipy convention: uppercase = intrinsic) + r = smc.Rotation.from_euler('ZXZ', [psi, theta, phi]) + + # From Euler angles in degrees, extrinsic + r = smc.Rotation.from_euler('xyz', [roll, pitch, yaw], degrees=True) + + # From rotation matrix + R = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) # 90° around x + r = smc.Rotation.from_matrix(R) + + # From quaternion [qx, qy, qz, qw] (scalar-last) + q = np.array([0, 0, np.sin(np.pi/4), np.cos(np.pi/4)]) # 90° around z + r = smc.Rotation.from_quat(q) + + # From rotation vector (axis × angle) + rotvec = np.array([0, 0, np.pi/2]) # 90° around z + r = smc.Rotation.from_rotvec(rotvec) + + # From axis and angle (simcoon-specific) + r = smc.Rotation.from_axis_angle(np.pi/2, 3) # 90° around z (axis 3) + + # Random rotation (uniform distribution) + r = smc.Rotation.random() + +Converting Between Representations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + r = smc.Rotation.from_euler('ZXZ', [0.5, 0.3, 0.7]) + + # To rotation matrix + R = r.as_matrix() # 3×3 numpy array + + # To quaternion + q = r.as_quat() # [qx, qy, qz, qw] + + # To Euler angles (scipy convention) + angles = r.as_euler('ZXZ') # [psi, theta, phi] + angles_deg = r.as_euler('ZXZ', degrees=True) + + # To rotation vector + rotvec = r.as_rotvec() + rotvec_deg = r.as_rotvec(degrees=True) + + # To Voigt rotation matrices (simcoon-specific) + QS = r.as_voigt_stress_rotation() # 6×6 for stress + QE = r.as_voigt_strain_rotation() # 6×6 for strain + +Applying Rotations +~~~~~~~~~~~~~~~~~~ + +**To 3D vectors:** + +.. code-block:: python + + r = smc.Rotation.from_axis_angle(np.pi/2, 3) + + v = np.array([1.0, 0.0, 0.0]) + v_rot = r.apply(v) # [0, 1, 0] + + # Inverse rotation + v_back = r.apply(v_rot, inverse=True) # [1, 0, 0] + +**To 3×3 tensors:** + +.. code-block:: python + + # Rotate a tensor: R · T · R^T + T = np.eye(3) + T_rot = r.apply_tensor(T) + +**To Voigt notation tensors:** + +.. code-block:: python + + # Stress vector [σ11, σ22, σ33, σ12, σ13, σ23] + sigma = np.array([100.0, 50.0, 25.0, 10.0, 5.0, 2.0]) + sigma_rot = r.apply_stress(sigma) + + # Strain vector [ε11, ε22, ε33, 2ε12, 2ε13, 2ε23] + epsilon = np.array([0.01, -0.005, -0.005, 0.002, 0.001, 0.0]) + epsilon_rot = r.apply_strain(epsilon) + + # Stiffness matrix (6×6) + L = smc.L_iso([210e9, 0.3], "Enu") + L_rot = r.apply_stiffness(L) + + # Compliance matrix (6×6) + M = smc.M_iso([210e9, 0.3], "Enu") + M_rot = r.apply_compliance(M) + + # Strain concentration tensor (6×6): A' = QE * A * QS^T + A_global = r.apply_strain_concentration(A_local) + + # Stress concentration tensor (6×6): B' = QS * B * QE^T + B_global = r.apply_stress_concentration(B_local) + +Composing Rotations +~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + r1 = smc.Rotation.from_axis_angle(np.pi/4, 3) # 45° around z + r2 = smc.Rotation.from_axis_angle(np.pi/6, 1) # 30° around x + + # Compose: apply r1 first, then r2 + r_combined = r2 * r1 + + # This is equivalent to: + v = np.array([1.0, 0.0, 0.0]) + v1 = r2.apply(r1.apply(v)) + v2 = r_combined.apply(v) + # v1 == v2 + + # Inverse + r_inv = r1.inv() + r_identity = r1 * r_inv # Should be identity + +Interpolating Rotations +~~~~~~~~~~~~~~~~~~~~~~~ + +SLERP (Spherical Linear Interpolation) is available via scipy's ``Slerp`` class: + +.. code-block:: python + + from scipy.spatial.transform import Slerp + + r_start = smc.Rotation.identity() + r_end = smc.Rotation.from_euler('ZXZ', [np.pi/2, np.pi/4, 0]) + + key_rots = smc.Rotation.concatenate([r_start, r_end]) + slerp = Slerp([0, 1], key_rots) + + # Interpolate along path + for t in np.linspace(0, 1, 11): + r_t = slerp(t) + print(f"t={t:.1f}: angle={np.degrees(r_t.magnitude()):.1f}°") + + # Or use the convenience method + r_mid = r_start.slerp_to(r_end, 0.5) + +Utility Methods +~~~~~~~~~~~~~~~ + +.. code-block:: python + + r = smc.Rotation.from_euler('ZXZ', [0.5, 0.3, 0.7]) + + # Get rotation magnitude (angle) + angle = r.magnitude() # radians + angle_deg = np.degrees(r.magnitude()) + + # Check if identity + is_id = r.is_identity() + + # Check equality (accounts for quaternion sign ambiguity) + r2 = smc.Rotation.from_euler('ZXZ', [0.5, 0.3, 0.7]) + are_equal = r.equals(r2, tol=1e-12) + +Scipy Integration +~~~~~~~~~~~~~~~~~ + +Since ``simcoon.Rotation`` inherits from ``scipy.spatial.transform.Rotation``, +all scipy features work directly: + +.. code-block:: python + + from scipy.spatial.transform import Rotation as R + + # simcoon.Rotation IS a scipy Rotation + r = smc.Rotation.from_axis_angle(np.pi/4, 3) + isinstance(r, R) # True + + # All scipy methods work directly + r.as_mrp() # Modified Rodrigues parameters + r.as_davenport(...) # Davenport angles + + # Batch operations + rots = smc.Rotation.random(100) + rots.mean() + + # RotationSpline + from scipy.spatial.transform import RotationSpline + times = [0, 1, 2, 3] + key_rots = smc.Rotation.random(4) + spline = RotationSpline(times, key_rots) + r_interp = spline(1.5) + +**Passing plain scipy Rotations to simcoon:** + +If you already have a ``scipy.spatial.transform.Rotation`` (e.g. from another +library), you can upgrade it to a ``simcoon.Rotation`` with ``from_scipy()``: + +.. code-block:: python + + from scipy.spatial.transform import Rotation as R + + # Created somewhere else, e.g. by a third-party library + scipy_rot = R.from_euler('z', 45, degrees=True) + + # Upgrade to simcoon.Rotation to unlock mechanics methods + r = smc.Rotation.from_scipy(scipy_rot) + L_rot = r.apply_stiffness(L) # now available + +Active vs Passive Rotations +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The ``active`` parameter on Voigt methods (``apply_stress``, ``apply_strain``, +``apply_stiffness``, ``apply_compliance``, ``as_voigt_stress_rotation``, ``as_voigt_strain_rotation``) controls the +rotation convention: + +- **active=True** (default): **Alibi** rotation — rotates the physical object + while the coordinate system stays fixed. +- **active=False**: **Alias** rotation — rotates the coordinate system while the + object stays fixed. This is equivalent to the inverse active rotation. + +.. code-block:: python + + r = smc.Rotation.from_axis_angle(np.pi/4, 3) + + # Active: rotate the stress tensor itself + sigma_active = r.apply_stress(sigma, active=True) + + # Passive: express the same stress in a rotated coordinate system + sigma_passive = r.apply_stress(sigma, active=False) + +Input Validation +~~~~~~~~~~~~~~~~ + +All mechanics methods that accept NumPy arrays validate the input shape and raise +``ValueError`` with a descriptive message if the dimensions are wrong: + +.. code-block:: python + + r = smc.Rotation.from_axis_angle(np.pi/4, 3) + + r.apply_stress(np.array([1.0, 2.0, 3.0])) + # ValueError: sigma must have 6 elements, got 3 + + r.apply_stiffness(np.eye(3)) + # ValueError: L must have shape (6, 6), got (3, 3) + +Expected input shapes: + +- ``apply``: 1D array, 3 elements (or Nx3 array for batch) +- ``apply_tensor``: 2D array, shape (3, 3) +- ``apply_stress``, ``apply_strain``: 1D array, 6 elements +- ``apply_stiffness``, ``apply_compliance``: 2D array, shape (6, 6) + +.. note:: + + Quaternions ``q`` and ``-q`` represent the same rotation. The ``equals()`` + method accounts for this antipodal equivalence. + +Euler Angle Conventions +----------------------- + +Simcoon uses scipy's Euler angle conventions: + +- **Uppercase** letters (``'ZXZ'``, ``'XYZ'``, etc.) → intrinsic rotations (rotating frame) +- **Lowercase** letters (``'zxz'``, ``'xyz'``, etc.) → extrinsic rotations (fixed frame) + +**Proper Euler angles** (axis sequence where first = last): + +- ``ZXZ``, ``ZYZ``, ``XYX``, ``XZX``, ``YXY``, ``YZY`` + +**Tait-Bryan angles** (all axes different): + +- ``XYZ`` (roll-pitch-yaw), ``XZY``, ``YXZ``, ``YZX``, ``ZXY``, ``ZYX`` + +.. code-block:: python + + # Intrinsic ZXZ + r1 = smc.Rotation.from_euler('ZXZ', [psi, theta, phi]) + + # Extrinsic xyz (same physical rotation as intrinsic ZYX) + r2 = smc.Rotation.from_euler('xyz', [roll, pitch, yaw]) + +Practical Examples +------------------ + +Example 1: Rotating Orthotropic Material Properties +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import simcoon as smc + import numpy as np + + # Define orthotropic material in local coordinates + # [E1, E2, E3, nu12, nu13, nu23, G12, G13, G23] + props = [150e9, 10e9, 10e9, 0.3, 0.3, 0.45, 5e9, 5e9, 3.5e9] + L_local = smc.L_ortho(props, "EnuG") + + # Fiber orientation: 45° rotation around z-axis + r = smc.Rotation.from_axis_angle(np.pi/4, 3) + + # Rotate stiffness to global coordinates + L_global = r.apply_stiffness(L_local) + + print("Local stiffness L11:", L_local[0, 0] / 1e9, "GPa") + print("Global stiffness L11:", L_global[0, 0] / 1e9, "GPa") + +Example 2: Stress Transformation in a Rotated Element +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import simcoon as smc + import numpy as np + + # Stress in global coordinates (Voigt notation) + sigma_global = np.array([100.0, 50.0, 0.0, 25.0, 0.0, 0.0]) # MPa + + # Element orientation (Euler angles) + psi, theta, phi = np.radians([30, 45, 60]) + + # Create rotation (scipy convention: uppercase = intrinsic) + r = smc.Rotation.from_euler('ZXZ', [psi, theta, phi]) + + # Transform stress to local element coordinates + sigma_local = r.inv().apply_stress(sigma_global) + + print("Global stress:", sigma_global) + print("Local stress:", sigma_local) + +Example 3: Averaging Orientations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import simcoon as smc + import numpy as np + + # Generate random orientations + rotations = smc.Rotation.random(10) + + # Use scipy's built-in mean + r_mean = rotations.mean() + + print("Mean rotation angle:", np.degrees(r_mean.magnitude()), "degrees") + +See Also +-------- + +- :doc:`doc_constitutive` - Constitutive matrices (L_iso, L_ortho, etc.) +- :doc:`doc_transfer` - Voigt notation conversions +- :doc:`doc_kinematics` - Finite strain kinematics +- :doc:`/cpp_api/simulation/rotation` - Full C++ API reference diff --git a/docs/continuum_mechanics/functions/index.rst b/docs/continuum_mechanics/functions/index.rst index 137edc73c..96ec73a72 100755 --- a/docs/continuum_mechanics/functions/index.rst +++ b/docs/continuum_mechanics/functions/index.rst @@ -6,6 +6,7 @@ Continuum Mechanics generic functions :caption: Contents: doc_constitutive + doc_rotation doc_contimech doc_criteria doc_transfer diff --git a/docs/cpp_api/simulation/index.rst b/docs/cpp_api/simulation/index.rst index a9322a848..a0a8b4617 100644 --- a/docs/cpp_api/simulation/index.rst +++ b/docs/cpp_api/simulation/index.rst @@ -2,7 +2,7 @@ Simulation ========== -The Simulation module provides tools for running simulations, including solvers, +The Simulation module provides tools for running simulations, including solvers, identification algorithms, and mathematical utilities. .. toctree:: @@ -10,6 +10,7 @@ identification algorithms, and mathematical utilities. :caption: Submodules: maths + rotation solver identification phase @@ -19,7 +20,23 @@ Overview This module contains simulation and numerical tools: -- **Maths**: Mathematical utilities (rotation, random numbers, solvers) +- **Maths**: Mathematical utilities (random numbers, statistics, solvers) +- **Rotation**: Comprehensive 3D rotation tools with ``Rotation`` class and Voigt notation support - **Solver**: Material point simulation solvers - **Identification**: Parameter identification algorithms - **Phase**: Phase management and properties + +What's New in Simcoon 2.0 +========================= + +The **Rotation module** has been significantly enhanced with: + +- New ``Rotation`` class inspired by ``scipy.spatial.transform.Rotation`` +- Unit quaternion internal representation for numerical stability +- Support for multiple Euler angle conventions (zxz, zyz, xyz, etc.) +- Intrinsic and extrinsic rotation modes +- SLERP interpolation for smooth orientation transitions +- Direct Voigt notation operations (stress, strain, stiffness, compliance) +- Full Python bindings with NumPy integration + +See :doc:`rotation` for complete documentation. diff --git a/docs/cpp_api/simulation/maths.rst b/docs/cpp_api/simulation/maths.rst index 3607ecde6..5ca721e9a 100644 --- a/docs/cpp_api/simulation/maths.rst +++ b/docs/cpp_api/simulation/maths.rst @@ -5,6 +5,40 @@ Mathematical Functions Mathematical utility functions for rotations, random number generation, statistics, and numerical solving. -.. doxygengroup:: maths +.. toctree:: + :maxdepth: 2 + + rotation + +Overview +======== + +The Maths module provides essential mathematical utilities for continuum mechanics +simulations: + +- **Rotation**: Comprehensive 3D rotation tools including the new ``Rotation`` class + and Voigt notation transformations (see :doc:`rotation`) +- **Random**: Random number generation utilities +- **Statistics**: Statistical functions +- **Numerical Solvers**: Root finding and optimization + +Random Number Generation +======================== + +.. doxygengroup:: random + :project: simcoon + :content-only: + +Statistics +========== + +.. doxygengroup:: stats + :project: simcoon + :content-only: + +Numerical Solvers +================= + +.. doxygengroup:: solvers :project: simcoon :content-only: diff --git a/docs/cpp_api/simulation/rotation.rst b/docs/cpp_api/simulation/rotation.rst new file mode 100644 index 000000000..324f7b10d --- /dev/null +++ b/docs/cpp_api/simulation/rotation.rst @@ -0,0 +1,476 @@ +================== +The Rotation Module +================== + +The Rotation module provides comprehensive tools for 3D rotations in continuum mechanics +applications. It provides a modern **Rotation class** (inspired by ``scipy.spatial.transform.Rotation``) +for all rotation operations including tensor rotation in Voigt notation. + +.. contents:: Table of Contents + :local: + :depth: 2 + +Overview +======== + +Simcoon 2.0 introduces a new ``Rotation`` class that encapsulates rotation operations using +unit quaternions as the internal representation. This design provides: + +- **Numerical stability**: Avoids gimbal lock issues present in Euler angle representations +- **Efficient composition**: Quaternion multiplication is faster than matrix multiplication +- **Multiple representations**: Easy conversion between quaternions, matrices, Euler angles, and rotation vectors +- **Voigt notation support**: Direct application to stress/strain tensors and stiffness/compliance matrices + +The C++ free functions remain available for direct operations on matrices and vectors. + +The Rotation Class +================== + +The ``Rotation`` class uses unit quaternions in **scalar-last convention** ``[qx, qy, qz, qw]`` +as its internal representation. + +C++ API +------- + +.. code-block:: cpp + + #include + + using namespace simcoon; + + // Create rotations + Rotation r1 = Rotation::identity(); + Rotation r2 = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + Rotation r3 = Rotation::from_axis_angle(M_PI/4, 3); // 45° around z-axis + Rotation r4 = Rotation::from_matrix(R); + Rotation r5 = Rotation::random(); + + // Convert to different representations + arma::mat::fixed<3,3> R = r2.as_matrix(); + arma::vec::fixed<4> q = r2.as_quat(); + arma::vec::fixed<3> euler = r2.as_euler("zxz"); + arma::vec::fixed<3> rotvec = r2.as_rotvec(); + + // Apply to vectors and tensors + arma::vec::fixed<3> v_rot = r2.apply(v); + arma::vec::fixed<6> sigma_rot = r2.apply_stress(sigma); + arma::mat::fixed<6,6> L_rot = r2.apply_stiffness(L); + + // Compose rotations + Rotation r_combined = r2 * r3; + Rotation r_inv = r2.inv(); + + // Interpolation + Rotation r_mid = r2.slerp(r3, 0.5); + +Python API +---------- + +.. code-block:: python + + import simcoon as smc + import numpy as np + + # Create rotations (scipy factory methods + simcoon's from_axis_angle) + r1 = smc.Rotation.identity() + r2 = smc.Rotation.from_euler('ZXZ', [0.5, 0.3, 0.7]) + r3 = smc.Rotation.from_axis_angle(np.pi/4, 3) + r4 = smc.Rotation.from_matrix(R) + r5 = smc.Rotation.random() + + # Convert to different representations + R = r2.as_matrix() + q = r2.as_quat() + euler = r2.as_euler('ZXZ') + rotvec = r2.as_rotvec() + + # Apply to vectors and tensors + v_rot = r2.apply(v) + sigma_rot = r2.apply_stress(sigma) + L_rot = r2.apply_stiffness(L) + + # Compose rotations + r_combined = r2 * r3 + r_inv = r2.inv() + + # Interpolation + r_mid = r2.slerp_to(r3, 0.5) + +Factory Methods +--------------- + +.. list-table:: Factory Methods + :widths: 30 70 + :header-rows: 1 + + * - Method + - Description + * - ``identity()`` + - Create an identity rotation (no rotation) + * - ``from_quat(quat)`` + - Create from quaternion [qx, qy, qz, qw] (scalar-last) + * - ``from_matrix(R)`` + - Create from 3x3 rotation matrix + * - ``from_euler(seq, angles, degrees)`` + - Create from Euler angles (scipy convention, e.g. ``'ZXZ'``) + * - ``from_rotvec(rotvec, degrees)`` + - Create from rotation vector (axis × angle) + * - ``from_axis_angle(angle, axis, degrees)`` + - Create from angle and axis (1=x, 2=y, 3=z) + * - ``random()`` + - Create uniformly distributed random rotation + +Euler Angle Conventions +----------------------- + +The ``from_euler`` and ``as_euler`` methods support multiple conventions: + +**Proper Euler angles** (first and last axis are the same): + +- ``zxz``, ``zyz``, ``xyx``, ``xzx``, ``yxy``, ``yzy`` + +**Tait-Bryan angles** (all three axes are different): + +- ``xyz``, ``xzy``, ``yxz``, ``yzx``, ``zxy``, ``zyx`` + +**Parameters:** + +- ``intrinsic`` (bool): If ``true`` (default), rotations are about the axes of the rotating + coordinate system. If ``false``, rotations are about the fixed coordinate system (extrinsic). +- ``degrees`` (bool): If ``true``, angles are in degrees. Default is ``false`` (radians). + +.. code-block:: cpp + + // Intrinsic ZXZ Euler angles + Rotation r1 = Rotation::from_euler(psi, theta, phi, "zxz", true, false); + + // Extrinsic XYZ (aerospace convention) + Rotation r2 = Rotation::from_euler(roll, pitch, yaw, "xyz", false, true); + +Conversion Methods +------------------ + +.. list-table:: Conversion Methods + :widths: 30 70 + :header-rows: 1 + + * - Method + - Description + * - ``as_quat()`` + - Return quaternion [qx, qy, qz, qw] + * - ``as_matrix()`` + - Return 3×3 rotation matrix + * - ``as_euler(conv, intrinsic, degrees)`` + - Return Euler angles [psi, theta, phi] + * - ``as_rotvec(degrees)`` + - Return rotation vector (axis × angle) + * - ``as_voigt_stress_rotation(active)`` + - Return 6×6 stress rotation matrix + * - ``as_voigt_strain_rotation(active)`` + - Return 6×6 strain rotation matrix + +Apply Methods +------------- + +**3D Vectors and Tensors:** + +.. list-table:: Apply Methods for 3D Objects + :widths: 30 70 + :header-rows: 1 + + * - Method + - Description + * - ``apply(v, inverse)`` + - Rotate a 3D vector + * - ``apply_tensor(m, inverse)`` + - Rotate a 3×3 tensor: R·m·R\ :sup:`T` + +**Voigt Notation (Continuum Mechanics):** + +.. list-table:: Apply Methods for Voigt Notation + :widths: 30 70 + :header-rows: 1 + + * - Method + - Description + * - ``apply_stress(sigma, active)`` + - Rotate 6-component stress vector + * - ``apply_strain(epsilon, active)`` + - Rotate 6-component strain vector + * - ``apply_stiffness(L, active)`` + - Rotate 6×6 stiffness matrix: QS·L·QS\ :sup:`T` + * - ``apply_compliance(M, active)`` + - Rotate 6×6 compliance matrix: QE·M·QE\ :sup:`T` + * - ``apply_strain_concentration(A, active)`` + - Rotate 6×6 strain concentration tensor: QE·A·QS\ :sup:`T` + * - ``apply_stress_concentration(B, active)`` + - Rotate 6×6 stress concentration tensor: QS·B·QE\ :sup:`T` + +The ``active`` parameter controls whether the rotation is **active** (alibi, rotating the object) +or **passive** (alias, rotating the coordinate system). + +Operations +---------- + +.. list-table:: Rotation Operations + :widths: 30 70 + :header-rows: 1 + + * - Method + - Description + * - ``inv()`` + - Return inverse rotation + * - ``magnitude()`` + - Return rotation angle in radians + * - ``r1 * r2`` + - Compose rotations (apply r2 first, then r1) + * - ``r1 *= r2`` + - Compose in-place + * - ``slerp_to(other, t)`` + - Spherical linear interpolation (t ∈ [0, 1]) + * - ``equals(other, tol)`` + - Check equality within tolerance + * - ``is_identity(tol)`` + - Check if identity rotation + +Rotation Free Functions +======================= + +The following C++ free functions provide direct rotation operations without creating a ``Rotation`` +object. They are used internally and available in the C++ API. + +Vector and Matrix Rotation +-------------------------- + +.. doxygenfunction:: simcoon::rotate_vec(const arma::vec&, const arma::mat&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_vec(const arma::vec&, const double&, const int&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_mat(const arma::mat&, const arma::mat&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_mat(const arma::mat&, const double&, const int&) + :project: simcoon + +Stress and Strain Rotation +-------------------------- + +.. doxygenfunction:: simcoon::rotate_stress(const arma::vec&, const double&, const int&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_stress(const arma::vec&, const arma::mat&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_strain(const arma::vec&, const double&, const int&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_strain(const arma::vec&, const arma::mat&, const bool&) + :project: simcoon + +Stiffness and Compliance Rotation +--------------------------------- + +.. doxygenfunction:: simcoon::rotate_stiffness(const arma::mat&, const double&, const int&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_stiffness(const arma::mat&, const arma::mat&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_compliance(const arma::mat&, const double&, const int&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_compliance(const arma::mat&, const arma::mat&, const bool&) + :project: simcoon + +Strain and Stress Concentration Tensor Rotation +------------------------------------------------ + +.. doxygenfunction:: simcoon::rotate_strain_concentration(const arma::mat&, const double&, const int&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_strain_concentration(const arma::mat&, const arma::mat&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_stress_concentration(const arma::mat&, const double&, const int&, const bool&) + :project: simcoon + +.. doxygenfunction:: simcoon::rotate_stress_concentration(const arma::mat&, const arma::mat&, const bool&) + :project: simcoon + +Examples +======== + +Example 1: Material Orientation +------------------------------- + +Rotating material properties from local to global coordinates: + +.. code-block:: python + + import simcoon as smc + import numpy as np + + # Define orthotropic stiffness in material coordinates + L_local = smc.L_ortho([E1, E2, E3, nu12, nu13, nu23, G12, G13, G23], "EnuG") + + # Material orientation: Euler angles (radians) + psi, theta, phi = 0.5, 0.3, 0.7 + + # Using Rotation class + r = smc.Rotation.from_euler('ZXZ', [psi, theta, phi]) + L_global = r.apply_stiffness(L_local) + +Example 2: Stress Transformation +-------------------------------- + +Transforming stress from global to local coordinates: + +.. code-block:: python + + import simcoon as smc + import numpy as np + + # Stress in global coordinates (Voigt notation) + sigma_global = np.array([100.0, 50.0, 25.0, 10.0, 5.0, 2.0]) + + # Create rotation from orientation + r = smc.Rotation.from_euler('ZXZ', [0.5, 0.3, 0.7]) + + # Transform to local coordinates (inverse rotation) + sigma_local = r.inv().apply_stress(sigma_global) + +Example 3: Rotation Composition +------------------------------- + +Combining multiple rotations: + +.. code-block:: cpp + + #include + + using namespace simcoon; + + // First rotation: 45° around z-axis + Rotation r1 = Rotation::from_axis_angle(M_PI/4, 3); + + // Second rotation: 30° around x-axis + Rotation r2 = Rotation::from_axis_angle(M_PI/6, 1); + + // Combined rotation: apply r1 first, then r2 + Rotation r_combined = r2 * r1; + + // Apply to vector + arma::vec::fixed<3> v = {1, 0, 0}; + arma::vec::fixed<3> v_rot = r_combined.apply(v); + + // This is equivalent to: + arma::vec::fixed<3> v_rot2 = r2.apply(r1.apply(v)); + +Example 4: Interpolating Orientations +------------------------------------- + +Smooth interpolation between two orientations using SLERP: + +.. code-block:: python + + import simcoon as smc + import numpy as np + from scipy.spatial.transform import Slerp + + # Start and end orientations + r_start = smc.Rotation.identity() + r_end = smc.Rotation.from_euler('ZXZ', [np.pi/2, np.pi/4, 0]) + + # Interpolate at 10 steps + key_rots = smc.Rotation.concatenate([r_start, r_end]) + slerp = Slerp([0, 1], key_rots) + for t in np.linspace(0, 1, 10): + r_t = slerp(t) + R = r_t.as_matrix() + print(f"t={t:.1f}: angle={np.degrees(r_t.magnitude()):.1f}°") + +Theory +====== + +Quaternion Representation +------------------------- + +The ``Rotation`` class uses unit quaternions in scalar-last convention: + +.. math:: + + q = [q_x, q_y, q_z, q_w] = [\mathbf{v}, w] + +where :math:`\|\mathbf{q}\| = 1`. + +A rotation by angle :math:`\theta` around unit axis :math:`\mathbf{u}` is represented as: + +.. math:: + + q = \left[\mathbf{u} \sin\frac{\theta}{2}, \cos\frac{\theta}{2}\right] + +Quaternion to Rotation Matrix +----------------------------- + +The rotation matrix is computed from the quaternion as: + +.. math:: + + R = \begin{bmatrix} + 1 - 2(q_y^2 + q_z^2) & 2(q_x q_y - q_z q_w) & 2(q_x q_z + q_y q_w) \\ + 2(q_x q_y + q_z q_w) & 1 - 2(q_x^2 + q_z^2) & 2(q_y q_z - q_x q_w) \\ + 2(q_x q_z - q_y q_w) & 2(q_y q_z + q_x q_w) & 1 - 2(q_x^2 + q_y^2) + \end{bmatrix} + +Voigt Notation +-------------- + +For stress and strain tensors, simcoon uses Voigt notation: + +.. math:: + + \boldsymbol{\sigma} = [\sigma_{11}, \sigma_{22}, \sigma_{33}, \sigma_{12}, \sigma_{13}, \sigma_{23}]^T + +.. math:: + + \boldsymbol{\varepsilon} = [\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33}, 2\varepsilon_{12}, 2\varepsilon_{13}, 2\varepsilon_{23}]^T + +The 6×6 rotation matrices :math:`Q_\sigma` (for stress) and :math:`Q_\varepsilon` (for strain) +are related to the 3×3 rotation matrix :math:`R` and satisfy: + +.. math:: + + \boldsymbol{\sigma}' = Q_\sigma \boldsymbol{\sigma} + +.. math:: + + \boldsymbol{\varepsilon}' = Q_\varepsilon \boldsymbol{\varepsilon} + +.. math:: + + L' = Q_\sigma L Q_\sigma^T + +.. math:: + + M' = Q_\varepsilon M Q_\varepsilon^T + +SLERP Interpolation +------------------- + +Spherical Linear Interpolation (SLERP) provides smooth interpolation between rotations: + +.. math:: + + \text{slerp}(q_0, q_1, t) = \frac{\sin((1-t)\Omega)}{\sin\Omega} q_0 + \frac{\sin(t\Omega)}{\sin\Omega} q_1 + +where :math:`\Omega = \arccos(q_0 \cdot q_1)` is the angle between the quaternions. + +See Also +======== + +- :doc:`/continuum_mechanics/functions/doc_kinematics` - Finite strain kinematics +- :doc:`/continuum_mechanics/functions/doc_constitutive` - Constitutive matrices +- :doc:`/continuum_mechanics/functions/doc_transfer` - Voigt/tensor conversions +- :doc:`/continuum_mechanics/functions/doc_rotation` - Python user documentation diff --git a/docs/cpp_api/simulation_overview.md b/docs/cpp_api/simulation_overview.md index 2d8378512..033d4550d 100644 --- a/docs/cpp_api/simulation_overview.md +++ b/docs/cpp_api/simulation_overview.md @@ -195,12 +195,12 @@ Rotation operations for objective stress integration: - Active/passive rotation conventions **Key Functions:** -- `fillR_euler()` - Construct rotation matrix from Euler angles (ZXZ, ZYZ conventions) -- `fillR_angle()` - Rotation matrix from axis-angle representation +- `Rotation::from_euler()` - Construct rotation from Euler angles (ZXZ, ZYZ, XYZ, etc.) +- `Rotation::from_axis_angle()` - Rotation from axis-angle representation - `rotate_strain()` - Rotate strain tensor (Voigt notation) - `rotate_stress()` - Rotate stress tensor (Voigt notation) -- `rotateL()` - Rotate stiffness matrix -- `rotateM()` - Rotate compliance matrix +- `rotate_stiffness()` - Rotate stiffness matrix +- `rotate_compliance()` - Rotate compliance matrix #### **lagrange.hpp** Lagrange multiplier methods for constrained problems: diff --git a/examples/continuum_mechanics/rotation.py b/examples/continuum_mechanics/rotation.py index 44a4c16ba..4c0150d80 100644 --- a/examples/continuum_mechanics/rotation.py +++ b/examples/continuum_mechanics/rotation.py @@ -1,13 +1,16 @@ """ Rotation library Examples ~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The simcoon ``Rotation`` class inherits from ``scipy.spatial.transform.Rotation``, +so all scipy factory methods and operations are available directly. Plain scipy +``Rotation`` objects can be upgraded with ``from_scipy()`` or passed directly to +C++ methods that accept rotations. """ import numpy as np -import matplotlib.pyplot as plt -from scipy.spatial.transform import Rotation as R +from scipy.spatial.transform import Rotation as ScipyRotation import simcoon as sim -import os # %% @@ -19,65 +22,58 @@ m = np.eye(3) angle = np.pi / 4 # 45 degrees axis = 2 # y-axis (1=x, 2=y, 3=z) -copy = True active = True # %% -# 1. Rotate a vector using a rotation matrix +# 1. Rotate a vector using the Rotation class # -# This example shows how to build a rotation matrix and rotate a vector. -# It uses :func:`simcoon.fillR_angle` and :func:`simcoon.rotate_vec_R`. - +# This example shows how to build a rotation from an angle and axis +# and rotate a vector using :meth:`simcoon.Rotation.apply`. -Rmat = sim.fillR_angle(angle, axis, active, copy) -v_rot1 = sim.rotate_vec_R(v, Rmat, copy) -print("Rotated vector (using R):", v_rot1) -print("Rotation matrix (using R):", Rmat) +rot = sim.Rotation.from_axis_angle(angle, axis) +Rmat = rot.as_matrix() +v_rot1 = rot.apply(v) +print("Rotated vector:", v_rot1) +print("Rotation matrix:\n", Rmat) # %% -# 2. Rotate a vector using angle/axis +# 2. Rotate a vector using Euler angles (scipy syntax) # -# This example shows how to rotate a vector using an angle and an axis directly. -# It uses :func:`simcoon.rotate_vec_angle`. - +# Uppercase letters indicate intrinsic rotations (rotating frame). +# Lowercase letters indicate extrinsic rotations (fixed frame). -v_rot2 = sim.rotate_vec_angle(v, angle, axis, copy) -print("Rotated vector (using angle/axis):", v_rot2) +rot2 = sim.Rotation.from_euler("ZXZ", [np.pi / 6, np.pi / 4, np.pi / 3]) +v_rot2 = rot2.apply(v) +print("Rotated vector:", v_rot2) # %% -# 3. Rotate a matrix using a rotation matrix +# 3. Rotate a matrix (3x3 tensor) # -# This example shows how to rotate a matrix using a rotation matrix. -# It uses :func:`simcoon.rotate_mat_R`. +# This example shows how to rotate a 3x3 tensor using +# :meth:`simcoon.Rotation.apply_tensor`. - -m_rot1 = sim.rotate_mat_R(m, Rmat, copy) -print("Rotated matrix (using R):\n", m_rot1) +m_rot1 = rot.apply_tensor(m) +print("Rotated matrix:\n", m_rot1) # %% -# 4. Rotate a matrix using angle/axis -# -# This example shows how to rotate a matrix using an angle and an axis directly. -# It uses :func:`simcoon.rotate_mat_angle`. - +# 4. Rotate a matrix using a rotation from Euler angles -m_rot2 = sim.rotate_mat_angle(m, angle, axis, copy) -print("Rotated matrix (using angle/axis):\n", m_rot2) +rot_euler = sim.Rotation.from_euler("ZXZ", [np.pi / 6, np.pi / 4, np.pi / 3]) +m_rot2 = rot_euler.apply_tensor(m) +print("Rotated matrix:\n", m_rot2) # %% # 5. Build a rotation matrix from Euler angles -# -# This example shows how to create a rotation matrix from Euler angles. -# It uses :func:`simcoon.fillR_euler`. psi, theta, phi = np.pi / 6, np.pi / 4, np.pi / 3 -R_euler = sim.fillR_euler(psi, theta, phi, active, "zxz", copy) -print("Rotation matrix from Euler angles (zxz):\n", R_euler) +r_euler = sim.Rotation.from_euler("ZXZ", [psi, theta, phi]) +R_euler = r_euler.as_matrix() +print("Rotation matrix from Euler angles (ZXZ):\n", R_euler) # %% @@ -88,81 +84,63 @@ # %% -# 6. Rotate a stress vector (single and batch) +# 6. Rotate a stress vector # -# This example uses :func:`simcoon.rotate_stress_angle`. +# This example uses :meth:`simcoon.Rotation.apply_stress`. stress = np.array([1, 2, 3, 4, 5, 6], dtype=float) -stress_batch = np.stack([stress, stress * 2], axis=1) -rot_stress1 = sim.rotate_stress_angle(stress, angle, axis, active, copy) -rot_stress2 = sim.rotate_stress_angle(stress_batch, angle, axis, active, copy) -print("Rotated stress (single):", rot_stress1) -print("Rotated stress (batch):\n", rot_stress2) +rot_stress1 = rot.apply_stress(stress, active) +print("Rotated stress:", rot_stress1) # %% -# 7. Rotate a strain vector (single and batch) +# 7. Rotate a strain vector # -# This example uses :func:`simcoon.rotate_strain_angle`. +# This example uses :meth:`simcoon.Rotation.apply_strain`. strain = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6], dtype=float) -strain_batch = np.stack([strain, strain * 2], axis=1) -rot_strain1 = sim.rotate_strain_angle(strain, angle, axis, active, copy) -rot_strain2 = sim.rotate_strain_angle(strain_batch, angle, axis, active, copy) -print("Rotated strain (single):", rot_strain1) -print("Rotated strain (batch):\n", rot_strain2) +rot_strain1 = rot.apply_strain(strain, active) +print("Rotated strain:", rot_strain1) # %% # 8. Rotate a stiffness matrix (L) # -# This example uses both :func:`simcoon.rotateL_angle` and -# :func:`simcoon.rotateL_R`. +# This example uses :meth:`simcoon.Rotation.apply_stiffness`. L6 = np.eye(6) -rotL1 = sim.rotateL_angle(L6, angle, axis, active, copy) -rotL2 = sim.rotateL_R(L6, Rmat, active, copy) -print("Rotated L (angle):\n", rotL1) -print("Rotated L (R):\n", rotL2) +rotL1 = rot.apply_stiffness(L6, active) +print("Rotated L:\n", rotL1) # %% # 9. Rotate a compliance matrix (M) # -# This example uses both :func:`simcoon.rotateM_angle` and -# :func:`simcoon.rotateM_R`. +# This example uses :meth:`simcoon.Rotation.apply_compliance`. M6 = np.eye(6) -rotM1 = sim.rotateM_angle(M6, angle, axis, active, copy) -rotM2 = sim.rotateM_R(M6, Rmat, active, copy) -print("Rotated M (angle):\n", rotM1) -print("Rotated M (R):\n", rotM2) +rotM1 = rot.apply_compliance(M6, active) +print("Rotated M:\n", rotM1) # %% -# 10. Rotate a strain concentration matrix (A) +# 10. Rotate a strain concentration tensor (A) # -# This example uses both :func:`simcoon.rotateA_angle` and -# :func:`simcoon.rotateA_R`. +# This example uses :meth:`simcoon.Rotation.apply_strain_concentration`. A6 = np.eye(6) -rotA1 = sim.rotateA_angle(A6, angle, axis, active, copy) -rotA2 = sim.rotateA_R(A6, Rmat, active, copy) -print("Rotated A (angle):\n", rotA1) -print("Rotated A (R):\n", rotA2) +rotA1 = rot.apply_strain_concentration(A6, active) +print("Rotated A:\n", rotA1) # %% -# 11. Rotate a stress concentration matrix (B) +# 11. Rotate a stress concentration tensor (B) # -# This example uses both :func:`simcoon.rotateB_angle` and -# :func:`simcoon.rotateB_R`. +# This example uses :meth:`simcoon.Rotation.apply_stress_concentration`. B6 = np.eye(6) -rotB1 = sim.rotateB_angle(B6, angle, axis, active, copy) -rotB2 = sim.rotateB_R(B6, Rmat, active, copy) -print("Rotated B (angle):\n", rotB1) -print("Rotated B (R):\n", rotB2) +rotB1 = rot.apply_stress_concentration(B6, active) +print("Rotated B:\n", rotB1) # %% @@ -184,14 +162,112 @@ print(x) alpha = np.pi / 4.0 +rot_z = sim.Rotation.from_axis_angle(alpha, 3) +L_rotate = rot_z.apply_stiffness(L) +print(np.array_str(L_rotate, suppress_small=True)) + -rot_matrix = np.array( - [[np.cos(alpha), -np.sin(alpha), 0.0], [np.sin(alpha), np.cos(alpha), 0], [0, 0, 1]] -) -# print(R) +# %% +# Scipy interoperability +# +# ``simcoon.Rotation`` inherits from ``scipy.spatial.transform.Rotation``, +# so every scipy feature is available. Plain scipy ``Rotation`` objects can +# be upgraded to ``simcoon.Rotation`` via ``from_scipy()`` to unlock the +# mechanics methods. -L_rotate = sim.rotateL_R(L, rot_matrix) -L_rotate_angle = sim.rotateL_angle(L, alpha, axis=3) -print(np.array_str(L_rotate, suppress_small=True)) -print(np.array_str(L_rotate_angle, suppress_small=True)) +# %% +# 13. simcoon.Rotation IS a scipy Rotation +# +# Type checks and scipy methods work directly. + +r = sim.Rotation.from_axis_angle(np.pi / 4, 3) +print("isinstance(r, ScipyRotation):", isinstance(r, ScipyRotation)) +print("Quaternion:", r.as_quat()) +print("Rotation vector:", r.as_rotvec()) +print("Euler ZXZ:", r.as_euler("ZXZ", degrees=True)) + + +# %% +# 14. Upgrade a plain scipy Rotation to simcoon +# +# Use ``from_scipy()`` to add mechanics capabilities to an existing scipy +# rotation (e.g. one returned by a third-party library). + +scipy_rot = ScipyRotation.from_euler("z", 45, degrees=True) +print("Type before:", type(scipy_rot)) + +smc_rot = sim.Rotation.from_scipy(scipy_rot) +print("Type after:", type(smc_rot)) + +# The quaternion is preserved exactly +print("Quaternions match:", np.allclose(scipy_rot.as_quat(), smc_rot.as_quat())) + +# Now mechanics methods are available +L_iso = sim.L_iso([70000, 0.3], "Enu") +L_rotated = smc_rot.apply_stiffness(L_iso) +print("Rotated stiffness L[0,0]:", L_rotated[0, 0]) + + +# %% +# 15. Verify scipy and simcoon produce the same rotation +# +# A simcoon rotation and a scipy rotation built from the same parameters +# give identical stiffness results. + +angle_test = np.pi / 6 +r_sim = sim.Rotation.from_axis_angle(angle_test, 3) +r_scipy = sim.Rotation.from_scipy(ScipyRotation.from_rotvec([0, 0, angle_test])) + +L_sim = r_sim.apply_stiffness(L_iso) +L_scipy = r_scipy.apply_stiffness(L_iso) +print("Stiffness matrices match:", np.allclose(L_sim, L_scipy)) + + +# %% +# 16. Compose simcoon and scipy rotations +# +# Because simcoon.Rotation inherits ``__mul__`` from scipy, composition +# with either type works directly. + +r1 = sim.Rotation.from_axis_angle(np.pi / 4, 3) # simcoon +r2 = ScipyRotation.from_euler("x", 30, degrees=True) # plain scipy + +r_composed = r1 * r2 # returns simcoon.Rotation +print("Composed type:", type(r_composed)) +print("Composed is simcoon Rotation:", isinstance(r_composed, sim.Rotation)) + + +# %% +# 17. Batch rotations and mean (scipy features) +# +# Generate random orientations and compute their mean — all scipy batch +# operations work out of the box. + +rots = sim.Rotation.random(100) +r_mean = rots.mean() +print("Mean rotation type:", type(r_mean)) +print("Mean rotation angle:", np.degrees(r_mean.magnitude()), "degrees") + + +# %% +# 18. Slerp interpolation +# +# Use scipy's ``Slerp`` class for multi-keyframe interpolation, or the +# convenience ``slerp_to()`` method for simple two-rotation interpolation. + +from scipy.spatial.transform import Slerp + +r_start = sim.Rotation.identity() +r_end = sim.Rotation.from_axis_angle(np.pi / 2, 3) + +# Via scipy Slerp class +key_rots = sim.Rotation.concatenate([r_start, r_end]) +slerp = Slerp([0, 1], key_rots) +for t in [0.0, 0.25, 0.5, 0.75, 1.0]: + r_t = slerp(t) + print(f" t={t:.2f}: angle = {np.degrees(r_t.magnitude()):.1f} deg") + +# Via convenience method +r_mid = r_start.slerp_to(r_end, 0.5) +print("Midpoint angle:", np.degrees(r_mid.magnitude()), "degrees") diff --git a/include/simcoon/Simulation/Maths/rotation.hpp b/include/simcoon/Simulation/Maths/rotation.hpp index ec3983ec0..022745998 100755 --- a/include/simcoon/Simulation/Maths/rotation.hpp +++ b/include/simcoon/Simulation/Maths/rotation.hpp @@ -1,50 +1,466 @@ /* This file is part of simcoon. - + simcoon is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. - + simcoon is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - + You should have received a copy of the GNU General Public License along with simcoon. If not, see . - + */ ///@file rotation.hpp -///@brief rotation of a Voigt tensor -///@version 1.0 +///@brief Rotation class and Voigt tensor rotation utilities. +///@version 2.0 #pragma once #include +#include namespace simcoon{ /** * @file rotation.hpp - * @brief Mathematical utility functions. + * @brief Rotation class and convenience free functions for rotating + * vectors, tensors, and Voigt-notation quantities. */ /** @addtogroup maths * @{ */ +// ========================================================================= +// Rotation Class +// ========================================================================= + +/** + * @class Rotation + * @brief A class representing 3D rotations using unit quaternions. + * + * Inspired by scipy.spatial.transform.Rotation, this class provides a unified + * interface for working with rotations in various representations (quaternion, + * rotation matrix, Euler angles, rotation vector) while maintaining high + * numerical precision and performance. + * + * Internal representation uses unit quaternions in scalar-last convention [qx, qy, qz, qw]. + * + * @details Example usage: + * @code + * // Create rotation from Euler angles + * Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + * + * // Apply to a vector + * arma::vec::fixed<3> v = {1, 0, 0}; + * arma::vec::fixed<3> v_rot = r.apply(v); + * + * // Compose rotations + * Rotation r2 = Rotation::from_axis_angle(M_PI/4, 3); + * Rotation r_combined = r * r2; + * + * // Apply to Voigt stress tensor + * arma::vec::fixed<6> sigma = {100, 50, 0, 25, 0, 0}; + * arma::vec::fixed<6> sigma_rot = r.apply_stress(sigma); + * @endcode + */ +class Rotation { +private: + arma::vec::fixed<4> _quat; // [qx, qy, qz, qw] - scalar-last convention, stack allocated + + // Private constructor from quaternion (assumes already normalized) + explicit Rotation(const arma::vec::fixed<4>& quat); + +public: + /** + * @brief Default constructor creating identity rotation. + */ + Rotation(); + + /** + * @brief Copy constructor. + */ + Rotation(const Rotation& other) = default; + + /** + * @brief Move constructor. + */ + Rotation(Rotation&& other) noexcept = default; + + /** + * @brief Copy assignment operator. + */ + Rotation& operator=(const Rotation& other) = default; + + /** + * @brief Move assignment operator. + */ + Rotation& operator=(Rotation&& other) noexcept = default; + + /** + * @brief Destructor. + */ + ~Rotation() = default; + + // ========================================================================= + // Factory Methods + // ========================================================================= + + /** + * @brief Create an identity rotation. + * @return Identity rotation (no rotation) + */ + static Rotation identity(); + + /** + * @brief Create rotation from a unit quaternion. + * @param quat Quaternion in [qx, qy, qz, qw] format (scalar-last) + * @return Rotation object + * @note The quaternion will be normalized if not already unit length + */ + static Rotation from_quat(const arma::vec::fixed<4>& quat); + + /** + * @brief Create rotation from a unit quaternion. + * @param quat Quaternion in [qx, qy, qz, qw] format (scalar-last), dynamic size + * @return Rotation object + */ + static Rotation from_quat(const arma::vec& quat); + + /** + * @brief Create rotation from a 3x3 rotation matrix. + * @param R 3x3 rotation matrix (must be orthogonal with det=1) + * @return Rotation object + */ + static Rotation from_matrix(const arma::mat::fixed<3,3>& R); + + /** + * @brief Create rotation from a 3x3 rotation matrix. + * @param R 3x3 rotation matrix (dynamic size) + * @return Rotation object + */ + static Rotation from_matrix(const arma::mat& R); + + /** + * @brief Create rotation from Euler angles. + * @param psi First Euler angle in radians (or degrees if degrees=true) + * @param theta Second Euler angle in radians (or degrees if degrees=true) + * @param phi Third Euler angle in radians (or degrees if degrees=true) + * @param conv Euler angle convention: "zxz", "zyz", "xyz", "xzy", "yxz", "yzx", "zxy", "zyx", + * "xyx", "xzx", "yxy", "yzy", or "user" + * @param intrinsic If true (default), intrinsic rotations (rotating frame); + * if false, extrinsic rotations (fixed frame) + * @param degrees If true, angles are in degrees; if false (default), radians + * @return Rotation object + */ + static Rotation from_euler(double psi, double theta, double phi, + const std::string& conv, + bool intrinsic = true, bool degrees = false); + + /** + * @brief Create rotation from Euler angles (vector input). + * @param angles 3-element vector of Euler angles + * @param conv Euler angle convention + * @param intrinsic If true, intrinsic rotations; if false, extrinsic rotations + * @param degrees If true, angles are in degrees + * @return Rotation object + */ + static Rotation from_euler(const arma::vec::fixed<3>& angles, + const std::string& conv, + bool intrinsic = true, bool degrees = false); + + /** + * @brief Create rotation from a rotation vector (axis-angle representation). + * @param rotvec Rotation vector where direction is axis and magnitude is angle + * @param degrees If true, magnitude is in degrees; if false (default), radians + * @return Rotation object + */ + static Rotation from_rotvec(const arma::vec::fixed<3>& rotvec, bool degrees = false); + + /** + * @brief Create rotation from a rotation vector (dynamic size). + * @param rotvec Rotation vector (3 elements) + * @param degrees If true, magnitude is in degrees + * @return Rotation object + */ + static Rotation from_rotvec(const arma::vec& rotvec, bool degrees = false); + + /** + * @brief Create rotation around a single axis. + * @param angle Rotation angle in radians (or degrees if degrees=true) + * @param axis Axis of rotation: 1=x, 2=y, 3=z + * @param degrees If true, angle is in degrees + * @return Rotation object + */ + static Rotation from_axis_angle(double angle, int axis, bool degrees = false); + + /** + * @brief Create a random rotation (uniformly distributed). + * @return Random rotation object + */ + static Rotation random(); + + // ========================================================================= + // Conversion Methods + // ========================================================================= + + /** + * @brief Get rotation as a unit quaternion. + * @return Quaternion in [qx, qy, qz, qw] format (scalar-last) + */ + inline arma::vec::fixed<4> as_quat() const { return _quat; } + + /** + * @brief Get rotation as a 3x3 rotation matrix. + * @return 3x3 orthogonal rotation matrix + */ + arma::mat::fixed<3,3> as_matrix() const; + + /** + * @brief Get rotation as Euler angles. + * @param conv Euler angle convention + * @param intrinsic If true, return intrinsic angles; if false, extrinsic + * @param degrees If true, return angles in degrees + * @return 3-element vector of Euler angles [psi, theta, phi] + */ + arma::vec::fixed<3> as_euler(const std::string& conv, + bool intrinsic = true, bool degrees = false) const; + + /** + * @brief Get rotation as a rotation vector. + * @param degrees If true, magnitude is in degrees + * @return Rotation vector (axis * angle) + */ + arma::vec::fixed<3> as_rotvec(bool degrees = false) const; + + /** + * @brief Get 6x6 rotation matrix for stress tensors in Voigt notation. + * @param active If true (default), active rotation; if false, passive + * @return 6x6 stress rotation matrix + */ + arma::mat::fixed<6,6> as_voigt_stress_rotation(bool active = true) const; + + /** + * @brief Get 6x6 rotation matrix for strain tensors in Voigt notation. + * @param active If true (default), active rotation; if false, passive + * @return 6x6 strain rotation matrix + */ + arma::mat::fixed<6,6> as_voigt_strain_rotation(bool active = true) const; + + // ========================================================================= + // Apply Methods (3D objects) + // ========================================================================= + + /** + * @brief Apply rotation to a 3D vector. + * @param v 3D vector to rotate + * @param inverse If true, apply inverse rotation + * @return Rotated vector + */ + arma::vec::fixed<3> apply(const arma::vec::fixed<3>& v, bool inverse = false) const; + + /** + * @brief Apply rotation to a 3D vector (dynamic size). + * @param v 3D vector to rotate + * @param inverse If true, apply inverse rotation + * @return Rotated vector + */ + arma::vec apply(const arma::vec& v, bool inverse = false) const; + + /** + * @brief Apply rotation to a 3x3 tensor (matrix). + * @param m 3x3 tensor to rotate + * @param inverse If true, apply inverse rotation + * @return Rotated tensor: R * m * R^T (or R^T * m * R for inverse) + */ + arma::mat::fixed<3,3> apply_tensor(const arma::mat::fixed<3,3>& m, bool inverse = false) const; + + /** + * @brief Apply rotation to a 3x3 tensor (dynamic size). + * @param m 3x3 tensor to rotate + * @param inverse If true, apply inverse rotation + * @return Rotated tensor + */ + arma::mat apply_tensor(const arma::mat& m, bool inverse = false) const; + + // ========================================================================= + // Apply Methods (Voigt notation) + // ========================================================================= + + /** + * @brief Apply rotation to a stress vector in Voigt notation. + * @param sigma 6-component stress vector [s11, s22, s33, s12, s13, s23] + * @param active If true (default), active rotation + * @return Rotated stress vector + */ + arma::vec::fixed<6> apply_stress(const arma::vec::fixed<6>& sigma, bool active = true) const; + + /** + * @brief Apply rotation to a stress vector (dynamic size). + * @param sigma 6-component stress vector + * @param active If true, active rotation + * @return Rotated stress vector + */ + arma::vec apply_stress(const arma::vec& sigma, bool active = true) const; + + /** + * @brief Apply rotation to a strain vector in Voigt notation. + * @param epsilon 6-component strain vector [e11, e22, e33, 2*e12, 2*e13, 2*e23] + * @param active If true (default), active rotation + * @return Rotated strain vector + */ + arma::vec::fixed<6> apply_strain(const arma::vec::fixed<6>& epsilon, bool active = true) const; + + /** + * @brief Apply rotation to a strain vector (dynamic size). + * @param epsilon 6-component strain vector + * @param active If true, active rotation + * @return Rotated strain vector + */ + arma::vec apply_strain(const arma::vec& epsilon, bool active = true) const; + + /** + * @brief Apply rotation to a 6x6 stiffness matrix. + * @param L 6x6 stiffness matrix in Voigt notation + * @param active If true (default), active rotation + * @return Rotated stiffness matrix: voigt_stress * L * voigt_stress^T + */ + arma::mat::fixed<6,6> apply_stiffness(const arma::mat::fixed<6,6>& L, bool active = true) const; + + /** + * @brief Apply rotation to a stiffness matrix (dynamic size). + * @param L 6x6 stiffness matrix + * @param active If true, active rotation + * @return Rotated stiffness matrix + */ + arma::mat apply_stiffness(const arma::mat& L, bool active = true) const; + + /** + * @brief Apply rotation to a 6x6 compliance matrix. + * @param M 6x6 compliance matrix in Voigt notation + * @param active If true (default), active rotation + * @return Rotated compliance matrix: voigt_strain * M * voigt_strain^T + */ + arma::mat::fixed<6,6> apply_compliance(const arma::mat::fixed<6,6>& M, bool active = true) const; + + /** + * @brief Apply rotation to a compliance matrix (dynamic size). + * @param M 6x6 compliance matrix + * @param active If true, active rotation + * @return Rotated compliance matrix + */ + arma::mat apply_compliance(const arma::mat& M, bool active = true) const; + + /** + * @brief Apply rotation to a 6x6 strain concentration tensor. + * @param A 6x6 strain concentration tensor in Voigt notation + * @param active If true (default), active rotation + * @return Rotated strain concentration tensor: voigt_strain * A * voigt_stress^T + */ + arma::mat::fixed<6,6> apply_strain_concentration(const arma::mat::fixed<6,6>& A, bool active = true) const; + + /** + * @brief Apply rotation to a strain concentration tensor (dynamic size). + * @param A 6x6 strain concentration tensor + * @param active If true, active rotation + * @return Rotated strain concentration tensor + */ + arma::mat apply_strain_concentration(const arma::mat& A, bool active = true) const; + + /** + * @brief Apply rotation to a 6x6 stress concentration tensor. + * @param B 6x6 stress concentration tensor in Voigt notation + * @param active If true (default), active rotation + * @return Rotated stress concentration tensor: voigt_stress * B * voigt_strain^T + */ + arma::mat::fixed<6,6> apply_stress_concentration(const arma::mat::fixed<6,6>& B, bool active = true) const; + + /** + * @brief Apply rotation to a stress concentration tensor (dynamic size). + * @param B 6x6 stress concentration tensor + * @param active If true, active rotation + * @return Rotated stress concentration tensor + */ + arma::mat apply_stress_concentration(const arma::mat& B, bool active = true) const; + + // ========================================================================= + // Operations + // ========================================================================= + + /** + * @brief Get the inverse (conjugate) rotation. + * @return Inverse rotation + */ + Rotation inv() const; + + /** + * @brief Get the magnitude (rotation angle) of this rotation. + * @param degrees If true, return in degrees + * @return Rotation angle in radians (or degrees) + */ + double magnitude(bool degrees = false) const; + + /** + * @brief Compose this rotation with another (this * other). + * @param other Rotation to compose with + * @return Composed rotation + * @details The composition represents applying 'other' first, then 'this' + */ + Rotation operator*(const Rotation& other) const; + + /** + * @brief Compose this rotation with another in-place. + * @param other Rotation to compose with + * @return Reference to this rotation + */ + Rotation& operator*=(const Rotation& other); + + /** + * @brief Spherical linear interpolation between this rotation and another. + * @param other Target rotation + * @param t Interpolation parameter in [0, 1] (0 = this, 1 = other) + * @return Interpolated rotation + */ + Rotation slerp(const Rotation& other, double t) const; + + /** + * @brief Check if this rotation equals another within tolerance. + * @param other Rotation to compare with + * @param tol Tolerance for comparison + * @return True if rotations are equal (or antipodal quaternions) + */ + inline bool equals(const Rotation& other, double tol = 1e-12) const { + // Quaternions q and -q represent the same rotation + double diff1 = arma::norm(_quat - other._quat); + double diff2 = arma::norm(_quat + other._quat); + return (diff1 < tol) || (diff2 < tol); + } + + /** + * @brief Check if this rotation is identity (no rotation). + * @param tol Tolerance for comparison + * @return True if this is the identity rotation + */ + inline bool is_identity(double tol = 1e-12) const { + // Identity quaternion is [0, 0, 0, 1] or [0, 0, 0, -1] + return std::abs(std::abs(_quat(3)) - 1.0) < tol; + } +}; + +// ========================================================================= +// Convenience Free Functions +// ========================================================================= + /** * @brief Rotates a 3D vector using a rotation matrix. * @param v The input 3D vector (arma::vec) * @param R The 3x3 rotation matrix (arma::mat) * @return The rotated vector (arma::vec) - * @details Example: - * @code - * vec v = {1., 0., 0.}; - * mat R = fillR(M_PI/4., 3); - * vec v_rot = rotate_vec(v, R); - * @endcode */ arma::vec rotate_vec(const arma::vec &v, const arma::mat &R); @@ -54,11 +470,6 @@ arma::vec rotate_vec(const arma::vec &v, const arma::mat &R); * @param angle The rotation angle in radians (double) * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) * @return The rotated vector (arma::vec) - * @details Example: - * @code - * vec v = {1., 0., 0.}; - * vec v_rot = rotate_vec(v, M_PI/4., 3); // Rotate 45° around z-axis - * @endcode */ arma::vec rotate_vec(const arma::vec &v, const double &angle, const int &axis); @@ -67,12 +478,6 @@ arma::vec rotate_vec(const arma::vec &v, const double &angle, const int &axis); * @param m The input 3x3 matrix (arma::mat) * @param R The 3x3 rotation matrix (arma::mat) * @return The rotated matrix (arma::mat) - * @details The rotation is performed as \f$ \mathbf{m}' = \mathbf{R} \cdot \mathbf{m} \cdot \mathbf{R}^T \f$ - * @code - * mat m = eye(3,3); - * mat R = fillR(M_PI/4., 3); - * mat m_rot = rotate_mat(m, R); - * @endcode */ arma::mat rotate_mat(const arma::mat &m, const arma::mat &R); @@ -86,81 +491,52 @@ arma::mat rotate_mat(const arma::mat &m, const arma::mat &R); arma::mat rotate_mat(const arma::mat &m, const double &angle, const int &axis); /** - * @brief Generates a 3x3 rotation matrix for rotation around a single axis. - * @param angle The rotation angle in radians (double) - * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) - * @param active If true (default), active (alibi) rotation; if false, passive (alias) rotation (bool) - * @return The 3x3 rotation matrix (arma::mat) - * @details Example: - * @code - * mat R = fillR(M_PI/4., 3); // 45° rotation around z-axis - * @endcode - */ -arma::mat fillR(const double &angle, const int &axis, const bool &active = true); - -/** - * @brief Generates a 3x3 rotation matrix using Euler angles. - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @param active If true (default), active (alibi) rotation; if false, passive (alias) rotation (bool) - * @param conv The Euler angle convention: "zxz" (default), "zyz", "xyz", etc. Use "user" for custom convention (std::string) - * @return The 3x3 rotation matrix (arma::mat) - * @details Example: - * @code - * mat R = fillR(M_PI/4., M_PI/6., M_PI/3., true, "zxz"); - * @endcode - */ -arma::mat fillR(const double &psi, const double &theta, const double &phi, const bool &active = true, const std::string &conv = "zxz"); - -/** - * @brief Generates a 6x6 rotation matrix for stress tensors in Voigt notation. + * @brief Rotates a stress vector in Voigt notation. + * @param sigma The 6-component stress vector (arma::vec) * @param angle The rotation angle in radians (double) * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) * @param active If true (default), active rotation (bool) - * @return The 6x6 rotation matrix for stress (arma::mat) + * @return The rotated stress vector (arma::vec) */ -arma::mat fillQS(const double &angle, const int &axis, const bool &active = true); +arma::vec rotate_stress(const arma::vec &sigma, const double &angle, const int &axis, const bool &active = true); /** - * @brief Generates a 6x6 rotation matrix for stress tensors from a 3x3 rotation matrix. + * @brief Rotates a stress vector using a 3x3 rotation matrix. + * @param sigma The 6-component stress vector (arma::vec) * @param R The 3x3 rotation matrix (arma::mat) * @param active If true (default), active rotation (bool) - * @return The 6x6 rotation matrix for stress (arma::mat) + * @return The rotated stress vector (arma::vec) */ -arma::mat fillQS(const arma::mat &R, const bool &active = true); - +arma::vec rotate_stress(const arma::vec &sigma, const arma::mat &R, const bool &active = true); + /** - * @brief Generates a 6x6 rotation matrix for strain tensors in Voigt notation. + * @brief Rotates a strain vector in Voigt notation. + * @param epsilon The 6-component strain vector (arma::vec) * @param angle The rotation angle in radians (double) * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) * @param active If true (default), active rotation (bool) - * @return The 6x6 rotation matrix for strain (arma::mat) + * @return The rotated strain vector (arma::vec) */ -arma::mat fillQE(const double &angle, const int &axis, const bool &active = true); +arma::vec rotate_strain(const arma::vec &epsilon, const double &angle, const int &axis, const bool &active = true); /** - * @brief Generates a 6x6 rotation matrix for strain tensors from a 3x3 rotation matrix. + * @brief Rotates a strain vector using a 3x3 rotation matrix. + * @param epsilon The 6-component strain vector (arma::vec) * @param R The 3x3 rotation matrix (arma::mat) * @param active If true (default), active rotation (bool) - * @return The 6x6 rotation matrix for strain (arma::mat) + * @return The rotated strain vector (arma::vec) */ -arma::mat fillQE(const arma::mat &R, const bool &active = true); +arma::vec rotate_strain(const arma::vec &epsilon, const arma::mat &R, const bool &active = true); /** - * @brief Rotates a 6x6 stiffness matrix. + * @brief Rotates a 6x6 stiffness matrix by an angle around a given axis. * @param L The 6x6 stiffness matrix in Voigt notation (arma::mat) * @param angle The rotation angle in radians (double) * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) * @param active If true (default), active rotation (bool) * @return The rotated 6x6 stiffness matrix (arma::mat) - * @details Example: - * @code - * mat L = L_iso(E, nu, "EnuG"); - * mat L_rot = rotateL(L, M_PI/4., 3); - * @endcode */ -arma::mat rotateL(const arma::mat &L, const double &angle, const int &axis, const bool &active = true); +arma::mat rotate_stiffness(const arma::mat &L, const double &angle, const int &axis, const bool &active = true); /** * @brief Rotates a 6x6 stiffness matrix using a 3x3 rotation matrix. @@ -169,17 +545,17 @@ arma::mat rotateL(const arma::mat &L, const double &angle, const int &axis, cons * @param active If true (default), active rotation (bool) * @return The rotated 6x6 stiffness matrix (arma::mat) */ -arma::mat rotateL(const arma::mat &L, const arma::mat &R, const bool &active = true); +arma::mat rotate_stiffness(const arma::mat &L, const arma::mat &R, const bool &active = true); /** - * @brief Rotates a 6x6 compliance matrix. + * @brief Rotates a 6x6 compliance matrix by an angle around a given axis. * @param M The 6x6 compliance matrix in Voigt notation (arma::mat) * @param angle The rotation angle in radians (double) * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) * @param active If true (default), active rotation (bool) * @return The rotated 6x6 compliance matrix (arma::mat) */ -arma::mat rotateM(const arma::mat &M, const double &angle, const int &axis, const bool &active = true); +arma::mat rotate_compliance(const arma::mat &M, const double &angle, const int &axis, const bool &active = true); /** * @brief Rotates a 6x6 compliance matrix using a 3x3 rotation matrix. @@ -188,210 +564,45 @@ arma::mat rotateM(const arma::mat &M, const double &angle, const int &axis, cons * @param active If true (default), active rotation (bool) * @return The rotated 6x6 compliance matrix (arma::mat) */ -arma::mat rotateM(const arma::mat &M, const arma::mat &R, const bool &active = true); - -/** - * @brief Rotates a 6x6 strain localization matrix A. - * @param A The 6x6 strain localization matrix (arma::mat) - * @param angle The rotation angle in radians (double) - * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) - * @param active If true (default), active rotation (bool) - * @return The rotated 6x6 strain localization matrix (arma::mat) - */ -arma::mat rotateA(const arma::mat &A, const double &angle, const int &axis, const bool &active = true); - -/** - * @brief Rotates a 6x6 strain localization matrix A using a 3x3 rotation matrix. - * @param A The 6x6 strain localization matrix (arma::mat) - * @param R The 3x3 rotation matrix (arma::mat) - * @param active If true (default), active rotation (bool) - * @return The rotated 6x6 strain localization matrix (arma::mat) - */ -arma::mat rotateA(const arma::mat &A, const arma::mat &R, const bool &active = true); +arma::mat rotate_compliance(const arma::mat &M, const arma::mat &R, const bool &active = true); /** - * @brief Rotates a 6x6 stress localization matrix B. - * @param B The 6x6 stress localization matrix (arma::mat) + * @brief Rotates a 6x6 strain concentration matrix by an angle around a given axis. + * @param A The 6x6 strain concentration matrix (arma::mat) * @param angle The rotation angle in radians (double) * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) * @param active If true (default), active rotation (bool) - * @return The rotated 6x6 stress localization matrix (arma::mat) + * @return The rotated 6x6 strain concentration matrix (arma::mat) */ -arma::mat rotateB(const arma::mat &B, const double &angle, const int &axis, const bool &active = true); +arma::mat rotate_strain_concentration(const arma::mat &A, const double &angle, const int &axis, const bool &active = true); /** - * @brief Rotates a 6x6 stress localization matrix B using a 3x3 rotation matrix. - * @param B The 6x6 stress localization matrix (arma::mat) + * @brief Rotates a 6x6 strain concentration matrix using a 3x3 rotation matrix. + * @param A The 6x6 strain concentration matrix (arma::mat) * @param R The 3x3 rotation matrix (arma::mat) * @param active If true (default), active rotation (bool) - * @return The rotated 6x6 stress localization matrix (arma::mat) - */ -arma::mat rotateB(const arma::mat &B, const arma::mat &R, const bool &active = true); - -/** - * @brief Rotates a stress vector in Voigt notation. - * @param sigma The 6-component stress vector (arma::vec) - * @param angle The rotation angle in radians (double) - * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) - * @param active If true (default), active rotation (bool) - * @return The rotated stress vector (arma::vec) - * @details Example: - * @code - * vec sigma = {100., 50., 0., 25., 0., 0.}; - * vec sigma_rot = rotate_stress(sigma, M_PI/4., 3); - * @endcode + * @return The rotated 6x6 strain concentration matrix (arma::mat) */ -arma::vec rotate_stress(const arma::vec &sigma, const double &angle, const int &axis, const bool &active = true); +arma::mat rotate_strain_concentration(const arma::mat &A, const arma::mat &R, const bool &active = true); /** - * @brief Rotates a stress vector using a 3x3 rotation matrix. - * @param sigma The 6-component stress vector (arma::vec) - * @param R The 3x3 rotation matrix (arma::mat) - * @param active If true (default), active rotation (bool) - * @return The rotated stress vector (arma::vec) - */ -arma::vec rotate_stress(const arma::vec &sigma, const arma::mat &R, const bool &active = true); - -/** - * @brief Rotates a strain vector in Voigt notation. - * @param epsilon The 6-component strain vector (arma::vec) + * @brief Rotates a 6x6 stress concentration matrix by an angle around a given axis. + * @param B The 6x6 stress concentration matrix (arma::mat) * @param angle The rotation angle in radians (double) * @param axis The axis of rotation: 1=x, 2=y, 3=z (int) * @param active If true (default), active rotation (bool) - * @return The rotated strain vector (arma::vec) + * @return The rotated 6x6 stress concentration matrix (arma::mat) */ -arma::vec rotate_strain(const arma::vec &epsilon, const double &angle, const int &axis, const bool &active = true); +arma::mat rotate_stress_concentration(const arma::mat &B, const double &angle, const int &axis, const bool &active = true); /** - * @brief Rotates a strain vector using a 3x3 rotation matrix. - * @param epsilon The 6-component strain vector (arma::vec) + * @brief Rotates a 6x6 stress concentration matrix using a 3x3 rotation matrix. + * @param B The 6x6 stress concentration matrix (arma::mat) * @param R The 3x3 rotation matrix (arma::mat) * @param active If true (default), active rotation (bool) - * @return The rotated strain vector (arma::vec) - */ -arma::vec rotate_strain(const arma::vec &epsilon, const arma::mat &R, const bool &active = true); - -/** - * @brief Rotates a strain tensor from local to global coordinates using Euler angles. - * @param E The 6-component strain vector in local coordinates (arma::vec) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 3x3 strain tensor in global coordinates (arma::mat) - */ -arma::mat rotate_l2g_strain(const arma::vec &E, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a strain tensor from global to local coordinates using Euler angles. - * @param E The 6-component strain vector in global coordinates (arma::vec) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 3x3 strain tensor in local coordinates (arma::mat) - */ -arma::mat rotate_g2l_strain(const arma::vec &E, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a stress tensor from local to global coordinates using Euler angles. - * @param sigma The 6-component stress vector in local coordinates (arma::vec) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 3x3 stress tensor in global coordinates (arma::mat) - */ -arma::mat rotate_l2g_stress(const arma::vec &sigma, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a stress tensor from global to local coordinates using Euler angles. - * @param sigma The 6-component stress vector in global coordinates (arma::vec) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 3x3 stress tensor in local coordinates (arma::mat) - */ -arma::mat rotate_g2l_stress(const arma::vec &sigma, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a stiffness matrix from local to global coordinates using Euler angles. - * @param L The 6x6 stiffness matrix in local coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 stiffness matrix in global coordinates (arma::mat) - */ -arma::mat rotate_l2g_L(const arma::mat &L, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a stiffness matrix from global to local coordinates using Euler angles. - * @param L The 6x6 stiffness matrix in global coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 stiffness matrix in local coordinates (arma::mat) - */ -arma::mat rotate_g2l_L(const arma::mat &L, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a strain localization matrix from local to global coordinates using Euler angles. - * @param A The 6x6 strain localization matrix in local coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 strain localization matrix in global coordinates (arma::mat) - */ -arma::mat rotate_l2g_A(const arma::mat &A, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a strain localization matrix from global to local coordinates using Euler angles. - * @param A The 6x6 strain localization matrix in global coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 strain localization matrix in local coordinates (arma::mat) - */ -arma::mat rotate_g2l_A(const arma::mat &A, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a stress localization matrix from local to global coordinates using Euler angles. - * @param B The 6x6 stress localization matrix in local coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 stress localization matrix in global coordinates (arma::mat) - */ -arma::mat rotate_l2g_B(const arma::mat &B, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a stress localization matrix from global to local coordinates using Euler angles. - * @param B The 6x6 stress localization matrix in global coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 stress localization matrix in local coordinates (arma::mat) - */ -arma::mat rotate_g2l_B(const arma::mat &B, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a compliance matrix from local to global coordinates using Euler angles. - * @param M The 6x6 compliance matrix in local coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 compliance matrix in global coordinates (arma::mat) - */ -arma::mat rotate_l2g_M(const arma::mat &M, const double &psi, const double &theta, const double &phi); - -/** - * @brief Rotates a compliance matrix from global to local coordinates using Euler angles. - * @param M The 6x6 compliance matrix in global coordinates (arma::mat) - * @param psi First Euler angle in radians (double) - * @param theta Second Euler angle in radians (double) - * @param phi Third Euler angle in radians (double) - * @return The 6x6 compliance matrix in local coordinates (arma::mat) + * @return The rotated 6x6 stress concentration matrix (arma::mat) */ -arma::mat rotate_g2l_M(const arma::mat &M, const double &psi, const double &theta, const double &phi); - - +arma::mat rotate_stress_concentration(const arma::mat &B, const arma::mat &R, const bool &active = true); /** @} */ // end of maths group diff --git a/pyproject.toml b/pyproject.toml index 20defb5a3..c65bfec0c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,6 +21,7 @@ maintainers = [ requires-python = ">=3.10" dependencies = [ "numpy>=2.0", + "scipy>=1.10", ] classifiers = [ "Development Status :: 5 - Production/Stable", diff --git a/python-setup/simcoon/__init__.py b/python-setup/simcoon/__init__.py index 266fd2118..09d94462c 100755 --- a/python-setup/simcoon/__init__.py +++ b/python-setup/simcoon/__init__.py @@ -1,5 +1,6 @@ from simcoon._core import * from simcoon.__version__ import __version__ +from simcoon.rotation import Rotation # override _CppRotation from star-import # Backward compatibility alias - simmit was the legacy module name from simcoon import _core as simmit diff --git a/python-setup/simcoon/rotation.py b/python-setup/simcoon/rotation.py new file mode 100644 index 000000000..dea3d04f1 --- /dev/null +++ b/python-setup/simcoon/rotation.py @@ -0,0 +1,310 @@ +""" +Rotation class extending scipy.spatial.transform.Rotation with mechanics operations. + +This module provides a ``Rotation`` class that inherits from scipy's ``Rotation``, +giving users access to all scipy rotation features (batch operations, ``mean()``, +``Slerp``, ``RotationSpline``, etc.) while adding simcoon's continuum-mechanics +methods (``apply_stress``, ``apply_stiffness``, etc.). +""" + +import numpy as np +from scipy.spatial.transform import Rotation as ScipyRotation + +from simcoon._core import _CppRotation + + +class Rotation(ScipyRotation): + """Rotation class extending scipy with continuum-mechanics operations. + + Inherits all scipy.spatial.transform.Rotation functionality (batch rotations, + ``mean()``, ``Slerp``, ``RotationSpline``, ``from_euler``, ``from_quat``, etc.) + and adds methods for rotating Voigt-notation stress, strain, stiffness, and + compliance tensors via the C++ backend. + + Examples + -------- + >>> import simcoon as smc + >>> import numpy as np + + >>> # All scipy factory methods return a simcoon Rotation + >>> r = smc.Rotation.from_rotvec([0, 0, np.pi/4]) + >>> isinstance(r, smc.Rotation) + True + + >>> # Rotate a stiffness matrix + >>> L = smc.L_iso([70000, 0.3], 'Enu') + >>> L_rot = r.apply_stiffness(L) + + >>> # Simcoon-specific: create rotation around a principal axis + >>> r = smc.Rotation.from_axis_angle(np.pi/4, 3) # 45 deg around z + """ + + # ------------------------------------------------------------------ + # Simcoon-specific factory methods + # ------------------------------------------------------------------ + + @classmethod + def from_axis_angle(cls, angle, axis, degrees=False): + """Create a rotation around a principal axis. + + Parameters + ---------- + angle : float + Rotation angle. + axis : int + Axis of rotation: 1 = x, 2 = y, 3 = z. + degrees : bool, optional + If True, *angle* is in degrees. Default is False (radians). + + Returns + ------- + Rotation + """ + if axis not in (1, 2, 3): + raise ValueError(f"axis must be 1, 2, or 3, got {axis}") + direction = np.zeros(3) + direction[axis - 1] = 1.0 + return cls.from_rotvec(direction * (angle if not degrees else np.radians(angle))) + + @classmethod + def from_scipy(cls, scipy_rot): + """Create a simcoon Rotation from a plain scipy Rotation. + + This is a convenience method for upgrading a + ``scipy.spatial.transform.Rotation`` to a ``simcoon.Rotation`` so + that the mechanics methods (``apply_stress``, ``apply_stiffness``, + etc.) become available. + + Parameters + ---------- + scipy_rot : scipy.spatial.transform.Rotation + A scipy Rotation object (single or batch). + + Returns + ------- + Rotation + """ + return cls.from_quat(scipy_rot.as_quat()) + + # ------------------------------------------------------------------ + # Internal helper + # ------------------------------------------------------------------ + + def _to_cpp(self): + """Convert to a _CppRotation for C++ method dispatch.""" + return _CppRotation.from_quat(self.as_quat()) + + # ------------------------------------------------------------------ + # Mechanics methods (delegate to _CppRotation) + # ------------------------------------------------------------------ + + def apply_stress(self, sigma, active=True): + """Apply rotation to a stress vector in Voigt notation. + + Parameters + ---------- + sigma : array_like + 6-component stress vector [s11, s22, s33, s12, s13, s23]. + active : bool, optional + If True (default), active rotation. + + Returns + ------- + numpy.ndarray + Rotated stress vector. + """ + return self._to_cpp().apply_stress(np.asarray(sigma, dtype=float), active).ravel() + + def apply_strain(self, epsilon, active=True): + """Apply rotation to a strain vector in Voigt notation. + + Parameters + ---------- + epsilon : array_like + 6-component strain vector [e11, e22, e33, 2*e12, 2*e13, 2*e23]. + active : bool, optional + If True (default), active rotation. + + Returns + ------- + numpy.ndarray + Rotated strain vector. + """ + return self._to_cpp().apply_strain(np.asarray(epsilon, dtype=float), active).ravel() + + def apply_stiffness(self, L, active=True): + """Apply rotation to a 6x6 stiffness matrix. + + Parameters + ---------- + L : array_like + 6x6 stiffness matrix in Voigt notation. + active : bool, optional + If True (default), active rotation. + + Returns + ------- + numpy.ndarray + Rotated 6x6 stiffness matrix. + """ + return self._to_cpp().apply_stiffness(np.asarray(L, dtype=float), active) + + def apply_compliance(self, M, active=True): + """Apply rotation to a 6x6 compliance matrix. + + Parameters + ---------- + M : array_like + 6x6 compliance matrix in Voigt notation. + active : bool, optional + If True (default), active rotation. + + Returns + ------- + numpy.ndarray + Rotated 6x6 compliance matrix. + """ + return self._to_cpp().apply_compliance(np.asarray(M, dtype=float), active) + + def apply_strain_concentration(self, A, active=True): + """Apply rotation to a 6x6 strain concentration tensor. + + Parameters + ---------- + A : array_like + 6x6 strain concentration tensor in Voigt notation. + active : bool, optional + If True (default), active rotation. + + Returns + ------- + numpy.ndarray + Rotated strain concentration tensor: QE * A * QS^T. + """ + return self._to_cpp().apply_strain_concentration(np.asarray(A, dtype=float), active) + + def apply_stress_concentration(self, B, active=True): + """Apply rotation to a 6x6 stress concentration tensor. + + Parameters + ---------- + B : array_like + 6x6 stress concentration tensor in Voigt notation. + active : bool, optional + If True (default), active rotation. + + Returns + ------- + numpy.ndarray + Rotated stress concentration tensor: QS * B * QE^T. + """ + return self._to_cpp().apply_stress_concentration(np.asarray(B, dtype=float), active) + + def apply_tensor(self, m, inverse=False): + """Apply rotation to a 3x3 tensor (matrix). + + Parameters + ---------- + m : array_like + 3x3 tensor to rotate. + inverse : bool, optional + If True, apply inverse rotation. Default is False. + + Returns + ------- + numpy.ndarray + Rotated tensor: R * m * R^T (or R^T * m * R for inverse). + """ + return self._to_cpp().apply_tensor(np.asarray(m, dtype=float), inverse) + + def as_voigt_stress_rotation(self, active=True): + """Get 6x6 rotation matrix for stress tensors in Voigt notation. + + Parameters + ---------- + active : bool, optional + If True (default), active rotation; if False, passive. + + Returns + ------- + numpy.ndarray + 6x6 stress rotation matrix (QS). + """ + return self._to_cpp().as_voigt_stress_rotation(active) + + def as_voigt_strain_rotation(self, active=True): + """Get 6x6 rotation matrix for strain tensors in Voigt notation. + + Parameters + ---------- + active : bool, optional + If True (default), active rotation; if False, passive. + + Returns + ------- + numpy.ndarray + 6x6 strain rotation matrix (QE). + """ + return self._to_cpp().as_voigt_strain_rotation(active) + + # ------------------------------------------------------------------ + # Compatibility helpers + # ------------------------------------------------------------------ + + def equals(self, other, tol=1e-12): + """Check if this rotation equals another within tolerance. + + Accounts for quaternion antipodal equivalence (q and -q represent + the same rotation). + + Parameters + ---------- + other : Rotation + Rotation to compare with. + tol : float, optional + Tolerance for comparison. Default is 1e-12. + + Returns + ------- + bool + """ + q1 = self.as_quat() + q2 = other.as_quat() + return np.allclose(q1, q2, atol=tol) or np.allclose(q1, -q2, atol=tol) + + def is_identity(self, tol=1e-12): + """Check if this rotation is identity (no rotation). + + Parameters + ---------- + tol : float, optional + Tolerance for comparison. Default is 1e-12. + + Returns + ------- + bool + """ + return self.magnitude() < tol + + def slerp_to(self, other, t): + """Spherical linear interpolation between this rotation and *other*. + + For full-featured interpolation over multiple keyframes, use + ``scipy.spatial.transform.Slerp`` instead. + + Parameters + ---------- + other : Rotation + Target rotation. + t : float + Interpolation parameter in [0, 1] (0 = self, 1 = other). + + Returns + ------- + Rotation + Interpolated rotation. + """ + from scipy.spatial.transform import Slerp + key_rots = type(self).concatenate([self, other]) + interp = Slerp([0.0, 1.0], key_rots) + return interp(t) diff --git a/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Maths/rotation.hpp b/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Maths/rotation.hpp old mode 100755 new mode 100644 index f97f961a8..100308d62 --- a/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Maths/rotation.hpp +++ b/simcoon-python-builder/include/simcoon/python_wrappers/Libraries/Maths/rotation.hpp @@ -2,72 +2,9 @@ #include #include -namespace simpy{ +namespace simpy { -//This function returns a rotated vector (3) according to a rotation matrix -pybind11::array_t rotate_vec_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool ©=true); +// Function to register the Rotation class with pybind11 +void register_rotation(pybind11::module_& m); -//This function returns a rotated vector (3) according to an angle and an axis -pybind11::array_t rotate_vec_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool ©=true); - -//This function returns a rotated matrix (3x3) according to a rotation matrix -pybind11::array_t rotate_mat_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool ©=true); - -//This function returns a rotated matrix (3x3) according to an angle and an axis -pybind11::array_t rotate_mat_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool ©=true); - -//This function returns the 3*3 rotation matrix according to an angle, an axis and depending if it is active or passive rotation -pybind11::array_t fillR_angle(const double &angle, const int &axis, const bool &active=true, const bool ©=true); - -//This function returns the 3*3 rotation matrix according to the three Euler angles, depending if it is active or passive rotation and the Euler convention (ex :"zxz") -pybind11::array_t fillR_euler(const double &psi, const double &theta, const double &phi, const bool &active=true, const std::string &conv="zxz", const bool ©=true); - -//This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'stress' from an angle and an axis -pybind11::array_t fillQS_angle(const double &angle, const int &axis, const bool &active=true, const bool ©=true); - -//This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'stress' from a rotation matrix -pybind11::array_t fillQS_R(const pybind11::array_t &R, const bool &active=true, const bool ©=true); - -//This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'strain' from an angle and an axis -pybind11::array_t fillQE_angle(const double &angle, const int &axis, const bool &active=true, const bool ©=true); - -//This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'strain' from a rotation matrix -pybind11::array_t fillQE_R(const pybind11::array_t &R, const bool &active=true, const bool ©=true); - -//These functions rotates a 6*6 stiffness arma::matrix (L) -pybind11::array_t rotateL_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool &active=true, const bool ©=true); -pybind11::array_t rotateL_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_l2g_L(const pybind11::array_t &input, const double &, const double &, const double &, const bool ©=true); -pybind11::array_t rotate_g2l_L(const pybind11::array_t &input, const double &, const double &, const double &, const bool ©=true); - -//These functions rotates a 6*6 compliance arma::matrix (M) -pybind11::array_t rotateM_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool &active=true, const bool ©=true); -pybind11::array_t rotateM_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_l2g_M(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); -pybind11::array_t rotate_g2l_M(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); - -//These functions rotates a 6*6 strain concentration (A) -pybind11::array_t rotateA_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool &active=true, const bool ©=true); -pybind11::array_t rotateA_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_l2g_A(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); -pybind11::array_t rotate_g2l_A(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); - -//These functions rotates a 6*6 stress concentration (B) -pybind11::array_t rotateB_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool &active=true, const bool ©=true); -pybind11::array_t rotateB_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_l2g_B(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); -pybind11::array_t rotate_g2l_B(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); - -//These functions rotates strain arma::vectors -pybind11::array_t rotate_strain_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_strain_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_l2g_strain(const pybind11::array_t &, const double &psi, const double &theta, const double &phi, const bool ©=true); -pybind11::array_t rotate_g2l_strain(const pybind11::array_t &, const double &psi, const double &theta, const double &phi, const bool ©=true); - -//These functions rotates stress arma::vectors -pybind11::array_t rotate_stress_angle(const pybind11::array_t &input, const double &angle, const int &axis, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_stress_R(const pybind11::array_t &input, const pybind11::array_t &R, const bool &active=true, const bool ©=true); -pybind11::array_t rotate_l2g_stress(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); -pybind11::array_t rotate_g2l_stress(const pybind11::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©=true); - -} //namespace simpy \ No newline at end of file +} // namespace simpy diff --git a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp old mode 100755 new mode 100644 index efa223b56..f3b0f3eab --- a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp +++ b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -10,307 +11,140 @@ using namespace std; using namespace arma; -namespace py=pybind11; +namespace py = pybind11; namespace simpy { -py::array_t rotate_vec_R(const py::array_t &input, const py::array_t &R, const bool ©) { - vec v_cpp = carma::arr_to_col(input); - mat R_cpp = carma::arr_to_mat(R); - vec rotated_v = simcoon::rotate_vec(v_cpp,R_cpp); - return carma::col_to_arr(rotated_v, copy); -} - -py::array_t rotate_vec_angle(const py::array_t &input, const double &angle, const int &axis, const bool ©) { - vec v_cpp = carma::arr_to_col(input); - vec rotated_v = simcoon::rotate_vec(v_cpp,angle,axis); - return carma::col_to_arr(rotated_v, copy); -} - -py::array_t rotate_mat_R(const py::array_t &input, const py::array_t &R, const bool ©) { - mat m_cpp = carma::arr_to_mat(input); - mat R_cpp = carma::arr_to_mat(R); - mat rotated_m = simcoon::rotate_mat(m_cpp,R_cpp); - return carma::mat_to_arr(rotated_m, copy); -} - -py::array_t rotate_mat_angle(const py::array_t &input, const double &angle, const int &axis, const bool ©) { - mat m_cpp = carma::arr_to_mat(input); - mat rotated_m = simcoon::rotate_mat(m_cpp,angle,axis); - return carma::mat_to_arr(rotated_m, copy); -} - -py::array_t fillR_angle(const double &angle, const int &axis, const bool &active, const bool ©) { - mat R = simcoon::fillR(angle,axis,active); - return carma::mat_to_arr(R, copy); -} - -//This function returns the 3*3 rotation matrix -py::array_t fillR_euler(const double &psi, const double &theta, const double &phi, const bool &active, const string &conv, const bool ©) { - mat R = simcoon::fillR(psi,theta,phi,active,conv); - return carma::mat_to_arr(R, copy); -} - -//This function returns the 6*6 rotation matrix of a vector of type 'stress' -py::array_t fillQS_angle(const double &angle, const int &axis, const bool &active, const bool ©) { - mat QS = simcoon::fillQS(angle,axis,active); - return carma::mat_to_arr(QS, copy); -} - -//This function returns the 6*6 rotation matrix of a vector of type 'stress' -py::array_t fillQS_R(const py::array_t &R, const bool &active, const bool ©) { - mat R_cpp = carma::arr_to_mat(R); - mat QS = simcoon::fillQS(R_cpp,active); - return carma::mat_to_arr(QS, copy); -} - -//This function returns the 6*6 rotation matrix of a vector of type 'strain' -py::array_t fillQE_angle(const double &angle, const int &axis, const bool &active, const bool ©) { - mat QE = simcoon::fillQE(angle,axis,active); - return carma::mat_to_arr(QE, copy); -} - -//This function returns the 6*6 rotation matrix of a vector of type 'strain' -py::array_t fillQE_R(const py::array_t &R, const bool &active, const bool ©) { - mat R_cpp = carma::arr_to_mat(R); - mat QE = simcoon::fillQE(R_cpp,active); - return carma::mat_to_arr(QE, copy); -} - -//This function rotates a 6*6 stiffness matrix (L) -py::array_t rotateL_angle(const py::array_t &input, const double &angle, const int & axis, const bool &active, const bool ©) { - mat L_cpp = carma::arr_to_mat(input); - mat rotated_L = simcoon::rotateL(L_cpp,angle,axis,active); - return carma::mat_to_arr(rotated_L, copy); -} - -//This function rotates a 6*6 stiffness matrix (L) -py::array_t rotateL_R(const py::array_t &input, const py::array_t &R, const bool &active, const bool ©) { - mat L_cpp = carma::arr_to_mat(input); - mat R_cpp = carma::arr_to_mat(R); - mat rotated_L = simcoon::rotateL(L_cpp,R_cpp,active); - return carma::mat_to_arr(rotated_L, copy); -} - -//This function rotates a 6*6 compliance matrix (M) -py::array_t rotateM_angle(const py::array_t &input, const double &angle, const int &axis, const bool &active, const bool ©) { - mat M_cpp = carma::arr_to_mat(input); - mat rotated_M = simcoon::rotateM(M_cpp,angle,axis,active); - return carma::mat_to_arr(rotated_M, copy); -} - -//This function rotates a 6*6 compliance matrix (M) -py::array_t rotateM_R(const py::array_t &input, const py::array_t &R, const bool &active, const bool ©) { - mat M_cpp = carma::arr_to_mat(input); - mat R_cpp = carma::arr_to_mat(R); - mat rotated_M = simcoon::rotateM(M_cpp,R_cpp,active); - return carma::mat_to_arr(rotated_M); -} - -//This function rotates a 6*6 strain concentration (A) -py::array_t rotateA_angle(const py::array_t &input, const double &angle, const int &axis, const bool &active, const bool ©) { - mat A_cpp = carma::arr_to_mat(input); - mat rotated_A = simcoon::rotateA(A_cpp,angle,axis,active); - return carma::mat_to_arr(rotated_A, copy); -} - -//This function rotates a 6*6 strain concentration (A) -py::array_t rotateA_R(const py::array_t &input, const py::array_t &R, const bool &active, const bool ©) { - mat A_cpp = carma::arr_to_mat(input); - mat R_cpp = carma::arr_to_mat(R); - mat rotated_A = simcoon::rotateA(A_cpp,R_cpp,active); - return carma::mat_to_arr(rotated_A, copy); -} - -//This function rotates a 6*6 stress concentration (B) -py::array_t rotateB_angle(const py::array_t &input, const double &angle, const int &axis, const bool &active, const bool ©) { - mat B_cpp = carma::arr_to_mat(input); - mat rotated_B = simcoon::rotateB(B_cpp,angle,axis,active); - return carma::mat_to_arr(rotated_B, copy); -} - -//This function rotates a 6*6 stress concentration (B) -py::array_t rotateB_R(const py::array_t &input, const py::array_t &R, const bool &active, const bool ©) { - mat B_cpp = carma::arr_to_mat(input); - mat R_cpp = carma::arr_to_mat(R); - mat rotated_B = simcoon::rotateB(B_cpp,R_cpp,active); - return carma::mat_to_arr(rotated_B, copy); -} - -//This function rotates stress vectors - Can be used with stack of arrays (vectorized) -py::array_t rotate_stress_angle(const py::array_t &input, const double &angle, const int &axis, const bool &active, const bool ©) { - if (input.ndim()==1){ - vec v = carma::arr_to_col(input); - vec rotated_stress = simcoon::rotate_stress(v,angle,axis,active); - return carma::col_to_arr(rotated_stress, copy); - } - else if (input.ndim() == 2){ - mat m = carma::arr_to_mat_view(input); - int nb_points = m.n_cols; - mat rotated_stress(6,nb_points); - for (int pt = 0; pt < nb_points; pt++) { - vec v = m.unsafe_col(pt); - rotated_stress.col(pt) = simcoon::rotate_stress(v,angle,axis,active); +namespace { + void validate_vector_size(const py::array_t& arr, size_t expected, const string& name) { + if (arr.ndim() != 1) { + throw invalid_argument(name + " must be a 1D array, got " + to_string(arr.ndim()) + "D"); } - return carma::mat_to_arr(rotated_stress, copy); - } - else{ - throw std::invalid_argument("input.ndim should be 1 or 2"); - } -} - -//This function rotates stress vectors - Can be used with stack of arrays (vectorized) -py::array_t rotate_stress_R(const py::array_t &input, const py::array_t &R, const bool &active, const bool ©) { - if (input.ndim()==1){ - vec v = carma::arr_to_col(input); - mat R_cpp = carma::arr_to_mat(R); - vec rotated_stress = simcoon::rotate_stress(v,R_cpp,active); - return carma::col_to_arr(rotated_stress, copy); - } - else if (input.ndim() == 2){ - mat m = carma::arr_to_mat_view(input); - cube R_cpp = carma::arr_to_cube_view(R); - int nb_points = m.n_cols; - mat rotated_stress(6,nb_points); - for (int pt = 0; pt < nb_points; pt++) { - vec v = m.unsafe_col(pt); - rotated_stress.col(pt) = simcoon::rotate_stress(v,R_cpp.slice(pt),active); + if (static_cast(arr.size()) != expected) { + throw invalid_argument(name + " must have " + to_string(expected) + + " elements, got " + to_string(arr.size())); } - return carma::mat_to_arr(rotated_stress, copy); - } - else{ - throw std::invalid_argument("input.ndim should be 1 or 2"); } -} -//This function rotates strain vectors - Can be used with stack of arrays (vectorized) -py::array_t rotate_strain_angle(const py::array_t &input, const double &angle, const int &axis, const bool &active, const bool ©) { - if (input.ndim()==1){ - vec v = carma::arr_to_col(input); - vec rotated_strain = simcoon::rotate_strain(v,angle,axis,active); - return carma::col_to_arr(rotated_strain, copy); - } - else if (input.ndim() == 2){ - mat m = carma::arr_to_mat_view(input); - int nb_points = m.n_cols; - mat rotated_strain(6,nb_points); - for (int pt = 0; pt < nb_points; pt++) { - vec v = m.unsafe_col(pt); - rotated_strain.col(pt) = simcoon::rotate_strain(v,angle,axis,active); + void validate_matrix_size(const py::array_t& arr, size_t rows, size_t cols, const string& name) { + if (arr.ndim() != 2) { + throw invalid_argument(name + " must be a 2D array, got " + to_string(arr.ndim()) + "D"); } - return carma::mat_to_arr(rotated_strain, copy); - } - else{ - throw std::invalid_argument("input.ndim should be 1 or 2"); - } -} - -//This function rotates strain vectors - Can be used with stack of arrays (vectorized) -py::array_t rotate_strain_R(const py::array_t &input, const py::array_t &R, const bool &active, const bool ©) { - if (input.ndim()==1){ - vec v = carma::arr_to_col(input); - mat R_cpp = carma::arr_to_mat(R); - vec rotated_strain = simcoon::rotate_strain(v,R_cpp,active); - return carma::col_to_arr(rotated_strain, copy); - } - else if (input.ndim() == 2){ - mat m = carma::arr_to_mat_view(input); - cube R_cpp = carma::arr_to_cube_view(R); - int nb_points = m.n_cols; - mat rotated_strain(6,nb_points); - for (int pt = 0; pt < nb_points; pt++) { - vec v = m.unsafe_col(pt); - rotated_strain.col(pt) = simcoon::rotate_strain(v,R_cpp.slice(pt),active); + auto shape = arr.shape(); + if (static_cast(shape[0]) != rows || static_cast(shape[1]) != cols) { + throw invalid_argument(name + " must have shape (" + to_string(rows) + ", " + + to_string(cols) + "), got (" + to_string(shape[0]) + ", " + to_string(shape[1]) + ")"); } - return carma::mat_to_arr(rotated_strain, copy); - } - else{ - throw std::invalid_argument("input.ndim should be 1 or 2"); } -} - -//This function rotates strain vectors from a local to global set of coordinates (using Euler angles) -py::array_t rotate_l2g_strain(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - vec E_cpp = carma::arr_to_col(input); - vec rotated_E = simcoon::rotate_l2g_strain(E_cpp,psi,theta,phi); - return carma::col_to_arr(rotated_E, copy); -} - -//This function rotates strain vectors from a global to local set of coordinates (using Euler angles) -py::array_t rotate_g2l_strain(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - vec E_cpp = carma::arr_to_col(input); - vec rotated_E = simcoon::rotate_g2l_strain(E_cpp,psi,theta,phi); - return carma::col_to_arr(rotated_E, copy); -} - -//This function rotates stress vectors from a local to global set of coordinates (using Euler angles) -py::array_t rotate_l2g_stress(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - vec S_cpp = carma::arr_to_col(input); - vec rotated_S = simcoon::rotate_l2g_strain(S_cpp,psi,theta,phi); - return carma::col_to_arr(rotated_S, copy); -} - -//This function rotates stress vectors from a global to local set of coordinates (using Euler angles) -py::array_t rotate_g2l_stress(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - vec S_cpp = carma::arr_to_col(input); - vec rotated_S = simcoon::rotate_g2l_stress(S_cpp,psi,theta,phi); - return carma::col_to_arr(rotated_S, copy); -} - -//This function rotates stiffness matrices (L) from a local to global set of coordinates (using Euler angles) -py::array_t rotate_l2g_L(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat L_cpp = carma::arr_to_mat(input); - mat rotated_L = simcoon::rotate_l2g_L(L_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_L, copy); -} - -//This function rotates stiffness matrices (L) from a global to local set of coordinates (using Euler angles) -py::array_t rotate_g2l_L(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat L_cpp = carma::arr_to_mat(input); - mat rotated_L = simcoon::rotate_g2l_L(L_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_L, copy); -} - -//This function rotates compliance matrices (M) from a local to global set of coordinates (using Euler angles) -py::array_t rotate_l2g_M(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat M_cpp = carma::arr_to_mat(input); - mat rotated_M = simcoon::rotate_l2g_M(M_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_M, copy); -} - -//This function rotates compliance matrices (M) from a global to local set of coordinates (using Euler angles) -py::array_t rotate_g2l_M(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat M_cpp = carma::arr_to_mat(input); - mat rotated_M = simcoon::rotate_g2l_M(M_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_M, copy); -} - -//This function rotates strain concentration matrices (A) from a local to global set of coordinates (using Euler angles) -py::array_t rotate_l2g_A(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat A_cpp = carma::arr_to_mat(input); - mat rotated_A = simcoon::rotate_l2g_A(A_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_A, copy); -} - -//This function rotates strain concentration matrices (A) from a global to local set of coordinates (using Euler angles) -py::array_t rotate_g2l_A(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat A_cpp = carma::arr_to_mat(input); - mat rotated_A = simcoon::rotate_g2l_A(A_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_A, copy); -} - -//This function rotates stress concentration matrices (B) from a local to global set of coordinates (using Euler angles) -py::array_t rotate_l2g_B(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat B_cpp = carma::arr_to_mat(input); - mat rotated_B = simcoon::rotate_l2g_B(B_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_B, copy); -} - -//This function rotates stress concentration matrices (B) from a global to local set of coordinates (using Euler angles) -py::array_t rotate_g2l_B(const py::array_t &input, const double &psi, const double &theta, const double &phi, const bool ©) { - mat B_cpp = carma::arr_to_mat(input); - mat rotated_B = simcoon::rotate_g2l_B(B_cpp,psi,theta,phi); - return carma::mat_to_arr(rotated_B, copy); -} - -} //namepsace simpy +} // anonymous namespace + +void register_rotation(py::module_& m) { + py::class_(m, "_CppRotation", + R"doc( + Internal C++ rotation backend using unit quaternions (scalar-last). + + End users should use ``simcoon.Rotation`` instead, which inherits from + ``scipy.spatial.transform.Rotation`` and delegates mechanics operations + to this class. + )doc") + + // The only factory method needed — Python Rotation._to_cpp() uses this + .def_static("from_quat", + [](py::array_t quat) { + validate_vector_size(quat, 4, "quat"); + vec q = carma::arr_to_col(quat); + return simcoon::Rotation::from_quat(q); + }, + py::arg("quat"), + "Create rotation from quaternion [qx, qy, qz, qw] (scalar-last)") + + // Voigt rotation matrices + .def("as_voigt_stress_rotation", + [](const simcoon::Rotation& self, bool active) { + return carma::mat_to_arr(mat(self.as_voigt_stress_rotation(active))); + }, + py::arg("active") = true, + "Get 6x6 rotation matrix for stress tensors in Voigt notation") + + .def("as_voigt_strain_rotation", + [](const simcoon::Rotation& self, bool active) { + return carma::mat_to_arr(mat(self.as_voigt_strain_rotation(active))); + }, + py::arg("active") = true, + "Get 6x6 rotation matrix for strain tensors in Voigt notation") + + // Apply methods — the core mechanics operations + .def("apply_tensor", + [](const simcoon::Rotation& self, py::array_t m, bool inverse) { + validate_matrix_size(m, 3, 3, "m"); + mat m_cpp = carma::arr_to_mat(m); + mat result = self.apply_tensor(m_cpp, inverse); + return carma::mat_to_arr(result); + }, + py::arg("m"), py::arg("inverse") = false, + "Apply rotation to a 3x3 tensor: R * m * R^T") + + .def("apply_stress", + [](const simcoon::Rotation& self, py::array_t sigma, bool active) { + validate_vector_size(sigma, 6, "sigma"); + vec sigma_cpp = carma::arr_to_col(sigma); + vec result = self.apply_stress(sigma_cpp, active); + return carma::col_to_arr(result); + }, + py::arg("sigma"), py::arg("active") = true, + "Apply rotation to a 6-component stress vector in Voigt notation") + + .def("apply_strain", + [](const simcoon::Rotation& self, py::array_t epsilon, bool active) { + validate_vector_size(epsilon, 6, "epsilon"); + vec epsilon_cpp = carma::arr_to_col(epsilon); + vec result = self.apply_strain(epsilon_cpp, active); + return carma::col_to_arr(result); + }, + py::arg("epsilon"), py::arg("active") = true, + "Apply rotation to a 6-component strain vector in Voigt notation") + + .def("apply_stiffness", + [](const simcoon::Rotation& self, py::array_t L, bool active) { + validate_matrix_size(L, 6, 6, "L"); + mat L_cpp = carma::arr_to_mat(L); + mat result = self.apply_stiffness(L_cpp, active); + return carma::mat_to_arr(result); + }, + py::arg("L"), py::arg("active") = true, + "Apply rotation to a 6x6 stiffness matrix") + + .def("apply_compliance", + [](const simcoon::Rotation& self, py::array_t M, bool active) { + validate_matrix_size(M, 6, 6, "M"); + mat M_cpp = carma::arr_to_mat(M); + mat result = self.apply_compliance(M_cpp, active); + return carma::mat_to_arr(result); + }, + py::arg("M"), py::arg("active") = true, + "Apply rotation to a 6x6 compliance matrix") + + .def("apply_strain_concentration", + [](const simcoon::Rotation& self, py::array_t A, bool active) { + validate_matrix_size(A, 6, 6, "A"); + mat A_cpp = carma::arr_to_mat(A); + mat result = self.apply_strain_concentration(A_cpp, active); + return carma::mat_to_arr(result); + }, + py::arg("A"), py::arg("active") = true, + "Apply rotation to a 6x6 strain concentration tensor") + + .def("apply_stress_concentration", + [](const simcoon::Rotation& self, py::array_t B, bool active) { + validate_matrix_size(B, 6, 6, "B"); + mat B_cpp = carma::arr_to_mat(B); + mat result = self.apply_stress_concentration(B_cpp, active); + return carma::mat_to_arr(result); + }, + py::arg("B"), py::arg("active") = true, + "Apply rotation to a 6x6 stress concentration tensor") + + ; +} + +} // namespace simpy diff --git a/simcoon-python-builder/src/python_wrappers/python_module.cpp b/simcoon-python-builder/src/python_wrappers/python_module.cpp index 04489d88b..670ebae83 100755 --- a/simcoon-python-builder/src/python_wrappers/python_module.cpp +++ b/simcoon-python-builder/src/python_wrappers/python_module.cpp @@ -208,41 +208,8 @@ PYBIND11_MODULE(_core, m) m.def("Eshelby", &Eshelby, "L"_a, "a1"_a = 1., "a2"_a = 1., "a3"_a = 1., "mp"_a = 50, "np"_a = 50, "copy"_a = true, simcoon_docs::Eshelby); m.def("T_II", &T_II, "L"_a, "a1"_a = 1., "a2"_a = 1., "a3"_a = 1., "mp"_a = 50, "np"_a = 50, "copy"_a = true, simcoon_docs::T_II); - // Register the rotation library - m.def("rotate_vec_R", &rotate_vec_R, "input"_a, "R"_a, "copy"_a = true, "This function returns a rotated vector (3) according to a rotation matrix"); - m.def("rotate_vec_angle", &rotate_vec_angle, "input"_a, "angle"_a, "axis"_a, "copy"_a = true, "This function returns a rotated vector (3) according to an angle and an axis"); - m.def("rotate_mat_R", &rotate_mat_R, "input"_a, "R"_a, "copy"_a = true, "This function returns a rotated matrix (3x3) according to a rotation matrix"); - m.def("rotate_mat_angle", &rotate_mat_angle, "input"_a, "angle"_a, "axis"_a, "copy"_a = true, "This function returns a rotated matrix (3x3) according to an angle and an axis"); - m.def("fillR_angle", &fillR_angle, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = true, "This function returns the 3*3 rotation matrix according to an angle, an axis and depending if it is active or passive rotation"); - m.def("fillR_euler", &fillR_euler, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "conv"_a = "zxz", "copy"_a = true, "This function returns the 3*3 rotation matrix according to the three Euler angles, depending if it is active or passive rotation and the Euler convention (ex :zxz)"); - m.def("fillQS_angle", &fillQS_angle, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = true, "This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'stress' from an angle and an axis"); - m.def("fillQS_R", &fillQS_R, "R"_a, "active"_a = true, "copy"_a = true, "This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'stress' from a rotation matrix"); - m.def("fillQE_angle", &fillQE_angle, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = true, "This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'strain' from an angle and an axis"); - m.def("fillQE_R", &fillQE_R, "R"_a, "active"_a = true, "copy"_a = true, "This function returns the 6*6 rotation arma::matrix of a arma::vector of type 'strain' from a rotation matrix"); - m.def("rotateL_angle", &rotateL_angle, "input"_a, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 stiffness matrix according to an angle and an axis"); - m.def("rotateL_R", &rotateL_R, "input"_a, "R"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 stiffness matrix according to a rotation matrix"); - m.def("rotate_l2g_L", &rotate_l2g_L, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 stiffness matrix from local to global frame"); - m.def("rotate_g2l_L", &rotate_g2l_L, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 stiffness matrix from global to local frame"); - m.def("rotateM_angle", &rotateM_angle, "input"_a, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 compliance matrix according to an angle and an axis"); - m.def("rotateM_R", &rotateM_R, "input"_a, "R"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 compliance matrix according to a rotation matrix"); - m.def("rotate_l2g_M", &rotate_l2g_M, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 compliance matrix from local to global frame"); - m.def("rotate_g2l_M", &rotate_g2l_M, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 compliance matrix from global to local frame"); - m.def("rotateA_angle", &rotateA_angle, "input"_a, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 strain concentration matrix according to an angle and an axis"); - m.def("rotateA_R", &rotateA_R, "input"_a, "R"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 strain concentration matrix according to a rotation matrix"); - m.def("rotate_l2g_A", &rotate_l2g_A, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 strain concentration matrix from local to global frame"); - m.def("rotate_g2l_A", &rotate_g2l_A, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 strain concentration matrix from global to local frame"); - m.def("rotateB_angle", &rotateB_angle, "input"_a, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 stress concentration matrix according to an angle and an axis"); - m.def("rotateB_R", &rotateB_R, "input"_a, "R"_a, "active"_a = true, "copy"_a = true, "Return the rotated 6*6 stress concentration matrix according to a rotation matrix"); - m.def("rotate_l2g_B", &rotate_l2g_B, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 stress concentration matrix from local to global frame"); - m.def("rotate_g2l_B", &rotate_g2l_B, "input"_a, "psi"_a, "theta"_a, "phi"_a, "active"_a = true, "Return the rotated 6*6 stress concentration matrix from global to local frame"); - m.def("rotate_strain_angle", &rotate_strain_angle, "input"_a, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = false, "Return the rotated strain matrix using voigt notations according to an angle and an axis"); - m.def("rotate_strain_R", &rotate_strain_R, "input"_a, "R"_a, "active"_a = true, "copy"_a = false, "Return the rotated strain matrix using voigt notations according to a rotation matrix"); - m.def("rotate_l2g_strain", &rotate_l2g_strain, "input"_a, "psi"_a, "theta"_a, "phi"_a, "copy"_a = true, "Return the rotated strain matrix using voigt notations from local to global frame"); - m.def("rotate_g2l_strain", &rotate_g2l_strain, "input"_a, "psi"_a, "theta"_a, "phi"_a, "copy"_a = true, "Return the rotated strain matrix using voigt notations from global to local frame"); - m.def("rotate_stress_angle", &rotate_stress_angle, "input"_a, "angle"_a, "axis"_a, "active"_a = true, "copy"_a = false, "Return the rotated stress matrix using voigt notations according to an angle and an axis"); - m.def("rotate_stress_R", &rotate_stress_R, "input"_a, "R"_a, "active"_a = true, "copy"_a = false, "Return the rotated stress matrix using voigt notations according to a rotation matrix"); - m.def("rotate_l2g_stress", &rotate_l2g_stress, "input"_a, "psi"_a, "theta"_a, "phi"_a, "copy"_a = true, "Return the rotated stress matrix using voigt notations from local to global frame"); - m.def("rotate_g2l_stress", &rotate_g2l_stress, "input"_a, "psi"_a, "theta"_a, "phi"_a, "copy"_a = true, "Return the rotated stress matrix using voigt notations from global to local frame"); + // Register the Rotation class + register_rotation(m); // Register the from-python converters for lagrange m.def("lagrange_exp", &lagrange_exp, "This function is used to determine an exponential Lagrange Multiplier (like contact in Abaqus)"); diff --git a/simcoon-python-builder/test/test_core/test_rotation.py b/simcoon-python-builder/test/test_core/test_rotation.py new file mode 100644 index 000000000..a9f43a50b --- /dev/null +++ b/simcoon-python-builder/test/test_core/test_rotation.py @@ -0,0 +1,293 @@ +"""Tests for the Rotation class and scipy interoperability.""" + +import numpy as np +import pytest +from scipy.spatial.transform import Rotation as ScipyRotation + +import simcoon as sim + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture +def angle(): + return np.pi / 4 + + +@pytest.fixture +def sim_rot(angle): + """A simcoon.Rotation (scipy subclass).""" + return sim.Rotation.from_axis_angle(angle, 3) + + +@pytest.fixture +def scipy_rot(angle): + """A plain scipy Rotation representing the same rotation.""" + return ScipyRotation.from_rotvec([0, 0, angle]) + + +@pytest.fixture +def L(): + """An isotropic stiffness matrix.""" + return sim.L_iso([70000, 0.3], "Enu") + + +@pytest.fixture +def M(): + """An isotropic compliance matrix.""" + return sim.M_iso([70000, 0.3], "Enu") + + +@pytest.fixture +def sigma(): + return np.array([100.0, 50.0, 0.0, 25.0, 0.0, 0.0]) + + +@pytest.fixture +def epsilon(): + return np.array([0.01, -0.005, -0.005, 0.002, 0.001, 0.0]) + + +# =================================================================== +# 1. Type identity +# =================================================================== + +class TestTypeIdentity: + """simcoon.Rotation IS a scipy Rotation.""" + + def test_isinstance_scipy(self, sim_rot): + assert isinstance(sim_rot, ScipyRotation) + + def test_isinstance_simcoon(self, sim_rot): + assert isinstance(sim_rot, sim.Rotation) + + def test_from_euler_returns_simcoon(self): + r = sim.Rotation.from_euler("ZXZ", [0.5, 0.3, 0.7]) + assert isinstance(r, sim.Rotation) + + def test_from_quat_returns_simcoon(self): + r = sim.Rotation.from_quat([0, 0, 0.38268343, 0.92387953]) + assert isinstance(r, sim.Rotation) + + def test_from_matrix_returns_simcoon(self): + R = np.eye(3) + r = sim.Rotation.from_matrix(R) + assert isinstance(r, sim.Rotation) + + def test_from_rotvec_returns_simcoon(self): + r = sim.Rotation.from_rotvec([0, 0, 0.5]) + assert isinstance(r, sim.Rotation) + + def test_identity_returns_simcoon(self): + r = sim.Rotation.identity() + assert isinstance(r, sim.Rotation) + + def test_random_returns_simcoon(self): + r = sim.Rotation.random() + assert isinstance(r, sim.Rotation) + + def test_composition_preserves_type(self, sim_rot): + r2 = sim.Rotation.from_axis_angle(np.pi / 6, 1) + assert isinstance(sim_rot * r2, sim.Rotation) + + def test_inv_preserves_type(self, sim_rot): + assert isinstance(sim_rot.inv(), sim.Rotation) + + +# =================================================================== +# 2. from_scipy conversion +# =================================================================== + +class TestFromScipy: + """Upgrade a plain scipy Rotation to simcoon.Rotation.""" + + def test_from_scipy_type(self, scipy_rot): + r = sim.Rotation.from_scipy(scipy_rot) + assert isinstance(r, sim.Rotation) + + def test_from_scipy_preserves_quaternion(self, scipy_rot): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose(r.as_quat(), scipy_rot.as_quat(), atol=1e-15) + + def test_from_scipy_preserves_matrix(self, scipy_rot): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose(r.as_matrix(), scipy_rot.as_matrix(), atol=1e-14) + + def test_from_scipy_mechanics_work(self, scipy_rot, L): + r = sim.Rotation.from_scipy(scipy_rot) + L_rot = r.apply_stiffness(L) + assert L_rot.shape == (6, 6) + + +# =================================================================== +# 3. Mechanics methods +# =================================================================== + +class TestMechanics: + """Mechanics methods on simcoon.Rotation delegate to C++ correctly.""" + + def test_apply_stress_shape(self, sim_rot, sigma): + result = sim_rot.apply_stress(sigma) + assert result.shape == (6,) + + def test_apply_strain_shape(self, sim_rot, epsilon): + result = sim_rot.apply_strain(epsilon) + assert result.shape == (6,) + + def test_apply_stiffness_shape(self, sim_rot, L): + result = sim_rot.apply_stiffness(L) + assert result.shape == (6, 6) + + def test_apply_compliance_shape(self, sim_rot, M): + result = sim_rot.apply_compliance(M) + assert result.shape == (6, 6) + + def test_apply_tensor_shape(self, sim_rot): + T = np.eye(3) + result = sim_rot.apply_tensor(T) + assert result.shape == (3, 3) + + def test_apply_strain_concentration_shape(self, sim_rot): + A = np.eye(6) + result = sim_rot.apply_strain_concentration(A) + assert result.shape == (6, 6) + + def test_apply_stress_concentration_shape(self, sim_rot): + B = np.eye(6) + result = sim_rot.apply_stress_concentration(B) + assert result.shape == (6, 6) + + def test_voigt_stress_rotation_shape(self, sim_rot): + QS = sim_rot.as_voigt_stress_rotation() + assert QS.shape == (6, 6) + + def test_voigt_strain_rotation_shape(self, sim_rot): + QE = sim_rot.as_voigt_strain_rotation() + assert QE.shape == (6, 6) + + def test_identity_stiffness_unchanged(self, L): + r = sim.Rotation.identity() + np.testing.assert_allclose(r.apply_stiffness(L), L, atol=1e-10) + + def test_roundtrip_stress(self, sim_rot, sigma): + sigma_rot = sim_rot.apply_stress(sigma) + sigma_back = sim_rot.apply_stress(sigma_rot, active=False) + np.testing.assert_allclose(sigma_back, sigma, atol=1e-10) + + def test_roundtrip_stiffness(self, sim_rot, L): + L_rot = sim_rot.apply_stiffness(L) + L_back = sim_rot.apply_stiffness(L_rot, active=False) + np.testing.assert_allclose(L_back, L, atol=1e-10) + + def test_stiffness_via_voigt(self, sim_rot, L): + """apply_stiffness(L) == QS @ L @ QS.T""" + QS = sim_rot.as_voigt_stress_rotation() + L_rot_direct = sim_rot.apply_stiffness(L) + L_rot_manual = QS @ L @ QS.T + np.testing.assert_allclose(L_rot_direct, L_rot_manual, atol=1e-10) + + +# =================================================================== +# 4. Mechanics from scipy-created rotation (via from_scipy) +# =================================================================== + +class TestMechanicsFromScipy: + """Verify that a scipy Rotation upgraded via from_scipy produces + identical results to a natively-created simcoon Rotation.""" + + def test_stiffness_matches(self, sim_rot, scipy_rot, L): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose( + r.apply_stiffness(L), sim_rot.apply_stiffness(L), atol=1e-10 + ) + + def test_compliance_matches(self, sim_rot, scipy_rot, M): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose( + r.apply_compliance(M), sim_rot.apply_compliance(M), atol=1e-10 + ) + + def test_stress_matches(self, sim_rot, scipy_rot, sigma): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose( + r.apply_stress(sigma), sim_rot.apply_stress(sigma), atol=1e-10 + ) + + def test_strain_matches(self, sim_rot, scipy_rot, epsilon): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose( + r.apply_strain(epsilon), sim_rot.apply_strain(epsilon), atol=1e-10 + ) + + def test_tensor_matches(self, sim_rot, scipy_rot): + T = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]) + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose( + r.apply_tensor(T), sim_rot.apply_tensor(T), atol=1e-10 + ) + + def test_voigt_stress_rotation_matches(self, sim_rot, scipy_rot): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose( + r.as_voigt_stress_rotation(), + sim_rot.as_voigt_stress_rotation(), + atol=1e-10, + ) + + def test_voigt_strain_rotation_matches(self, sim_rot, scipy_rot): + r = sim.Rotation.from_scipy(scipy_rot) + np.testing.assert_allclose( + r.as_voigt_strain_rotation(), + sim_rot.as_voigt_strain_rotation(), + atol=1e-10, + ) + + +# =================================================================== +# 5. Utility methods +# =================================================================== + +class TestUtilities: + + def test_equals_same(self, sim_rot): + r2 = sim.Rotation.from_axis_angle(np.pi / 4, 3) + assert sim_rot.equals(r2) + + def test_equals_different(self, sim_rot): + r2 = sim.Rotation.from_axis_angle(np.pi / 3, 3) + assert not sim_rot.equals(r2) + + def test_equals_with_scipy(self, sim_rot, scipy_rot): + assert sim_rot.equals(scipy_rot) + + def test_is_identity_true(self): + r = sim.Rotation.identity() + assert r.is_identity() + + def test_is_identity_false(self, sim_rot): + assert not sim_rot.is_identity() + + def test_slerp_to_endpoints(self, sim_rot): + r_id = sim.Rotation.identity() + r0 = r_id.slerp_to(sim_rot, 0.0) + r1 = r_id.slerp_to(sim_rot, 1.0) + np.testing.assert_allclose(r0.as_quat(), r_id.as_quat(), atol=1e-12) + np.testing.assert_allclose( + np.abs(r1.as_quat() @ sim_rot.as_quat()), 1.0, atol=1e-12 + ) + + def test_slerp_to_with_scipy_target(self, scipy_rot): + r_id = sim.Rotation.identity() + r_mid = r_id.slerp_to(scipy_rot, 0.5) + assert isinstance(r_mid, sim.Rotation) + + def test_from_axis_angle_degrees(self): + r_rad = sim.Rotation.from_axis_angle(np.pi / 4, 3) + r_deg = sim.Rotation.from_axis_angle(45, 3, degrees=True) + np.testing.assert_allclose(r_rad.as_quat(), r_deg.as_quat(), atol=1e-12) + + def test_from_axis_angle_invalid_axis(self): + with pytest.raises(ValueError, match="axis must be 1, 2, or 3"): + sim.Rotation.from_axis_angle(0.5, 4) diff --git a/src/Continuum_mechanics/Functions/recovery_props.cpp b/src/Continuum_mechanics/Functions/recovery_props.cpp index c77041af4..41bf049ee 100755 --- a/src/Continuum_mechanics/Functions/recovery_props.cpp +++ b/src/Continuum_mechanics/Functions/recovery_props.cpp @@ -58,29 +58,29 @@ void check_symetries(const mat &L, std::string &umat_type, int &axis, vec &props //Check reflexions //reflexion around the z axis - symmetry plane xy Q = { {1,0,0}, {0,1,0}, {0,0,-1} }; - L_test = rotateL(L_sym, Q); + L_test = rotate_stiffness(L_sym, Q); check_sym(2,0) = norm(L_sym - L_test,2); //reflexion around the y axis - symmetry plane xz Q = { {1,0,0}, {0,-1,0}, {0,0,1} }; - L_test = rotateL(L_sym, Q); + L_test = rotate_stiffness(L_sym, Q); check_sym(1,0) = norm(L_sym - L_test,2); //reflexion around the x axis - symmetry plane yz Q = { {-1,0,0}, {0,1,0}, {0,0,1} }; - L_test = rotateL(L_sym, Q); + L_test = rotate_stiffness(L_sym, Q); check_sym(0,0) = norm(L_sym - L_test,2); //Check 90deg rotations //90 rotation around the z axis Q = { {0,1,0}, {-1,0,0}, {0,0,1} }; - L_test = rotateL(L_sym, Q); + L_test = rotate_stiffness(L_sym, Q); check_sym(2,1) = norm(L_sym - L_test,2); //90 rotation around the y axis Q = { {0,0,1}, {0,1,0}, {-1,0,0} }; - L_test = rotateL(L_sym, Q); + L_test = rotate_stiffness(L_sym, Q); check_sym(1,1) = norm(L_sym - L_test,2); //90 rotation around the x axis Q = { {1,0,0}, {0,0,1}, {0,-1,0} }; - L_test = rotateL(L_sym, Q); + L_test = rotate_stiffness(L_sym, Q); check_sym(0,1) = norm(L_sym - L_test,2); //Check equality between constants diff --git a/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp b/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp index d8d2fdc3f..3928b1964 100755 --- a/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp +++ b/src/Continuum_mechanics/Homogenization/ellipsoid_multi.cpp @@ -112,35 +112,38 @@ ellipsoid_multi::~ellipsoid_multi() {} void ellipsoid_multi::fillS_loc(const mat& Lt_m, const ellipsoid &ell) //------------------------------------- { - mat Ltm_local_geom = rotate_g2l_L(Lt_m, ell.psi_geom, ell.theta_geom, ell.phi_geom); + Rotation rot = Rotation::from_euler(ell.psi_geom, ell.theta_geom, ell.phi_geom, "zxz"); + mat Ltm_local_geom = rot.apply_stiffness(Lt_m, false); S_loc = Eshelby(Ltm_local_geom, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); } - + //------------------------------------- void ellipsoid_multi::fillP_loc(const mat& Lt_m, const ellipsoid &ell) //------------------------------------- { - mat Ltm_local_geom = rotate_g2l_L(Lt_m, ell.psi_geom, ell.theta_geom, ell.phi_geom); + Rotation rot = Rotation::from_euler(ell.psi_geom, ell.theta_geom, ell.phi_geom, "zxz"); + mat Ltm_local_geom = rot.apply_stiffness(Lt_m, false); P_loc = T_II(Ltm_local_geom, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); } - + //------------------------------------- void ellipsoid_multi::fillT(const mat& Lt_m, const mat& Lt, const ellipsoid &ell) //This method correspond to the classical Eshelby method //------------------------------------- { - mat Lt_m_local_geom = rotate_g2l_L(Lt_m, ell.psi_geom, ell.theta_geom, ell.phi_geom); + Rotation rot = Rotation::from_euler(ell.psi_geom, ell.theta_geom, ell.phi_geom, "zxz"); + mat Lt_m_local_geom = rot.apply_stiffness(Lt_m, false); S_loc = Eshelby(Lt_m_local_geom, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); - mat Lt_local_geom = rotate_g2l_L(Lt, ell.psi_geom, ell.theta_geom, ell.phi_geom); - + mat Lt_local_geom = rot.apply_stiffness(Lt, false); + try { T_loc = inv(eye(6,6) + S_loc*inv(Lt_m_local_geom)*(Lt_local_geom - Lt_m_local_geom)); } catch (const std::runtime_error &e) { cerr << "Error in inv: " << e.what() << endl; throw simcoon::exception_inv("Error in inv function inside ellipsoid_multi::fillT."); - } - T = rotate_l2g_A(T_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); + } + T = rot.apply_strain_concentration(T_loc, true); } //------------------------------------- @@ -148,18 +151,19 @@ void ellipsoid_multi::fillT_iso(const mat& Lt_m, const mat& Lt, const ellipsoid //This method corresponf to the isotropization method //------------------------------------- { + Rotation rot = Rotation::from_euler(ell.psi_geom, ell.theta_geom, ell.phi_geom, "zxz"); mat Lt_m_iso = Isotropize(Lt_m); S_loc = Eshelby(Lt_m_iso, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); - mat Lt_local_geom = rotate_g2l_L(Lt, ell.psi_geom, ell.theta_geom, ell.phi_geom); - + mat Lt_local_geom = rot.apply_stiffness(Lt, false); + try { T_loc = inv(eye(6,6) + S_loc*inv(Lt_m_iso)*(Lt_local_geom - Lt_m_iso)); } catch (const std::runtime_error &e) { cerr << "Error in inv: " << e.what() << endl; throw simcoon::exception_inv("Error in inv function inside ellipsoid_multi::fillT_is."); - } - - T = rotate_l2g_A(T_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); + } + + T = rot.apply_strain_concentration(T_loc, true); } //------------------------------------- @@ -168,20 +172,21 @@ void ellipsoid_multi::fillT_mec_in(const mat& L_m, const mat& L, const ellipsoid //the interaction tensors T for the elastic and the inelastic part. //------------------------------------- { - mat L_m_local_geom = rotate_g2l_L(L_m, ell.psi_geom, ell.theta_geom, ell.phi_geom); + Rotation rot = Rotation::from_euler(ell.psi_geom, ell.theta_geom, ell.phi_geom, "zxz"); + mat L_m_local_geom = rot.apply_stiffness(L_m, false); S_loc = Eshelby(L_m_local_geom, ell.a1, ell.a2, ell.a3, x, wx, y, wy, mp, np); - mat L_local_geom = rotate_g2l_L(L, ell.psi_geom, ell.theta_geom, ell.phi_geom); - + mat L_local_geom = rot.apply_stiffness(L, false); + try { T_loc = inv(eye(6,6) + S_loc*inv(L_m_local_geom)*(L_local_geom - L_m_local_geom)); - T_in_loc = (eye(6,6)-T_loc)*inv(L_m_local_geom - L_local_geom); + T_in_loc = (eye(6,6)-T_loc)*inv(L_m_local_geom - L_local_geom); } catch (const std::runtime_error &e) { cerr << "Error in inv: " << e.what() << endl; throw simcoon::exception_inv("Error in inv function inside ellipsoid_multi::fillT_mec_in."); } - T = rotate_l2g_A(T_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); - T_in = rotate_l2g_M(T_in_loc, ell.psi_geom, ell.theta_geom, ell.phi_geom); + T = rot.apply_strain_concentration(T_loc, true); + T_in = rot.apply_compliance(T_in_loc, true); } //---------------------------------------------------------------------- diff --git a/src/Continuum_mechanics/Micromechanics/schemes.cpp b/src/Continuum_mechanics/Micromechanics/schemes.cpp index e91c7c86d..eaab628d3 100755 --- a/src/Continuum_mechanics/Micromechanics/schemes.cpp +++ b/src/Continuum_mechanics/Micromechanics/schemes.cpp @@ -386,8 +386,9 @@ void dE_Periodic_Layer(phase_characteristics &phase, const int &nbiter) { sv_r = std::dynamic_pointer_cast(r.sptr_sv_global); lay_multi = std::dynamic_pointer_cast(r.sptr_multi); lay = std::dynamic_pointer_cast(r.sptr_shape); - Lt_loc = rotate_g2l_L(sv_r->Lt, lay->psi_geom, lay->theta_geom, lay->phi_geom); - + Rotation rot_layer = Rotation::from_euler(lay->psi_geom, lay->theta_geom, lay->phi_geom, "zxz"); + Lt_loc = rot_layer.inv().apply_stiffness(sv_r->Lt); + lay_multi->Dnn(0,0) = Lt_loc(0,0); lay_multi->Dnn(0,1) = Lt_loc(0,3); lay_multi->Dnn(0,2) = Lt_loc(0,4); @@ -397,8 +398,8 @@ void dE_Periodic_Layer(phase_characteristics &phase, const int &nbiter) { lay_multi->Dnn(2,0) = Lt_loc(4,0); lay_multi->Dnn(2,1) = Lt_loc(4,3); lay_multi->Dnn(2,2) = Lt_loc(4,4); - - mat sigma_local = rotate_g2l_stress(sv_r->sigma, lay->psi_geom, lay->theta_geom, lay->phi_geom); + + vec sigma_local = rot_layer.inv().apply_stress(sv_r->sigma); lay_multi->sigma_hat(0) = sigma_local(0); lay_multi->sigma_hat(1) = sigma_local(3); lay_multi->sigma_hat(2) = sigma_local(4); @@ -449,11 +450,12 @@ void dE_Periodic_Layer(phase_characteristics &phase, const int &nbiter) { } lay = std::dynamic_pointer_cast(r.sptr_shape); - mat dEtot_local = zeros(6); + vec dEtot_local = zeros(6); dEtot_local(0) = lay_multi->dzdx1(0); dEtot_local(3) = lay_multi->dzdx1(1); dEtot_local(4) = lay_multi->dzdx1(2); - mat dEtot_global = rotate_l2g_strain(dEtot_local, lay->psi_geom, lay->theta_geom, lay->phi_geom); + Rotation rot_l2g = Rotation::from_euler(lay->psi_geom, lay->theta_geom, lay->phi_geom, "zxz"); + vec dEtot_global = rot_l2g.apply_strain(dEtot_local); sv_r->DEtot += dEtot_global; } @@ -475,7 +477,10 @@ void Lt_Periodic_Layer(phase_characteristics &phase) { sv_r = std::dynamic_pointer_cast(r.sptr_sv_global); lay_multi = std::dynamic_pointer_cast(r.sptr_multi); lay = std::dynamic_pointer_cast(r.sptr_shape); - Lt_loc = rotate_g2l_L(sv_r->Lt, lay->psi_geom, lay->theta_geom, lay->phi_geom); + { + Rotation rot = Rotation::from_euler(lay->psi_geom, lay->theta_geom, lay->phi_geom, "zxz"); + Lt_loc = rot.inv().apply_stiffness(sv_r->Lt); + } lay_multi->Dnn(0,0) = Lt_loc(0,0); lay_multi->Dnn(0,1) = Lt_loc(0,3); @@ -558,7 +563,10 @@ void Lt_Periodic_Layer(phase_characteristics &phase) { A_loc(4,4) += lay_multi->dXn(2,2); A_loc(4,5) += lay_multi->dXt(2,2); - lay_multi->A = rotate_l2g_A(A_loc, lay->psi_geom, lay->theta_geom, lay->phi_geom); + { + Rotation rot = Rotation::from_euler(lay->psi_geom, lay->theta_geom, lay->phi_geom, "zxz"); + lay_multi->A = rot.apply_strain_concentration(A_loc); + } } } diff --git a/src/Simulation/Maths/rotation.cpp b/src/Simulation/Maths/rotation.cpp index 333e3be01..c349efed0 100755 --- a/src/Simulation/Maths/rotation.cpp +++ b/src/Simulation/Maths/rotation.cpp @@ -1,25 +1,26 @@ /* This file is part of simcoon. - + simcoon is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. - + simcoon is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - + You should have received a copy of the GNU General Public License along with simcoon. If not, see . */ ///@file rotation.cpp -///@brief rotation of a Voigt tensor -///@version 1.0 +///@brief Rotation class and Voigt tensor rotation utilities. +///@version 2.0 -#include +#include #include +#include #include #include #include @@ -30,37 +31,77 @@ using namespace arma; namespace simcoon{ -vec rotate_vec(const vec &v, const mat &DR) { - return DR*v; +// ============================================================================= +// Anonymous namespace: internal helpers +// ============================================================================= + +namespace { + +// Convert axis number (1,2,3) to axis index (0,1,2) +inline int axis_to_index(int axis) { + if (axis < 1 || axis > 3) { + throw invalid_argument("Axis must be 1 (x), 2 (y), or 3 (z)"); + } + return axis - 1; } - -vec rotate_vec(const vec &v, const double &alpha, const int &axis) { - assert(axis>0); - mat DR = fillR(alpha, axis); - return rotate_vec(v, DR); + +// Parse Euler convention string to get axis sequence +// Returns axes as 0=x, 1=y, 2=z +void parse_euler_convention(const string& conv, int& axis1, int& axis2, int& axis3) { + if (conv.length() != 3) { + throw invalid_argument("Euler convention must be 3 characters (e.g., 'zxz', 'xyz')"); + } + + auto char_to_axis = [](char c) -> int { + switch (c) { + case 'x': case 'X': return 0; + case 'y': case 'Y': return 1; + case 'z': case 'Z': return 2; + default: + throw invalid_argument("Invalid axis character in Euler convention"); + } + }; + + axis1 = char_to_axis(conv[0]); + axis2 = char_to_axis(conv[1]); + axis3 = char_to_axis(conv[2]); } -mat rotate_mat(const mat &m, const mat &DR) { - return trans(DR)*m*DR; +// Normalize a quaternion to unit length +vec::fixed<4> normalize_quat(const vec::fixed<4>& q) { + double n = norm(q); + if (n < simcoon::iota) { + return {0.0, 0.0, 0.0, 1.0}; + } + return q / n; } -mat rotate_mat(const mat &m, const double &alpha, const int &axis) { - mat DR = fillR(alpha, axis); - return rotate_mat(m, DR); -} - -mat fillR(const double &alpha, const int &axis, const bool &active) { - - double act = 0; - if (active) - act = 1.; - else - act = -1.; - +// Hamilton product of two quaternions (scalar-last convention) +vec::fixed<4> quat_multiply(const vec::fixed<4>& q1, const vec::fixed<4>& q2) { + double x1 = q1(0), y1 = q1(1), z1 = q1(2), w1 = q1(3); + double x2 = q2(0), y2 = q2(1), z2 = q2(2), w2 = q2(3); + + return { + w1*x2 + x1*w2 + y1*z2 - z1*y2, + w1*y2 - x1*z2 + y1*w2 + z1*x2, + w1*z2 + x1*y2 - y1*x2 + z1*w2, + w1*w2 - x1*x2 - y1*y2 - z1*z2 + }; +} + +// Quaternion conjugate (inverse for unit quaternion) +vec::fixed<4> quat_conjugate(const vec::fixed<4>& q) { + return {-q(0), -q(1), -q(2), q(3)}; +} + +// Internal: Generate a 3x3 rotation matrix for rotation around a single axis +mat fillR_internal(const double &alpha, const int &axis, const bool &active) { + double act = active ? 1. : -1.; + mat R = zeros(3,3); double c = cos(alpha); double s = act*sin(alpha); - + switch(axis) { case 1: { R = { {1,0,0}, {0,c,-s}, {0,s,c} }; @@ -75,16 +116,16 @@ mat fillR(const double &alpha, const int &axis, const bool &active) { break; } default: { - cout << "Please choose the axis 1,2 or 3\n"; + throw invalid_argument("Please choose the axis 1, 2 or 3"); } } return R; } -mat fillR(const double &psi, const double &theta, const double &phi, const bool &active, const std::string &conv) { - +// Internal: Generate a 3x3 rotation matrix using Euler angles +mat fillR_euler_internal(const double &psi, const double &theta, const double &phi, const bool &active, const std::string &conv) { mat R = zeros(3,3); - + if(conv == "zxz") { double c1 = cos(psi); double s1 = sin(psi); @@ -112,752 +153,823 @@ mat fillR(const double &psi, const double &theta, const double &phi, const bool R = { {c3*c2*c1-s3*s1,c3*s1+c2*c1*s3,-c1*s2}, {-c1*s3-c3*c2*s1,c3*c1-c2*s3*s1,s2*s1}, {c3*s2,s3*s2,c2} }; } else if(conv == "user") { - - mat R1 = fillR(psi, simcoon::axis_psi, active); - mat R2 = fillR(theta, simcoon::axis_theta, active); - mat R3 = fillR(phi, simcoon::axis_phi, active); - + mat R1 = fillR_internal(psi, simcoon::axis_psi, active); + mat R2 = fillR_internal(theta, simcoon::axis_theta, active); + mat R3 = fillR_internal(phi, simcoon::axis_phi, active); R = R3*R2*R1; } else { - cout << "error in Simulation/Maths/rotation.cpp: please provide a consistent convention for Euler rotation, i.e. zxz or zyz. The by-default convention has been utilized (as defined in parameter.hpp)" << endl; - mat R1 = fillR(psi, simcoon::axis_psi, active); - mat R2 = fillR(theta, simcoon::axis_theta, active); - mat R3 = fillR(phi, simcoon::axis_phi, active); - - R = R3*R2*R1; + throw std::invalid_argument("error in Simulation/Maths/rotation.cpp: invalid Euler convention '" + conv + "'. Supported conventions: zxz, zyz, xyz, xzy, yxz, yzx, zxy, zyx, xyx, xzx, yxy, yzy, user"); } - + return R; } - -mat fillQS(const double &alpha, const int &axis, const bool &active) { - double act = 0; - if (active) { - act = 1.; +// Internal: Generate 6x6 stress Voigt rotation matrix from angle/axis +mat fill_voigt_stress(const double &alpha, const int &axis, const bool &active) { + double act = active ? 1. : -1.; + + mat voigt_stress = zeros(6,6); + double c = cos(alpha); + double s = act*sin(alpha); + + switch(axis) { + case 1: { + voigt_stress(0,0) = 1.; + voigt_stress(1,1) = c*c; + voigt_stress(1,2) = s*s; + voigt_stress(1,5) = -2.*s*c; + voigt_stress(2,1) = s*s; + voigt_stress(2,2) = c*c; + voigt_stress(2,5) = 2.*s*c; + voigt_stress(3,3) = c; + voigt_stress(3,4) = -s; + voigt_stress(4,3) = s; + voigt_stress(4,4) = c; + voigt_stress(5,1) = s*c; + voigt_stress(5,2) = -s*c; + voigt_stress(5,5) = c*c-s*s; + break; + } + case 2: { + voigt_stress(0,0) = c*c; + voigt_stress(0,2) = s*s; + voigt_stress(0,4) = 2*c*s; + voigt_stress(1,1) = 1.; + voigt_stress(2,0) = s*s; + voigt_stress(2,2) = c*c; + voigt_stress(2,4) = -2*c*s; + voigt_stress(3,3) = c; + voigt_stress(3,5) = s; + voigt_stress(4,0) = -c*s; + voigt_stress(4,2) = c*s; + voigt_stress(4,4) = c*c-s*s; + voigt_stress(5,3) = -s; + voigt_stress(5,5) = c; + break; + } + case 3: { + voigt_stress(0,0) = c*c; + voigt_stress(0,1) = s*s; + voigt_stress(0,3) = -2*s*c; + voigt_stress(1,0) = s*s; + voigt_stress(1,1) = c*c; + voigt_stress(1,3) = 2*s*c; + voigt_stress(2,2) = 1.; + voigt_stress(3,0) = s*c; + voigt_stress(3,1) = -s*c; + voigt_stress(3,3) = c*c-s*s; + voigt_stress(4,4) = c; + voigt_stress(4,5) = -s; + voigt_stress(5,4) = s; + voigt_stress(5,5) = c; + break; + } + default: { + throw invalid_argument("Please choose the axis 1, 2 or 3"); + } } - else { - act = -1.; - } - - mat QS = zeros(6,6); - double c = cos(alpha); - double s = act*sin(alpha); - - switch(axis) { - case 1: { - QS(0,0) = 1.; - QS(0,1) = 0.; - QS(0,2) = 0.; - QS(0,3) = 0.; - QS(0,4) = 0.; - QS(0,5) = 0.; - QS(1,0) = 0.; - QS(1,1) = c*c; - QS(1,2) = s*s; - QS(1,3) = 0.; - QS(1,4) = 0.; - QS(1,5) = -2.*s*c; - QS(2,0) = 0.; - QS(2,1) = s*s; - QS(2,2) = c*c; - QS(2,3) = 0.; - QS(2,4) = 0.; - QS(2,5) = 2.*s*c; - QS(3,0) = 0.; - QS(3,1) = 0.; - QS(3,2) = 0.; - QS(3,3) = c; - QS(3,4) = -s; - QS(3,5) = 0.; - QS(4,0) = 0.; - QS(4,1) = 0.; - QS(4,2) = 0.; - QS(4,3) = s; - QS(4,4) = c; - QS(4,5) = 0.; - QS(5,0) = 0.; - QS(5,1) = s*c; - QS(5,2) = -s*c; - QS(5,3) = 0.; - QS(5,4) = 0.; - QS(5,5) = c*c-s*s; - - break; - } - case 2: { - QS(0,0) = c*c; - QS(0,1) = 0.; - QS(0,2) = s*s; - QS(0,3) = 0.; - QS(0,4) = 2*c*s; - QS(0,5) = 0.; - QS(1,0) = 0.; - QS(1,1) = 1.; - QS(1,2) = 0.; - QS(1,3) = 0.; - QS(1,4) = 0.; - QS(1,5) = 0.; - QS(2,0) = s*s; - QS(2,1) = 0.; - QS(2,2) = c*c; - QS(2,3) = 0.; - QS(2,4) = -2*c*s; - QS(2,5) = 0.; - QS(3,0) = 0.; - QS(3,1) = 0.; - QS(3,2) = 0.; - QS(3,3) = c; - QS(3,4) = 0.; - QS(3,5) = s; - QS(4,0) = -c*s; - QS(4,1) = 0.; - QS(4,2) = c*s; - QS(4,3) = 0.; - QS(4,4) = c*c-s*s; - QS(4,5) = 0.; - QS(5,0) = 0.; - QS(5,1) = 0.; - QS(5,2) = 0.; - QS(5,3) = -s; - QS(5,4) = 0.; - QS(5,5) = c; - - break; - - } - case 3: { - QS(0,0) = c*c; - QS(0,1) = s*s; - QS(0,2) = 0.; - QS(0,3) = -2*s*c; - QS(0,4) = 0.; - QS(0,5) = 0.; - QS(1,0) = s*s; - QS(1,1) = c*c; - QS(1,2) = 0.; - QS(1,3) = 2*s*c; - QS(1,4) = 0.; - QS(1,5) = 0.; - QS(2,0) = 0.; - QS(2,1) = 0.; - QS(2,2) = 1.; - QS(2,3) = 0.; - QS(2,4) = 0.; - QS(2,5) = 0.; - QS(3,0) = s*c; - QS(3,1) = -s*c; - QS(3,2) = 0.; - QS(3,3) = c*c-s*s; - QS(3,4) = 0.; - QS(3,5) = 0.; - QS(4,0) = 0.; - QS(4,1) = 0.; - QS(4,2) = 0.; - QS(4,3) = 0.; - QS(4,4) = c; - QS(4,5) = -s; - QS(5,0) = 0.; - QS(5,1) = 0.; - QS(5,2) = 0.; - QS(5,3) = 0.; - QS(5,4) = s; - QS(5,5) = c; - - break; - } - default: { - cout << "Please choose the axis 1,2 or 3\n"; - } - } - return QS; -} - -mat fillQS(const mat &DR, const bool &active) { - - double a = 0.; - double d = 0.; - double g = 0.; - double b = 0.; - double e = 0.; - double h = 0.; - double c = 0.; - double f = 0.; - double i = 0.; - + return voigt_stress; +} + +// Internal: Generate 6x6 stress Voigt rotation matrix from rotation matrix +mat fill_voigt_stress(const mat &DR, const bool &active) { + double a, b, c, d, e, f, g, h, i; + if(active) { - a = DR(0,0); - b = DR(0,1); - c = DR(0,2); - d = DR(1,0); - e = DR(1,1); - f = DR(1,2); - g = DR(2,0); - h = DR(2,1); - i = DR(2,2); - } - else{ - a = DR(0,0); - d = DR(0,1); - g = DR(0,2); - b = DR(1,0); - e = DR(1,1); - h = DR(1,2); - c = DR(2,0); - f = DR(2,1); - i = DR(2,2); - } - - mat QS= zeros(6,6); - QS(0,0) = a*a; - QS(0,1) = b*b; - QS(0,2) = c*c; - QS(0,3) = 2.*a*b; - QS(0,4) = 2.*a*c; - QS(0,5) = 2.*b*c; - QS(1,0) = d*d; - QS(1,1) = e*e; - QS(1,2) = f*f; - QS(1,3) = 2.*d*e; - QS(1,4) = 2.*d*f; - QS(1,5) = 2.*e*f; - QS(2,0) = g*g; - QS(2,1) = h*h; - QS(2,2) = i*i; - QS(2,3) = 2.*g*h; - QS(2,4) = 2.*g*i; - QS(2,5) = 2.*h*i; - QS(3,0) = a*d; - QS(3,1) = b*e; - QS(3,2) = c*f; - QS(3,3) = d*b+a*e; - QS(3,4) = d*c+a*f; - QS(3,5) = e*c+b*f; - QS(4,0) = a*g; - QS(4,1) = b*h; - QS(4,2) = c*i; - QS(4,3) = g*b+a*h; - QS(4,4) = g*c+a*i; - QS(4,5) = h*c+b*i; - QS(5,0) = d*g; - QS(5,1) = e*h; - QS(5,2) = f*i; - QS(5,3) = g*e+d*h; - QS(5,4) = g*f+d*i; - QS(5,5) = h*f+e*i; - - return QS; -} - -mat fillQE(const double &alpha, const int &axis, const bool &active) { - - double act = 0; - if (active) { - act = 1.; - } - else { - act = -1.; - } - - mat QE = zeros(6,6); + a = DR(0,0); b = DR(0,1); c = DR(0,2); + d = DR(1,0); e = DR(1,1); f = DR(1,2); + g = DR(2,0); h = DR(2,1); i = DR(2,2); + } else { + a = DR(0,0); d = DR(0,1); g = DR(0,2); + b = DR(1,0); e = DR(1,1); h = DR(1,2); + c = DR(2,0); f = DR(2,1); i = DR(2,2); + } + + mat voigt_stress = zeros(6,6); + voigt_stress(0,0) = a*a; voigt_stress(0,1) = b*b; voigt_stress(0,2) = c*c; + voigt_stress(0,3) = 2.*a*b; voigt_stress(0,4) = 2.*a*c; voigt_stress(0,5) = 2.*b*c; + voigt_stress(1,0) = d*d; voigt_stress(1,1) = e*e; voigt_stress(1,2) = f*f; + voigt_stress(1,3) = 2.*d*e; voigt_stress(1,4) = 2.*d*f; voigt_stress(1,5) = 2.*e*f; + voigt_stress(2,0) = g*g; voigt_stress(2,1) = h*h; voigt_stress(2,2) = i*i; + voigt_stress(2,3) = 2.*g*h; voigt_stress(2,4) = 2.*g*i; voigt_stress(2,5) = 2.*h*i; + voigt_stress(3,0) = a*d; voigt_stress(3,1) = b*e; voigt_stress(3,2) = c*f; + voigt_stress(3,3) = d*b+a*e; voigt_stress(3,4) = d*c+a*f; voigt_stress(3,5) = e*c+b*f; + voigt_stress(4,0) = a*g; voigt_stress(4,1) = b*h; voigt_stress(4,2) = c*i; + voigt_stress(4,3) = g*b+a*h; voigt_stress(4,4) = g*c+a*i; voigt_stress(4,5) = h*c+b*i; + voigt_stress(5,0) = d*g; voigt_stress(5,1) = e*h; voigt_stress(5,2) = f*i; + voigt_stress(5,3) = g*e+d*h; voigt_stress(5,4) = g*f+d*i; voigt_stress(5,5) = h*f+e*i; + + return voigt_stress; +} + +// Internal: Generate 6x6 strain Voigt rotation matrix from angle/axis +mat fill_voigt_strain(const double &alpha, const int &axis, const bool &active) { + double act = active ? 1. : -1.; + + mat voigt_strain = zeros(6,6); double c = cos(alpha); double s = act*sin(alpha); - - switch(axis) { - case 1: { - QE(0,0) = 1.; - QE(0,1) = 0.; - QE(0,2) = 0.; - QE(0,3) = 0.; - QE(0,4) = 0.; - QE(0,5) = 0.; - QE(1,0) = 0.; - QE(1,1) = c*c; - QE(1,2) = s*s; - QE(1,3) = 0.; - QE(1,4) = 0.; - QE(1,5) = -s*c; - QE(2,0) = 0.; - QE(2,1) = s*s; - QE(2,2) = c*c; - QE(2,3) = 0.; - QE(2,4) = 0.; - QE(2,5) = s*c; - QE(3,0) = 0.; - QE(3,1) = 0.; - QE(3,2) = 0.; - QE(3,3) = c; - QE(3,4) = -s; - QE(3,5) = 0.; - QE(4,0) = 0.; - QE(4,1) = 0.; - QE(4,2) = 0.; - QE(4,3) = s; - QE(4,4) = c; - QE(4,5) = 0.; - QE(5,0) = 0.; - QE(5,1) = 2.*s*c; - QE(5,2) = -2.*s*c; - QE(5,3) = 0.; - QE(5,4) = 0.; - QE(5,5) = c*c-s*s; - break; - } - case 2: { - QE(0,0) = c*c; - QE(0,1) = 0.; - QE(0,2) = s*s; - QE(0,3) = 0.; - QE(0,4) = c*s; - QE(0,5) = 0.; - QE(1,0) = 0.; - QE(1,1) = 1.; - QE(1,2) = 0.; - QE(1,3) = 0.; - QE(1,4) = 0.; - QE(1,5) = 0.; - QE(2,0) = s*s; - QE(2,1) = 0.; - QE(2,2) = c*c; - QE(2,3) = 0.; - QE(2,4) = -c*s; - QE(2,5) = 0.; - QE(3,0) = 0.; - QE(3,1) = 0.; - QE(3,2) = 0.; - QE(3,3) = c; - QE(3,4) = 0.; - QE(3,5) = s; - QE(4,0) = -2.*c*s; - QE(4,1) = 0.; - QE(4,2) = 2.*c*s; - QE(4,3) = 0.; - QE(4,4) = c*c-s*s; - QE(4,5) = 0.; - QE(5,0) = 0.; - QE(5,1) = 0.; - QE(5,2) = 0.; - QE(5,3) = -s; - QE(5,4) = 0.; - QE(5,5) = c; - break; - } - case 3: { - QE(0,0) = c*c; - QE(0,1) = s*s; - QE(0,2) = 0.; - QE(0,3) = -s*c; - QE(0,4) = 0.; - QE(0,5) = 0.; - QE(1,0) = s*s; - QE(1,1) = c*c; - QE(1,2) = 0.; - QE(1,3) = s*c; - QE(1,4) = 0.; - QE(1,5) = 0.; - QE(2,0) = 0.; - QE(2,1) = 0.; - QE(2,2) = 1.; - QE(2,3) = 0.; - QE(2,4) = 0.; - QE(2,5) = 0.; - QE(3,0) = 2.*s*c; - QE(3,1) = -2.*s*c; - QE(3,2) = 0.; - QE(3,3) = c*c-s*s; - QE(3,4) = 0.; - QE(3,5) = 0.; - QE(4,0) = 0.; - QE(4,1) = 0.; - QE(4,2) = 0.; - QE(4,3) = 0.; - QE(4,4) = c; - QE(4,5) = -s; - QE(5,0) = 0.; - QE(5,1) = 0.; - QE(5,2) = 0.; - QE(5,3) = 0.; - QE(5,4) = s; - QE(5,5) = c; - break; - } - default: { - cout << "Please choose the axis 1,2 or 3\n"; - } - } - return QE; -} - -mat fillQE(const mat &DR, const bool &active) { - - double a = 0.; - double d = 0.; - double g = 0.; - double b = 0.; - double e = 0.; - double h = 0.; - double c = 0.; - double f = 0.; - double i = 0.; - + + switch(axis) { + case 1: { + voigt_strain(0,0) = 1.; + voigt_strain(1,1) = c*c; + voigt_strain(1,2) = s*s; + voigt_strain(1,5) = -s*c; + voigt_strain(2,1) = s*s; + voigt_strain(2,2) = c*c; + voigt_strain(2,5) = s*c; + voigt_strain(3,3) = c; + voigt_strain(3,4) = -s; + voigt_strain(4,3) = s; + voigt_strain(4,4) = c; + voigt_strain(5,1) = 2.*s*c; + voigt_strain(5,2) = -2.*s*c; + voigt_strain(5,5) = c*c-s*s; + break; + } + case 2: { + voigt_strain(0,0) = c*c; + voigt_strain(0,2) = s*s; + voigt_strain(0,4) = c*s; + voigt_strain(1,1) = 1.; + voigt_strain(2,0) = s*s; + voigt_strain(2,2) = c*c; + voigt_strain(2,4) = -c*s; + voigt_strain(3,3) = c; + voigt_strain(3,5) = s; + voigt_strain(4,0) = -2.*c*s; + voigt_strain(4,2) = 2.*c*s; + voigt_strain(4,4) = c*c-s*s; + voigt_strain(5,3) = -s; + voigt_strain(5,5) = c; + break; + } + case 3: { + voigt_strain(0,0) = c*c; + voigt_strain(0,1) = s*s; + voigt_strain(0,3) = -s*c; + voigt_strain(1,0) = s*s; + voigt_strain(1,1) = c*c; + voigt_strain(1,3) = s*c; + voigt_strain(2,2) = 1.; + voigt_strain(3,0) = 2.*s*c; + voigt_strain(3,1) = -2.*s*c; + voigt_strain(3,3) = c*c-s*s; + voigt_strain(4,4) = c; + voigt_strain(4,5) = -s; + voigt_strain(5,4) = s; + voigt_strain(5,5) = c; + break; + } + default: { + throw invalid_argument("Please choose the axis 1, 2 or 3"); + } + } + return voigt_strain; +} + +// Internal: Generate 6x6 strain Voigt rotation matrix from rotation matrix +mat fill_voigt_strain(const mat &DR, const bool &active) { + double a, b, c, d, e, f, g, h, i; + if(active) { - a = DR(0,0); - b = DR(0,1); - c = DR(0,2); - d = DR(1,0); - e = DR(1,1); - f = DR(1,2); - g = DR(2,0); - h = DR(2,1); - i = DR(2,2); - } - else{ - a = DR(0,0); - d = DR(0,1); - g = DR(0,2); - b = DR(1,0); - e = DR(1,1); - h = DR(1,2); - c = DR(2,0); - f = DR(2,1); - i = DR(2,2); - } - - mat QE= zeros(6,6); - QE(0,0) = a*a; - QE(0,1) = b*b; - QE(0,2) = c*c; - QE(0,3) = a*b; - QE(0,4) = a*c; - QE(0,5) = b*c; - QE(1,0) = d*d; - QE(1,1) = e*e; - QE(1,2) = f*f; - QE(1,3) = d*e; - QE(1,4) = d*f; - QE(1,5) = e*f; - QE(2,0) = g*g; - QE(2,1) = h*h; - QE(2,2) = i*i; - QE(2,3) = g*h; - QE(2,4) = g*i; - QE(2,5) = h*i; - QE(3,0) = 2.*a*d; - QE(3,1) = 2.*b*e; - QE(3,2) = 2.*c*f; - QE(3,3) = d*b+a*e; - QE(3,4) = d*c+a*f; - QE(3,5) = e*c+b*f; - QE(4,0) = 2.*a*g; - QE(4,1) = 2.*b*h; - QE(4,2) = 2.*c*i; - QE(4,3) = g*b+a*h; - QE(4,4) = g*c+a*i; - QE(4,5) = h*c+b*i; - QE(5,0) = 2.*d*g; - QE(5,1) = 2.*e*h; - QE(5,2) = 2.*f*i; - QE(5,3) = g*e+d*h; - QE(5,4) = g*f+d*i; - QE(5,5) = h*f+e*i; - - return QE; -} - - -//To rotate a stiffness matrix (6,6) -mat rotateL(const mat &L, const double &alpha, const int &axis, const bool &active) { - - mat QS = fillQS(alpha, axis, active); - return QS*(L*trans(QS)); -} - -mat rotateL(const mat &L, const mat &DR, const bool &active) { - - mat QS = fillQS(DR, active); - return QS*(L*trans(QS)); -} - - -//To rotate a compliance matrix (6,6) -mat rotateM(const mat &M, const double &alpha, const int &axis, const bool &active) { - - mat QE = fillQE(alpha, axis, active); - return QE*(M*trans(QE)); -} - -mat rotateM(const mat &M, const mat &DR, const bool &active) { - - mat QE = fillQE(DR, active); - return QE*(M*trans(QE)); -} - - -//To rotate an interaction matrix of type A (strain) (6,6) -mat rotateA(const mat &A, const double &alpha, const int &axis, const bool &active) { - - mat QE = fillQE(alpha, axis, active); - mat QS = fillQS(alpha, axis, active); - - return QE*(A*trans(QS)); -} - -mat rotateA(const mat &A, const mat &DR, const bool &active) { - - mat QE = fillQE(DR, active); - mat QS = fillQS(DR, active); - - return QE*(A*trans(QS)); -} - -//To rotate an interaction matrix of type B (stress) (6,6) -mat rotateB(const mat &B, const double &alpha, const int &axis, const bool &active) { - - mat QE = fillQE(alpha, axis, active); - mat QS = fillQS(alpha, axis, active); - - return QS*(B*trans(QE)); -} - -mat rotateB(const mat &B, const mat &DR, const bool &active) { - - mat QE = fillQE(DR, active); - mat QS = fillQS(DR, active); - - return QS*(B*trans(QE)); -} - -//To rotate a stress vector (6) -vec rotate_stress(const vec &V, const double &alpha, const int &axis, const bool &active) { + a = DR(0,0); b = DR(0,1); c = DR(0,2); + d = DR(1,0); e = DR(1,1); f = DR(1,2); + g = DR(2,0); h = DR(2,1); i = DR(2,2); + } else { + a = DR(0,0); d = DR(0,1); g = DR(0,2); + b = DR(1,0); e = DR(1,1); h = DR(1,2); + c = DR(2,0); f = DR(2,1); i = DR(2,2); + } - mat QS = fillQS(alpha, axis, active); - return QS*V; + mat voigt_strain = zeros(6,6); + voigt_strain(0,0) = a*a; voigt_strain(0,1) = b*b; voigt_strain(0,2) = c*c; + voigt_strain(0,3) = a*b; voigt_strain(0,4) = a*c; voigt_strain(0,5) = b*c; + voigt_strain(1,0) = d*d; voigt_strain(1,1) = e*e; voigt_strain(1,2) = f*f; + voigt_strain(1,3) = d*e; voigt_strain(1,4) = d*f; voigt_strain(1,5) = e*f; + voigt_strain(2,0) = g*g; voigt_strain(2,1) = h*h; voigt_strain(2,2) = i*i; + voigt_strain(2,3) = g*h; voigt_strain(2,4) = g*i; voigt_strain(2,5) = h*i; + voigt_strain(3,0) = 2.*a*d; voigt_strain(3,1) = 2.*b*e; voigt_strain(3,2) = 2.*c*f; + voigt_strain(3,3) = d*b+a*e; voigt_strain(3,4) = d*c+a*f; voigt_strain(3,5) = e*c+b*f; + voigt_strain(4,0) = 2.*a*g; voigt_strain(4,1) = 2.*b*h; voigt_strain(4,2) = 2.*c*i; + voigt_strain(4,3) = g*b+a*h; voigt_strain(4,4) = g*c+a*i; voigt_strain(4,5) = h*c+b*i; + voigt_strain(5,0) = 2.*d*g; voigt_strain(5,1) = 2.*e*h; voigt_strain(5,2) = 2.*f*i; + voigt_strain(5,3) = g*e+d*h; voigt_strain(5,4) = g*f+d*i; voigt_strain(5,5) = h*f+e*i; + + return voigt_strain; +} + +} // anonymous namespace + +// ============================================================================= +// Rotation Class: Constructors +// ============================================================================= + +Rotation::Rotation() : _quat({0.0, 0.0, 0.0, 1.0}) { +} + +Rotation::Rotation(const vec::fixed<4>& quat) : _quat(quat) { +} + +// ============================================================================= +// Rotation Class: Factory Methods +// ============================================================================= + +Rotation Rotation::identity() { + return Rotation(); +} + +Rotation Rotation::from_quat(const vec::fixed<4>& quat) { + return Rotation(normalize_quat(quat)); +} + +Rotation Rotation::from_quat(const vec& quat) { + if (quat.n_elem != 4) { + throw invalid_argument("Quaternion must have 4 elements"); + } + vec::fixed<4> q_fixed = {quat(0), quat(1), quat(2), quat(3)}; + return from_quat(q_fixed); +} + +Rotation Rotation::from_matrix(const mat::fixed<3,3>& R) { + double trace = R(0,0) + R(1,1) + R(2,2); + vec::fixed<4> q; + + if (trace > 0) { + double s = 0.5 / sqrt(trace + 1.0); + q(3) = 0.25 / s; + q(0) = (R(2,1) - R(1,2)) * s; + q(1) = (R(0,2) - R(2,0)) * s; + q(2) = (R(1,0) - R(0,1)) * s; + } else if (R(0,0) > R(1,1) && R(0,0) > R(2,2)) { + double s = 2.0 * sqrt(1.0 + R(0,0) - R(1,1) - R(2,2)); + q(3) = (R(2,1) - R(1,2)) / s; + q(0) = 0.25 * s; + q(1) = (R(0,1) + R(1,0)) / s; + q(2) = (R(0,2) + R(2,0)) / s; + } else if (R(1,1) > R(2,2)) { + double s = 2.0 * sqrt(1.0 + R(1,1) - R(0,0) - R(2,2)); + q(3) = (R(0,2) - R(2,0)) / s; + q(0) = (R(0,1) + R(1,0)) / s; + q(1) = 0.25 * s; + q(2) = (R(1,2) + R(2,1)) / s; + } else { + double s = 2.0 * sqrt(1.0 + R(2,2) - R(0,0) - R(1,1)); + q(3) = (R(1,0) - R(0,1)) / s; + q(0) = (R(0,2) + R(2,0)) / s; + q(1) = (R(1,2) + R(2,1)) / s; + q(2) = 0.25 * s; + } + + return Rotation(normalize_quat(q)); +} + +Rotation Rotation::from_matrix(const mat& R) { + if (R.n_rows != 3 || R.n_cols != 3) { + throw invalid_argument("Rotation matrix must be 3x3"); + } + mat::fixed<3,3> R_fixed; + R_fixed = R; + return from_matrix(R_fixed); +} + +Rotation Rotation::from_euler(double psi, double theta, double phi, + const string& conv, bool intrinsic, bool degrees) { + if (degrees) { + psi *= simcoon::pi / 180.0; + theta *= simcoon::pi / 180.0; + phi *= simcoon::pi / 180.0; + } + + if (conv == "user") { + mat R = fillR_euler_internal(psi, theta, phi, true, "user"); + mat::fixed<3,3> R_fixed; + R_fixed = R; + return from_matrix(R_fixed); + } + + int axis1, axis2, axis3; + parse_euler_convention(conv, axis1, axis2, axis3); + + auto axis_quat = [](double angle, int axis) -> vec::fixed<4> { + double half_angle = angle / 2.0; + double s = sin(half_angle); + double c = cos(half_angle); + vec::fixed<4> q = {0, 0, 0, c}; + q(axis) = s; + return q; + }; + + vec::fixed<4> q1 = axis_quat(psi, axis1); + vec::fixed<4> q2 = axis_quat(theta, axis2); + vec::fixed<4> q3 = axis_quat(phi, axis3); + + vec::fixed<4> result; + if (intrinsic) { + result = quat_multiply(quat_multiply(q3, q2), q1); + } else { + result = quat_multiply(quat_multiply(q1, q2), q3); + } + + return Rotation(normalize_quat(result)); +} + +Rotation Rotation::from_euler(const vec::fixed<3>& angles, const string& conv, + bool intrinsic, bool degrees) { + return from_euler(angles(0), angles(1), angles(2), conv, intrinsic, degrees); +} + +Rotation Rotation::from_rotvec(const vec::fixed<3>& rotvec, bool degrees) { + double angle = norm(rotvec); + + if (degrees) { + angle *= simcoon::pi / 180.0; + } + + if (angle < simcoon::iota) { + return identity(); + } + + vec::fixed<3> axis = rotvec / (degrees ? norm(rotvec) * 180.0 / simcoon::pi : norm(rotvec)); + double half_angle = angle / 2.0; + double s = sin(half_angle); + double c = cos(half_angle); + + return Rotation({axis(0)*s, axis(1)*s, axis(2)*s, c}); +} + +Rotation Rotation::from_rotvec(const vec& rotvec, bool degrees) { + if (rotvec.n_elem != 3) { + throw invalid_argument("Rotation vector must have 3 elements"); + } + vec::fixed<3> rv_fixed = {rotvec(0), rotvec(1), rotvec(2)}; + return from_rotvec(rv_fixed, degrees); +} + +Rotation Rotation::from_axis_angle(double angle, int axis, bool degrees) { + if (degrees) { + angle *= simcoon::pi / 180.0; + } + + int ax_idx = axis_to_index(axis); + double half_angle = angle / 2.0; + double s = sin(half_angle); + double c = cos(half_angle); + + vec::fixed<4> q = {0, 0, 0, c}; + q(ax_idx) = s; + + return Rotation(q); +} + +Rotation Rotation::random() { + vec u = randu(3); + + double sqrt1_u0 = sqrt(1.0 - u(0)); + double sqrt_u0 = sqrt(u(0)); + double two_pi_u1 = 2.0 * simcoon::pi * u(1); + double two_pi_u2 = 2.0 * simcoon::pi * u(2); + + vec::fixed<4> q = { + sqrt1_u0 * sin(two_pi_u1), + sqrt1_u0 * cos(two_pi_u1), + sqrt_u0 * sin(two_pi_u2), + sqrt_u0 * cos(two_pi_u2) + }; + + return Rotation(q); +} + +// ============================================================================= +// Rotation Class: Conversion Methods +// ============================================================================= + +mat::fixed<3,3> Rotation::as_matrix() const { + double qx = _quat(0), qy = _quat(1), qz = _quat(2), qw = _quat(3); + + double qx2 = qx * qx, qy2 = qy * qy, qz2 = qz * qz; + double qxqy = qx * qy, qxqz = qx * qz, qxqw = qx * qw; + double qyqz = qy * qz, qyqw = qy * qw, qzqw = qz * qw; + + mat::fixed<3,3> R; + R(0,0) = 1.0 - 2.0*(qy2 + qz2); + R(0,1) = 2.0*(qxqy - qzqw); + R(0,2) = 2.0*(qxqz + qyqw); + R(1,0) = 2.0*(qxqy + qzqw); + R(1,1) = 1.0 - 2.0*(qx2 + qz2); + R(1,2) = 2.0*(qyqz - qxqw); + R(2,0) = 2.0*(qxqz - qyqw); + R(2,1) = 2.0*(qyqz + qxqw); + R(2,2) = 1.0 - 2.0*(qx2 + qy2); + + return R; +} + +vec::fixed<3> Rotation::as_euler(const string& conv, bool intrinsic, bool degrees) const { + mat::fixed<3,3> R = as_matrix(); + + int axis1, axis2, axis3; + + if (conv == "user") { + axis1 = axis_psi - 1; + axis2 = axis_theta - 1; + axis3 = axis_phi - 1; + } else { + parse_euler_convention(conv, axis1, axis2, axis3); + } + + mat::fixed<3,3> R_work = intrinsic ? R : R.t(); + if (!intrinsic) { + swap(axis1, axis3); + } + + vec::fixed<3> angles; + + bool proper_euler = (axis1 == axis3); + + int i = axis1; + int j = axis2; + int k = proper_euler ? (3 - i - j) : axis3; + + double sign = ((j - i + 3) % 3 == 1) ? 1.0 : -1.0; + + if (proper_euler) { + double sy = sqrt(R_work(i,j)*R_work(i,j) + R_work(i,k)*R_work(i,k)); + + if (sy > simcoon::iota) { + angles(0) = atan2(R_work(i,j), sign*R_work(i,k)); + angles(1) = atan2(sy, R_work(i,i)); + angles(2) = atan2(R_work(j,i), -sign*R_work(k,i)); + } else { + angles(0) = atan2(-sign*R_work(j,k), R_work(j,j)); + angles(1) = atan2(sy, R_work(i,i)); + angles(2) = 0.0; + } + } else { + double cy = sqrt(R_work(i,i)*R_work(i,i) + R_work(j,i)*R_work(j,i)); + + if (cy > simcoon::iota) { + angles(0) = atan2(sign*R_work(k,j), R_work(k,k)); + angles(1) = atan2(-sign*R_work(k,i), cy); + angles(2) = atan2(sign*R_work(j,i), R_work(i,i)); + } else { + angles(0) = atan2(-sign*R_work(j,k), R_work(j,j)); + angles(1) = atan2(-sign*R_work(k,i), cy); + angles(2) = 0.0; + } + } + + if (!intrinsic) { + swap(angles(0), angles(2)); + } + + if (degrees) { + angles *= 180.0 / simcoon::pi; + } + + return angles; +} + +vec::fixed<3> Rotation::as_rotvec(bool degrees) const { + double qw = _quat(3); + vec::fixed<3> qv = {_quat(0), _quat(1), _quat(2)}; + + double sin_half = norm(qv); + + if (sin_half < simcoon::iota) { + vec::fixed<3> result = {0.0, 0.0, 0.0}; + return result; + } + + double angle = 2.0 * atan2(sin_half, qw); + + if (angle > simcoon::pi) { + angle -= 2.0 * simcoon::pi; + } else if (angle < -simcoon::pi) { + angle += 2.0 * simcoon::pi; + } + + vec::fixed<3> axis = qv / sin_half; + + if (degrees) { + angle *= 180.0 / simcoon::pi; + } + + return axis * angle; +} + +mat::fixed<6,6> Rotation::as_voigt_stress_rotation(bool active) const { + mat R = as_matrix(); + mat vs_dyn = fill_voigt_stress(R, active); + mat::fixed<6,6> vs; + vs = vs_dyn; + return vs; +} + +mat::fixed<6,6> Rotation::as_voigt_strain_rotation(bool active) const { + mat R = as_matrix(); + mat ve_dyn = fill_voigt_strain(R, active); + mat::fixed<6,6> ve; + ve = ve_dyn; + return ve; +} + +// ============================================================================= +// Rotation Class: Apply Methods (3D objects) +// ============================================================================= + +vec::fixed<3> Rotation::apply(const vec::fixed<3>& v, bool inverse) const { + double qx = _quat(0), qy = _quat(1), qz = _quat(2), qw = _quat(3); + + if (inverse) { + qx = -qx; qy = -qy; qz = -qz; + } + + vec::fixed<3> qv = {qx, qy, qz}; + vec::fixed<3> uv = cross(qv, v); + vec::fixed<3> uuv = cross(qv, uv); + + return v + 2.0 * (qw * uv + uuv); +} + +vec Rotation::apply(const vec& v, bool inverse) const { + if (v.n_elem != 3) { + throw invalid_argument("Vector must have 3 elements"); + } + vec::fixed<3> v_fixed = {v(0), v(1), v(2)}; + vec::fixed<3> result = apply(v_fixed, inverse); + return vec(result); +} + +mat::fixed<3,3> Rotation::apply_tensor(const mat::fixed<3,3>& m, bool inverse) const { + mat::fixed<3,3> R = as_matrix(); + if (inverse) { + return R.t() * m * R; + } else { + return R * m * R.t(); + } +} + +mat Rotation::apply_tensor(const mat& m, bool inverse) const { + if (m.n_rows != 3 || m.n_cols != 3) { + throw invalid_argument("Tensor must be 3x3"); + } + mat::fixed<3,3> m_fixed; + m_fixed = m; + mat::fixed<3,3> result = apply_tensor(m_fixed, inverse); + return mat(result); +} + +// ============================================================================= +// Rotation Class: Apply Methods (Voigt notation) +// ============================================================================= + +vec::fixed<6> Rotation::apply_stress(const vec::fixed<6>& sigma, bool active) const { + mat::fixed<6,6> vs = as_voigt_stress_rotation(active); + vec::fixed<6> result; + result = vs * sigma; + return result; +} + +vec Rotation::apply_stress(const vec& sigma, bool active) const { + if (sigma.n_elem != 6) { + throw invalid_argument("Stress vector must have 6 elements"); + } + mat R = as_matrix(); + mat vs = fill_voigt_stress(R, active); + return vs * sigma; +} + +vec::fixed<6> Rotation::apply_strain(const vec::fixed<6>& epsilon, bool active) const { + mat::fixed<6,6> ve = as_voigt_strain_rotation(active); + vec::fixed<6> result; + result = ve * epsilon; + return result; +} + +vec Rotation::apply_strain(const vec& epsilon, bool active) const { + if (epsilon.n_elem != 6) { + throw invalid_argument("Strain vector must have 6 elements"); + } + mat R = as_matrix(); + mat ve = fill_voigt_strain(R, active); + return ve * epsilon; +} + +mat::fixed<6,6> Rotation::apply_stiffness(const mat::fixed<6,6>& L, bool active) const { + mat::fixed<6,6> vs = as_voigt_stress_rotation(active); + mat::fixed<6,6> result; + result = vs * L * vs.t(); + return result; +} + +mat Rotation::apply_stiffness(const mat& L, bool active) const { + if (L.n_rows != 6 || L.n_cols != 6) { + throw invalid_argument("Stiffness matrix must be 6x6"); + } + mat R = as_matrix(); + mat vs = fill_voigt_stress(R, active); + return vs * (L * trans(vs)); +} + +mat::fixed<6,6> Rotation::apply_compliance(const mat::fixed<6,6>& M, bool active) const { + mat::fixed<6,6> ve = as_voigt_strain_rotation(active); + mat::fixed<6,6> result; + result = ve * M * ve.t(); + return result; +} + +mat Rotation::apply_compliance(const mat& M, bool active) const { + if (M.n_rows != 6 || M.n_cols != 6) { + throw invalid_argument("Compliance matrix must be 6x6"); + } + mat R = as_matrix(); + mat ve = fill_voigt_strain(R, active); + return ve * (M * trans(ve)); +} + +mat::fixed<6,6> Rotation::apply_strain_concentration(const mat::fixed<6,6>& A, bool active) const { + mat::fixed<6,6> ve = as_voigt_strain_rotation(active); + mat::fixed<6,6> vs = as_voigt_stress_rotation(active); + mat::fixed<6,6> result; + result = ve * mat(A) * vs.t(); + return result; +} + +mat Rotation::apply_strain_concentration(const mat& A, bool active) const { + if (A.n_rows != 6 || A.n_cols != 6) { + throw invalid_argument("Strain concentration tensor must be 6x6"); + } + mat R = as_matrix(); + mat ve = fill_voigt_strain(R, active); + mat vs = fill_voigt_stress(R, active); + return ve * (A * trans(vs)); +} + +mat::fixed<6,6> Rotation::apply_stress_concentration(const mat::fixed<6,6>& B, bool active) const { + mat::fixed<6,6> vs = as_voigt_stress_rotation(active); + mat::fixed<6,6> ve = as_voigt_strain_rotation(active); + mat::fixed<6,6> result; + result = vs * mat(B) * ve.t(); + return result; +} + +mat Rotation::apply_stress_concentration(const mat& B, bool active) const { + if (B.n_rows != 6 || B.n_cols != 6) { + throw invalid_argument("Stress concentration tensor must be 6x6"); + } + mat R = as_matrix(); + mat vs = fill_voigt_stress(R, active); + mat ve = fill_voigt_strain(R, active); + return vs * (B * trans(ve)); +} + +// ============================================================================= +// Rotation Class: Operations +// ============================================================================= + +Rotation Rotation::inv() const { + return Rotation(quat_conjugate(_quat)); +} + +double Rotation::magnitude(bool degrees) const { + double angle = 2.0 * acos(min(abs(_quat(3)), 1.0)); + + if (degrees) { + angle *= 180.0 / simcoon::pi; + } + + return angle; +} + +Rotation Rotation::operator*(const Rotation& other) const { + return Rotation(normalize_quat(quat_multiply(_quat, other._quat))); +} + +Rotation& Rotation::operator*=(const Rotation& other) { + _quat = normalize_quat(quat_multiply(_quat, other._quat)); + return *this; +} + +Rotation Rotation::slerp(const Rotation& other, double t) const { + vec::fixed<4> q1 = _quat; + vec::fixed<4> q2 = other._quat; + + double dot = arma::dot(q1, q2); + + if (dot < 0.0) { + q2 = -q2; + dot = -dot; + } + + dot = min(dot, 1.0); + + vec::fixed<4> result; + + if (dot > 0.9995) { + result = q1 + t * (q2 - q1); + } else { + double theta_0 = acos(dot); + double theta = theta_0 * t; + double sin_theta = sin(theta); + double sin_theta_0 = sin(theta_0); + + double s0 = cos(theta) - dot * sin_theta / sin_theta_0; + double s1 = sin_theta / sin_theta_0; + + result = s0 * q1 + s1 * q2; + } + + return Rotation(normalize_quat(result)); +} + +// ============================================================================= +// Convenience Free Functions +// ============================================================================= + +vec rotate_vec(const vec &v, const mat &DR) { + return DR*v; +} + +vec rotate_vec(const vec &v, const double &alpha, const int &axis) { + assert(axis>0); + Rotation rot = Rotation::from_axis_angle(alpha, axis); + return rot.apply(v); +} + +mat rotate_mat(const mat &m, const mat &DR) { + return trans(DR)*m*DR; +} + +mat rotate_mat(const mat &m, const double &alpha, const int &axis) { + Rotation rot = Rotation::from_axis_angle(alpha, axis); + return rot.apply_tensor(m); +} + +vec rotate_stress(const vec &V, const double &alpha, const int &axis, const bool &active) { + mat vs = fill_voigt_stress(alpha, axis, active); + return vs*V; } -//To rotate a stress vector (6) vec rotate_stress(const vec &V, const mat &DR, const bool &active) { - - mat QS = fillQS(DR, active); - return QS*V; + mat vs = fill_voigt_stress(DR, active); + return vs*V; } - -//To rotate a strain vector (6) -vec rotate_strain(const vec &V, const double &alpha, const int &axis, const bool &active) { - mat QE = fillQE(alpha, axis, active); - return QE*V; +vec rotate_strain(const vec &V, const double &alpha, const int &axis, const bool &active) { + mat ve = fill_voigt_strain(alpha, axis, active); + return ve*V; } vec rotate_strain(const vec &V, const mat &DR, const bool &active) { - - mat QE = fillQE(DR, active); - return QE*V; -} - - -//To rotate from local to global a strain tensor (6) using Euler angles -mat rotate_l2g_strain(const vec &E, const double &psi, const double &theta, const double &phi) { - - mat E_temp = E; - if(fabs(phi) > simcoon::iota) { - E_temp = rotate_strain(E_temp, -phi, simcoon::axis_phi, false); - } - if(fabs(theta) > simcoon::iota) { - E_temp = rotate_strain(E_temp, -theta, simcoon::axis_theta, false); - } - if(fabs(psi) > simcoon::iota) { - E_temp = rotate_strain(E_temp, -psi, simcoon::axis_psi, false); - } - - return E_temp; -} - -//To rotate from global to local a strain tensor (6) using Euler angles -mat rotate_g2l_strain(const vec &E, const double &psi, const double &theta, const double &phi) { - - mat E_temp = E; - if(fabs(psi) > simcoon::iota) { - E_temp = rotate_strain(E_temp, psi, simcoon::axis_psi, false); - } - if(fabs(theta) > simcoon::iota) { - E_temp = rotate_strain(E_temp, theta, simcoon::axis_theta, false); - } - if(fabs(phi) > simcoon::iota) { - E_temp = rotate_strain(E_temp, phi, simcoon::axis_phi, false); - } - - return E_temp; -} - -//To rotate from local to global a stress tensor (6) -mat rotate_l2g_stress(const vec &S, const double &psi, const double &theta, const double &phi) { - - mat S_temp = S; - if(fabs(phi) > simcoon::iota) { - S_temp = rotate_stress(S_temp, -phi, simcoon::axis_phi, false); - } - if(fabs(theta) > simcoon::iota) { - S_temp = rotate_stress(S_temp, -theta, simcoon::axis_theta, false); - } - if(fabs(psi) > simcoon::iota) { - S_temp = rotate_stress(S_temp, -psi, simcoon::axis_psi, false); - } - - return S_temp; -} - -//To rotate from global to local a stress tensor (6) -mat rotate_g2l_stress(const vec &S, const double &psi, const double &theta, const double &phi) { - - mat S_temp = S; - if(fabs(psi) > simcoon::iota) { - S_temp = rotate_stress(S_temp, psi, simcoon::axis_psi, false); - } - if(fabs(theta) > simcoon::iota) { - S_temp = rotate_stress(S_temp, theta, simcoon::axis_theta, false); - } - if(fabs(phi) > simcoon::iota) { - S_temp = rotate_stress(S_temp, phi, simcoon::axis_phi, false); - } - - return S_temp; -} - -//To rotate from local to global a stiffness matrix (6,6) -mat rotate_l2g_L(const mat &Lt, const double &psi, const double &theta, const double &phi) { - - mat Lt_temp = Lt; - if(fabs(phi) > simcoon::iota) { - Lt_temp = rotateL(Lt_temp, -phi, simcoon::axis_phi, false); - } - if(fabs(theta) > simcoon::iota) { - Lt_temp = rotateL(Lt_temp, -theta, simcoon::axis_theta, false); - } - if(fabs(psi) > simcoon::iota) { - Lt_temp = rotateL(Lt_temp, -psi, simcoon::axis_psi, false); - } - - return Lt_temp; -} - -//To rotate from global to local a stiffness matrix (6,6) -mat rotate_g2l_L(const mat &Lt, const double &psi, const double &theta, const double &phi) { - - mat Lt_temp = Lt; - if(fabs(psi) > simcoon::iota) { - Lt_temp = rotateL(Lt_temp, psi, simcoon::axis_psi, false); - } - if(fabs(theta) > simcoon::iota) { - Lt_temp = rotateL(Lt_temp, theta, simcoon::axis_theta, false); - } - if(fabs(phi) > simcoon::iota) { - Lt_temp = rotateL(Lt_temp, phi, simcoon::axis_phi, false); - } - - return Lt_temp; -} - -//To rotate from local to global a localisation matrix (6,6) -mat rotate_l2g_M(const mat &M, const double &psi, const double &theta, const double &phi) { - - mat M_temp = M; - if(fabs(phi) > simcoon::iota) { - M_temp = rotateM(M_temp, -phi, simcoon::axis_phi, false); - } - if(fabs(theta) > simcoon::iota) { - M_temp = rotateM(M_temp, -theta, simcoon::axis_theta, false); - } - if(fabs(psi) > simcoon::iota) { - M_temp = rotateM(M_temp, -psi, simcoon::axis_psi, false); - } - - return M_temp; -} - -//To rotate from global to local a localisation matrix (6,6) -mat rotate_g2l_M(const mat &M, const double &psi, const double &theta, const double &phi) { - - mat M_temp = M; - if(fabs(psi) > simcoon::iota) { - M_temp = rotateM(M_temp, psi, simcoon::axis_psi, false); - } - if(fabs(theta) > simcoon::iota) { - M_temp = rotateM(M_temp, theta, simcoon::axis_theta, false); - } - if(fabs(phi) > simcoon::iota) { - M_temp = rotateM(M_temp, phi, simcoon::axis_phi, false); - } - - return M_temp; -} - -//To rotate from local to global a strain localisation matrix (6,6) -mat rotate_l2g_A(const mat &A, const double &psi, const double &theta, const double &phi) { - - mat A_temp = A; - if(fabs(phi) > simcoon::iota) { - A_temp = rotateA(A_temp, -phi, simcoon::axis_phi, false); - } - if(fabs(theta) > simcoon::iota) { - A_temp = rotateA(A_temp, -theta, simcoon::axis_theta, false); - } - if(fabs(psi) > simcoon::iota) { - A_temp = rotateA(A_temp, -psi, simcoon::axis_psi, false); - } - - return A_temp; -} - -//To rotate from global to local a strain localisation matrix (6,6) -mat rotate_g2l_A(const mat &A, const double &psi, const double &theta, const double &phi) { - - mat A_temp = A; - if(fabs(psi) > simcoon::iota) { - A_temp = rotateA(A_temp, psi, simcoon::axis_psi, false); - } - if(fabs(theta) > simcoon::iota) { - A_temp = rotateA(A_temp, theta, simcoon::axis_theta, false); - } - if(fabs(phi) > simcoon::iota) { - A_temp = rotateA(A_temp, phi, simcoon::axis_phi, false); - } - - return A_temp; -} - -//To rotate from local to global a stress localisation matrix (6,6) -mat rotate_l2g_B(const mat &B, const double &psi, const double &theta, const double &phi) { - - mat B_temp = B; - if(fabs(phi) > simcoon::iota) { - B_temp = rotateB(B_temp, -phi, simcoon::axis_phi, false); - } - if(fabs(theta) > simcoon::iota) { - B_temp = rotateB(B_temp, -theta, simcoon::axis_theta, false); - } - if(fabs(psi) > simcoon::iota) { - B_temp = rotateB(B_temp, -psi, simcoon::axis_psi, false); - } - - return B_temp; -} - -//To rotate from global to local a stress localisation matrix (6,6) -mat rotate_g2l_B(const mat &B, const double &psi, const double &theta, const double &phi) { - - mat B_temp = B; - if(fabs(psi) > simcoon::iota) { - B_temp = rotateB(B_temp, psi, simcoon::axis_psi, false); - } - if(fabs(theta) > simcoon::iota) { - B_temp = rotateB(B_temp, theta, simcoon::axis_theta, false); - } - if(fabs(phi) > simcoon::iota) { - B_temp = rotateB(B_temp, phi, simcoon::axis_phi, false); - } - - return B_temp; -} - + mat ve = fill_voigt_strain(DR, active); + return ve*V; +} + +mat rotate_stiffness(const mat &L, const double &alpha, const int &axis, const bool &active) { + mat vs = fill_voigt_stress(alpha, axis, active); + return vs*(L*trans(vs)); +} + +mat rotate_stiffness(const mat &L, const mat &DR, const bool &active) { + mat vs = fill_voigt_stress(DR, active); + return vs*(L*trans(vs)); +} + +mat rotate_compliance(const mat &M, const double &alpha, const int &axis, const bool &active) { + mat ve = fill_voigt_strain(alpha, axis, active); + return ve*(M*trans(ve)); +} + +mat rotate_compliance(const mat &M, const mat &DR, const bool &active) { + mat ve = fill_voigt_strain(DR, active); + return ve*(M*trans(ve)); +} + +mat rotate_strain_concentration(const mat &A, const double &alpha, const int &axis, const bool &active) { + mat ve = fill_voigt_strain(alpha, axis, active); + mat vs = fill_voigt_stress(alpha, axis, active); + return ve*(A*trans(vs)); +} + +mat rotate_strain_concentration(const mat &A, const mat &DR, const bool &active) { + mat ve = fill_voigt_strain(DR, active); + mat vs = fill_voigt_stress(DR, active); + return ve*(A*trans(vs)); +} + +mat rotate_stress_concentration(const mat &B, const double &alpha, const int &axis, const bool &active) { + mat ve = fill_voigt_strain(alpha, axis, active); + mat vs = fill_voigt_stress(alpha, axis, active); + return vs*(B*trans(ve)); +} + +mat rotate_stress_concentration(const mat &B, const mat &DR, const bool &active) { + mat ve = fill_voigt_strain(DR, active); + mat vs = fill_voigt_stress(DR, active); + return vs*(B*trans(ve)); +} + } //namespace simcoon diff --git a/src/Simulation/Phase/state_variables.cpp b/src/Simulation/Phase/state_variables.cpp index 5481c5d38..7b8730357 100755 --- a/src/Simulation/Phase/state_variables.cpp +++ b/src/Simulation/Phase/state_variables.cpp @@ -474,7 +474,7 @@ arma::mat state_variables::Biot_stress_start() state_variables& state_variables::rotate_l2g(const state_variables& sv, const double &psi, const double &theta, const double &phi) //---------------------------------------------------------------------- { - + Etot = sv.Etot; DEtot = sv.DEtot; etot = sv.etot; @@ -493,66 +493,31 @@ state_variables& state_variables::rotate_l2g(const state_variables& sv, const do DR = sv.DR; T = sv.T; DT = sv.DT; - + nstatev = sv.nstatev; statev = sv.statev; statev_start = sv.statev_start; - - if(fabs(phi) > simcoon::iota) { - Etot = rotate_strain(Etot, -phi, simcoon::axis_phi); - DEtot = rotate_strain(DEtot, -phi, simcoon::axis_phi); - etot = rotate_strain(etot, -phi, simcoon::axis_phi); - Detot = rotate_strain(Detot, -phi, simcoon::axis_phi); - PKII = rotate_stress(PKII, -phi, simcoon::axis_phi); - PKII_start = rotate_stress(PKII_start, -phi, simcoon::axis_phi); - tau = rotate_stress(tau, -phi, simcoon::axis_phi); - tau_start = rotate_stress(tau_start, -phi, simcoon::axis_phi); - sigma = rotate_stress(sigma, -phi, simcoon::axis_phi); - sigma_start = rotate_stress(sigma_start, -phi, simcoon::axis_phi); - F0 = rotate_mat(F0, -phi, simcoon::axis_phi); - F1 = rotate_mat(F1, -phi, simcoon::axis_phi); - U0 = rotate_mat(U0, -phi, simcoon::axis_phi); - U1 = rotate_mat(U1, -phi, simcoon::axis_phi); - R = rotate_mat(R, -phi, simcoon::axis_phi); - DR = rotate_mat(DR, -phi, simcoon::axis_phi); - } - if(fabs(theta) > simcoon::iota) { - Etot = rotate_strain(Etot, -theta, simcoon::axis_theta); - DEtot = rotate_strain(DEtot, -theta, simcoon::axis_theta); - etot = rotate_strain(etot, -theta, simcoon::axis_theta); - Detot = rotate_strain(Detot, -theta, simcoon::axis_theta); - PKII = rotate_stress(PKII, -theta, simcoon::axis_theta); - PKII_start = rotate_stress(PKII_start, -theta, simcoon::axis_theta); - tau = rotate_stress(tau, -theta, simcoon::axis_theta); - tau_start = rotate_stress(tau_start, -theta, simcoon::axis_theta); - sigma = rotate_stress(sigma, -theta, simcoon::axis_theta); - sigma_start = rotate_stress(sigma_start, -theta, simcoon::axis_theta); - F0 = rotate_mat(F0, -theta, simcoon::axis_theta); - F1 = rotate_mat(F1, -theta, simcoon::axis_theta); - U0 = rotate_mat(U0, -theta, simcoon::axis_theta); - U1 = rotate_mat(U1, -theta, simcoon::axis_theta); - R = rotate_mat(R, -theta, simcoon::axis_theta); - DR = rotate_mat(DR, -theta, simcoon::axis_theta); - } - if(fabs(psi) > simcoon::iota) { - Etot = rotate_strain(Etot, -psi, simcoon::axis_psi); - DEtot = rotate_strain(DEtot, -psi, simcoon::axis_psi); - etot = rotate_strain(etot, -psi, simcoon::axis_psi); - Detot = rotate_strain(Detot, -psi, simcoon::axis_psi); - PKII = rotate_stress(PKII, -psi, simcoon::axis_psi); - PKII_start = rotate_stress(PKII_start, -psi, simcoon::axis_psi); - tau = rotate_stress(tau, -psi, simcoon::axis_psi); - tau_start = rotate_stress(tau_start, -psi, simcoon::axis_psi); - sigma = rotate_stress(sigma, -psi, simcoon::axis_psi); - sigma_start = rotate_stress(sigma_start, -psi, simcoon::axis_psi); - F0 = rotate_mat(F0, -psi, simcoon::axis_psi); - F1 = rotate_mat(F1, -psi, simcoon::axis_psi); - U0 = rotate_mat(U0, -psi, simcoon::axis_psi); - U1 = rotate_mat(U1, -psi, simcoon::axis_psi); - R = rotate_mat(R, -psi, simcoon::axis_psi); - DR = rotate_mat(DR, -psi, simcoon::axis_psi); - } - + + Rotation rot = Rotation::from_euler(psi, theta, phi, "zxz"); + if (!rot.is_identity()) { + Etot = rot.apply_strain(Etot); + DEtot = rot.apply_strain(DEtot); + etot = rot.apply_strain(etot); + Detot = rot.apply_strain(Detot); + PKII = rot.apply_stress(PKII); + PKII_start = rot.apply_stress(PKII_start); + tau = rot.apply_stress(tau); + tau_start = rot.apply_stress(tau_start); + sigma = rot.apply_stress(sigma); + sigma_start = rot.apply_stress(sigma_start); + F0 = rot.apply_tensor(F0); + F1 = rot.apply_tensor(F1); + U0 = rot.apply_tensor(U0); + U1 = rot.apply_tensor(U1); + R = rot.apply_tensor(R); + DR = rot.apply_tensor(DR); + } + return *this; } @@ -575,69 +540,34 @@ state_variables& state_variables::rotate_g2l(const state_variables& sv, const do R = sv.R; DR = sv.DR; U0 = sv.U0; - U1 = sv.U1; + U1 = sv.U1; T = sv.T; DT = sv.DT; - + nstatev = sv.nstatev; statev = sv.statev; - statev_start = sv.statev_start; - - if(fabs(psi) > simcoon::iota) { - Etot = rotate_strain(Etot, psi, simcoon::axis_psi); - DEtot = rotate_strain(DEtot, psi, simcoon::axis_psi); - etot = rotate_strain(etot, psi, simcoon::axis_psi); - Detot = rotate_strain(Detot, psi, simcoon::axis_psi); - PKII = rotate_stress(PKII, psi, simcoon::axis_psi); - PKII_start = rotate_stress(PKII_start, psi, simcoon::axis_psi); - tau = rotate_stress(tau, psi, simcoon::axis_psi); - tau_start = rotate_stress(tau_start, psi, simcoon::axis_psi); - sigma = rotate_stress(sigma, psi, simcoon::axis_psi); - sigma_start = rotate_stress(sigma_start, psi, simcoon::axis_psi); - F0 = rotate_mat(F0, psi, simcoon::axis_psi); - F1 = rotate_mat(F1, psi, simcoon::axis_psi); - U0 = rotate_mat(U0, psi, simcoon::axis_psi); - U1 = rotate_mat(U1, psi, simcoon::axis_psi); - R = rotate_mat(R, psi, simcoon::axis_psi); - DR = rotate_mat(DR, psi, simcoon::axis_psi); - } - if(fabs(theta) > simcoon::iota) { - Etot = rotate_strain(Etot, theta, simcoon::axis_theta); - DEtot = rotate_strain(DEtot, theta, simcoon::axis_theta); - etot = rotate_strain(etot, theta, simcoon::axis_theta); - Detot = rotate_strain(Detot, theta, simcoon::axis_theta); - PKII = rotate_stress(PKII, theta, simcoon::axis_theta); - PKII_start = rotate_stress(PKII_start, theta, simcoon::axis_theta); - tau = rotate_stress(tau, theta, simcoon::axis_theta); - tau_start = rotate_stress(tau_start, theta, simcoon::axis_theta); - sigma = rotate_stress(sigma, theta, simcoon::axis_theta); - sigma_start = rotate_stress(sigma_start, theta, simcoon::axis_theta); - F0 = rotate_mat(F0, theta, simcoon::axis_theta); - F1 = rotate_mat(F1, theta, simcoon::axis_theta); - U0 = rotate_mat(U0, theta, simcoon::axis_theta); - U1 = rotate_mat(U1, theta, simcoon::axis_theta); - R = rotate_mat(R, theta, simcoon::axis_theta); - DR = rotate_mat(DR, theta, simcoon::axis_theta); - } - if(fabs(phi) > simcoon::iota) { - Etot = rotate_strain(Etot, phi, simcoon::axis_phi); - DEtot = rotate_strain(DEtot, phi, simcoon::axis_phi); - etot = rotate_strain(etot, phi, simcoon::axis_phi); - Detot = rotate_strain(Detot, phi, simcoon::axis_phi); - PKII = rotate_stress(PKII, phi, simcoon::axis_phi); - PKII_start = rotate_stress(PKII_start, phi, simcoon::axis_phi); - tau = rotate_stress(tau, phi, simcoon::axis_phi); - tau_start = rotate_stress(tau_start, phi, simcoon::axis_phi); - sigma = rotate_stress(sigma, phi, simcoon::axis_phi); - sigma_start = rotate_stress(sigma_start, phi, simcoon::axis_phi); - F0 = rotate_mat(F0, phi, simcoon::axis_phi); - F1 = rotate_mat(F1, phi, simcoon::axis_phi); - U0 = rotate_mat(U0, phi, simcoon::axis_phi); - U1 = rotate_mat(U1, phi, simcoon::axis_phi); - R = rotate_mat(R, phi, simcoon::axis_phi); - DR = rotate_mat(DR, phi, simcoon::axis_phi); + statev_start = sv.statev_start; + + Rotation rot = Rotation::from_euler(psi, theta, phi, "zxz").inv(); + if (!rot.is_identity()) { + Etot = rot.apply_strain(Etot); + DEtot = rot.apply_strain(DEtot); + etot = rot.apply_strain(etot); + Detot = rot.apply_strain(Detot); + PKII = rot.apply_stress(PKII); + PKII_start = rot.apply_stress(PKII_start); + tau = rot.apply_stress(tau); + tau_start = rot.apply_stress(tau_start); + sigma = rot.apply_stress(sigma); + sigma_start = rot.apply_stress(sigma_start); + F0 = rot.apply_tensor(F0); + F1 = rot.apply_tensor(F1); + U0 = rot.apply_tensor(U0); + U1 = rot.apply_tensor(U1); + R = rot.apply_tensor(R); + DR = rot.apply_tensor(DR); } - + return *this; } diff --git a/src/Simulation/Phase/state_variables_M.cpp b/src/Simulation/Phase/state_variables_M.cpp index 996947c44..e4f00ba0a 100755 --- a/src/Simulation/Phase/state_variables_M.cpp +++ b/src/Simulation/Phase/state_variables_M.cpp @@ -237,37 +237,26 @@ void state_variables_M::set_start(const int &corate_type) state_variables_M& state_variables_M::rotate_l2g(const state_variables_M& sv, const double &psi, const double &theta, const double &phi) //---------------------------------------------------------------------- { - + state_variables::rotate_l2g(sv, psi, theta, phi); - + sigma_in = sv.sigma_in; sigma_in_start = sv.sigma_in_start; - + Wm = sv.Wm; Wm_start = sv.Wm_start; - + L = sv.L; Lt = sv.Lt; - - if(fabs(phi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, -phi, simcoon::axis_phi); - sigma_in_start = rotate_stress(sigma_in_start, -phi, simcoon::axis_phi); - L = rotateL(L, -phi, simcoon::axis_phi); - Lt = rotateL(Lt, -phi, simcoon::axis_phi); - } - if(fabs(theta) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, -theta, simcoon::axis_theta); - sigma_in_start = rotate_stress(sigma_in_start, -theta, simcoon::axis_theta); - Lt = rotateL(Lt, -theta, simcoon::axis_theta); - L = rotateL(L, -theta, simcoon::axis_theta); - } - if(fabs(psi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, -psi, simcoon::axis_psi); - sigma_in_start = rotate_stress(sigma_in_start, -psi, simcoon::axis_psi); - Lt = rotateL(Lt, -psi, simcoon::axis_psi); - L = rotateL(L, -psi, simcoon::axis_psi); - } - + + Rotation rot = Rotation::from_euler(psi, theta, phi, "zxz"); + if (!rot.is_identity()) { + sigma_in = rot.apply_stress(sigma_in); + sigma_in_start = rot.apply_stress(sigma_in_start); + L = rot.apply_stiffness(L); + Lt = rot.apply_stiffness(Lt); + } + return *this; } @@ -277,35 +266,24 @@ state_variables_M& state_variables_M::rotate_g2l(const state_variables_M& sv, co { state_variables::rotate_g2l(sv, psi, theta, phi); - + sigma_in = sv.sigma_in; sigma_in_start = sv.sigma_in_start; - + Wm = sv.Wm; Wm_start = sv.Wm_start; L = sv.L; Lt = sv.Lt; - - if(fabs(psi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, psi, simcoon::axis_psi); - sigma_in_start = rotate_stress(sigma_in_start, psi, simcoon::axis_psi); - L = rotateL(L, psi, simcoon::axis_psi); - Lt = rotateL(Lt, psi, simcoon::axis_psi); - } - if(fabs(theta) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, theta, simcoon::axis_theta); - sigma_in_start = rotate_stress(sigma_in_start, theta, simcoon::axis_theta); - L = rotateL(L, theta, simcoon::axis_theta); - Lt = rotateL(Lt, theta, simcoon::axis_theta); - } - if(fabs(phi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, phi, simcoon::axis_phi); - sigma_in_start = rotate_stress(sigma_in_start, phi, simcoon::axis_phi); - L = rotateL(L, phi, simcoon::axis_phi); - Lt = rotateL(Lt, phi, simcoon::axis_phi); + + Rotation rot = Rotation::from_euler(psi, theta, phi, "zxz").inv(); + if (!rot.is_identity()) { + sigma_in = rot.apply_stress(sigma_in); + sigma_in_start = rot.apply_stress(sigma_in_start); + L = rot.apply_stiffness(L); + Lt = rot.apply_stiffness(Lt); } - + return *this; } diff --git a/src/Simulation/Phase/state_variables_T.cpp b/src/Simulation/Phase/state_variables_T.cpp index f19b9912a..d86b414c1 100755 --- a/src/Simulation/Phase/state_variables_T.cpp +++ b/src/Simulation/Phase/state_variables_T.cpp @@ -320,32 +320,19 @@ state_variables_T& state_variables_T::rotate_l2g(const state_variables_T& sv, co Wt = sv.Wt; Wm_start = sv.Wm_start; Wt_start = sv.Wt_start; - - if(fabs(phi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, -phi, simcoon::axis_phi); - sigma_in_start = rotate_stress(sigma_in_start, -phi, simcoon::axis_phi); - dSdE = rotateL(dSdE, -phi, simcoon::axis_phi); - dSdEt = rotateL(dSdEt, -phi, simcoon::axis_phi); - dSdT = rotate_stress(dSdT, -phi, simcoon::axis_phi); - drdE = rotate_strain(drdE, -phi, simcoon::axis_phi); - } - if(fabs(theta) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, -theta, simcoon::axis_theta); - sigma_in_start = rotate_stress(sigma_in_start, -theta, simcoon::axis_theta); - dSdE = rotateL(dSdE, -theta, simcoon::axis_theta); - dSdEt = rotateL(dSdEt, -theta, simcoon::axis_theta); - dSdT = rotate_stress(dSdT, -theta, simcoon::axis_theta); - drdE = rotate_strain(drdE, -theta, simcoon::axis_theta); - } - if(fabs(psi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, -psi, simcoon::axis_psi); - sigma_in_start = rotate_stress(sigma_in_start, -psi, simcoon::axis_psi); - dSdE = rotateL(dSdE, -psi, simcoon::axis_psi); - dSdEt = rotateL(dSdEt, -psi, simcoon::axis_psi); - dSdT = rotate_stress(dSdT, -psi, simcoon::axis_psi); - drdE = rotate_strain(drdE, -psi, simcoon::axis_psi); - } - + + Rotation rot = Rotation::from_euler(psi, theta, phi, "zxz"); + if (!rot.is_identity()) { + sigma_in = rot.apply_stress(sigma_in); + sigma_in_start = rot.apply_stress(sigma_in_start); + dSdE = rot.apply_stiffness(dSdE); + dSdEt = rot.apply_stiffness(dSdEt); + vec dSdT_v = rot.apply_stress(vec(dSdT.as_col())); + dSdT = dSdT_v.as_row(); + vec drdE_v = rot.apply_strain(vec(drdE.as_col())); + drdE = drdE_v.as_row(); + } + return *this; } @@ -355,10 +342,10 @@ state_variables_T& state_variables_T::rotate_g2l(const state_variables_T& sv, co { state_variables::rotate_g2l(sv, psi, theta, phi); - + sigma_in = sv.sigma_in; sigma_in_start = sv.sigma_in_start; - + dSdE = sv.dSdE; dSdEt = sv.dSdEt; dSdT = sv.dSdT; @@ -370,33 +357,19 @@ state_variables_T& state_variables_T::rotate_g2l(const state_variables_T& sv, co Wt = sv.Wt; Wm_start = sv.Wm_start; Wt_start = sv.Wt_start; - - if(fabs(psi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, psi, simcoon::axis_psi); - sigma_in_start = rotate_stress(sigma_in_start, psi, simcoon::axis_psi); - dSdE = rotateL(dSdE, psi, simcoon::axis_psi); - dSdEt = rotateL(dSdEt, psi, simcoon::axis_psi); - dSdT = rotate_stress(dSdT, psi, simcoon::axis_psi); - drdE = rotate_strain(drdE, psi, simcoon::axis_psi); - - } - if(fabs(theta) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, theta, simcoon::axis_theta); - sigma_in_start = rotate_stress(sigma_in_start, theta, simcoon::axis_theta); - dSdE = rotateL(dSdE, theta, simcoon::axis_theta); - dSdEt = rotateL(dSdEt, theta, simcoon::axis_theta); - dSdT = rotate_stress(dSdT, theta, simcoon::axis_theta); - drdE = rotate_strain(drdE, theta, simcoon::axis_theta); - } - if(fabs(phi) > simcoon::iota) { - sigma_in = rotate_stress(sigma_in, phi, simcoon::axis_phi); - sigma_in_start = rotate_stress(sigma_in_start, phi, simcoon::axis_phi); - dSdE = rotateL(dSdE, phi, simcoon::axis_phi); - dSdEt = rotateL(dSdEt, phi, simcoon::axis_phi); - dSdT = rotate_stress(dSdT, phi, simcoon::axis_phi); - drdE = rotate_strain(drdE, phi, simcoon::axis_phi); + + Rotation rot = Rotation::from_euler(psi, theta, phi, "zxz").inv(); + if (!rot.is_identity()) { + sigma_in = rot.apply_stress(sigma_in); + sigma_in_start = rot.apply_stress(sigma_in_start); + dSdE = rot.apply_stiffness(dSdE); + dSdEt = rot.apply_stiffness(dSdEt); + vec dSdT_v = rot.apply_stress(vec(dSdT.as_col())); + dSdT = dSdT_v.as_row(); + vec drdE_v = rot.apply_strain(vec(drdE.as_col())); + drdE = drdE_v.as_row(); } - + return *this; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ef6bf848d..42fb4fb1a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -27,6 +27,7 @@ set(TEST_SRCS test/Libraries/Maths/Tlagrange.cpp test/Libraries/Maths/Tnum_solve.cpp test/Libraries/Maths/Trotation.cpp + test/Libraries/Maths/Trotation_class.cpp test/Libraries/Maths/Tstats.cpp # Libraries - Phase diff --git a/test/Libraries/Maths/Trotation.cpp b/test/Libraries/Maths/Trotation.cpp index 4d8e3bceb..f888a7a34 100755 --- a/test/Libraries/Maths/Trotation.cpp +++ b/test/Libraries/Maths/Trotation.cpp @@ -46,10 +46,10 @@ TEST(Trotation, basic_rotation) mat R_zxz = R3*R2*R1; mat R_zyz = R3*R2b*R1; - mat R_zxz_a = fillR(psi, theta, phi, true, "zxz"); - mat R_zxz_p = fillR(psi, theta, phi, false, "zxz"); - mat R_zyz_a = fillR(psi, theta, phi, true, "zyz"); - mat R_zyz_p = fillR(psi, theta, phi, false, "zyz"); + mat R_zxz_a = mat(Rotation::from_euler(psi, theta, phi, "zxz", true).as_matrix()); + mat R_zxz_p = mat(Rotation::from_euler(psi, theta, phi, "zxz", true).as_matrix()).t(); + mat R_zyz_a = mat(Rotation::from_euler(psi, theta, phi, "zyz", true).as_matrix()); + mat R_zyz_p = mat(Rotation::from_euler(psi, theta, phi, "zyz", true).as_matrix()).t(); cout << "R_zxz" << R_zxz << endl; cout << "R_zxz_a" << R_zxz_a << endl; @@ -100,19 +100,19 @@ TEST(Trotation, rotation) //test of rotate A mat S_c = Eshelby_cylinder(0.12); - mat S_c1 = rotateA(S_c, psi, 3); - mat S_c2 = rotateA(S_c, theta, 1); - mat S_c3 = rotateA(S_c, phi, 3); + mat S_c1 = rotate_strain_concentration(S_c, psi, 3); + mat S_c2 = rotate_strain_concentration(S_c, theta, 1); + mat S_c3 = rotate_strain_concentration(S_c, phi, 3); - mat S_c11 = rotateA(S_c,R1); - mat S_c22 = rotateA(S_c,R2); - mat S_c33 = rotateA(S_c,R3); + mat S_c11 = rotate_strain_concentration(S_c,R1); + mat S_c22 = rotate_strain_concentration(S_c,R2); + mat S_c33 = rotate_strain_concentration(S_c,R3); - mat S_cR = rotateA(S_c, R1); - S_cR = rotateA(S_cR, R2); - S_cR = rotateA(S_cR, R3); + mat S_cR = rotate_strain_concentration(S_c, R1); + S_cR = rotate_strain_concentration(S_cR, R2); + S_cR = rotate_strain_concentration(S_cR, R3); - mat S_cRR = rotateA(S_c,R); + mat S_cRR = rotate_strain_concentration(S_c,R); EXPECT_LT( norm(S_c33-S_c3,2),1.E-9); EXPECT_LT( norm(S_c22-S_c2,2),1.E-9); @@ -129,8 +129,8 @@ TEST(Trotation, rotation) assert(axis_theta == 1); assert(axis_phi == 3); - mat Rp2 = fillR(psi,theta,phi); - mat Rp3 = fillR(psi,theta,phi, true, "zxz"); + mat Rp2 = mat(Rotation::from_euler(psi,theta,phi, "user", true).as_matrix()); + mat Rp3 = mat(Rotation::from_euler(psi,theta,phi, "zxz", true).as_matrix()); mat a1 = rotate_mat(a, R); mat a2 = rotate_mat(a, Rp2); diff --git a/test/Libraries/Maths/Trotation_class.cpp b/test/Libraries/Maths/Trotation_class.cpp new file mode 100644 index 000000000..0a7cd4cb1 --- /dev/null +++ b/test/Libraries/Maths/Trotation_class.cpp @@ -0,0 +1,685 @@ +/* This file is part of simcoon. + + simcoon is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + simcoon is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with simcoon. If not, see . + + */ + +///@file Trotation_class.cpp +///@brief Test for the Rotation class +///@version 1.0 + +#include +#include + +#include +#include + +using namespace std; +using namespace arma; +using namespace simcoon; + +// ============================================================================= +// Identity and Basic Tests +// ============================================================================= + +TEST(TRotationClass, identity) +{ + Rotation r = Rotation::identity(); + + EXPECT_TRUE(r.is_identity()); + + // Identity quaternion should be [0, 0, 0, 1] + vec::fixed<4> q = r.as_quat(); + EXPECT_NEAR(q(0), 0.0, 1.E-12); + EXPECT_NEAR(q(1), 0.0, 1.E-12); + EXPECT_NEAR(q(2), 0.0, 1.E-12); + EXPECT_NEAR(q(3), 1.0, 1.E-12); + + // Identity matrix should be eye(3,3) + mat::fixed<3,3> R = r.as_matrix(); + EXPECT_LT(norm(R - eye(3,3), "fro"), 1.E-12); +} + +TEST(TRotationClass, default_constructor) +{ + Rotation r; + EXPECT_TRUE(r.is_identity()); +} + +// ============================================================================= +// Factory Methods +// ============================================================================= + +TEST(TRotationClass, from_quat) +{ + // Test with a known quaternion (90 degrees around z-axis) + double angle = simcoon::pi / 2.0; + double s = sin(angle / 2.0); + double c = cos(angle / 2.0); + vec::fixed<4> q = {0, 0, s, c}; + + Rotation r = Rotation::from_quat(q); + vec::fixed<4> q_out = r.as_quat(); + + // Should be normalized and equal + EXPECT_NEAR(norm(q_out), 1.0, 1.E-12); + EXPECT_LT(norm(q - q_out), 1.E-12); +} + +TEST(TRotationClass, from_matrix) +{ + // Create a rotation matrix using Rotation class from_euler + double psi = 0.5; + double theta = 0.3; + double phi = 0.7; + + Rotation r_ref = Rotation::from_euler(psi, theta, phi, "zxz"); + mat::fixed<3,3> R_expected = r_ref.as_matrix(); + + Rotation r = Rotation::from_matrix(R_expected); + mat::fixed<3,3> R_actual = r.as_matrix(); + + EXPECT_LT(norm(R_expected - R_actual, "fro"), 1.E-9); +} + +TEST(TRotationClass, from_euler_zxz) +{ + double psi = 23.0 * (simcoon::pi / 180.0); + double theta = 42.0 * (simcoon::pi / 180.0); + double phi = 165.0 * (simcoon::pi / 180.0); + + // Build reference from individual axis rotations (intrinsic zxz) + Rotation rz1 = Rotation::from_axis_angle(psi, 3); + Rotation rx = Rotation::from_axis_angle(theta, 1); + Rotation rz2 = Rotation::from_axis_angle(phi, 3); + Rotation r_ref = rz2 * rx * rz1; + + // Create rotation using from_euler + Rotation r = Rotation::from_euler(psi, theta, phi, "zxz", true, false); + + EXPECT_TRUE(r.equals(r_ref, 1.E-9)); +} + +TEST(TRotationClass, from_euler_zyz) +{ + double psi = 23.0 * (simcoon::pi / 180.0); + double theta = 42.0 * (simcoon::pi / 180.0); + double phi = 165.0 * (simcoon::pi / 180.0); + + // Build reference from individual axis rotations (intrinsic zyz) + Rotation rz1 = Rotation::from_axis_angle(psi, 3); + Rotation ry = Rotation::from_axis_angle(theta, 2); + Rotation rz2 = Rotation::from_axis_angle(phi, 3); + Rotation r_ref = rz2 * ry * rz1; + + Rotation r = Rotation::from_euler(psi, theta, phi, "zyz", true, false); + + EXPECT_TRUE(r.equals(r_ref, 1.E-9)); +} + +TEST(TRotationClass, from_euler_degrees) +{ + // Test with degrees flag + Rotation r1 = Rotation::from_euler(90.0, 0.0, 0.0, "zxz", true, true); + Rotation r2 = Rotation::from_euler(simcoon::pi/2.0, 0.0, 0.0, "zxz", true, false); + + EXPECT_TRUE(r1.equals(r2, 1.E-9)); +} + +TEST(TRotationClass, from_axis_angle) +{ + // Test rotation around z-axis + double angle = simcoon::pi / 4.0; // 45 degrees + + Rotation r = Rotation::from_axis_angle(angle, 3, false); + mat::fixed<3,3> R_actual = r.as_matrix(); + + // Build expected rotation matrix for z-axis rotation + mat R_expected = {{cos(angle), -sin(angle), 0}, {sin(angle), cos(angle), 0}, {0, 0, 1}}; + + EXPECT_LT(norm(R_expected - R_actual, "fro"), 1.E-9); +} + +TEST(TRotationClass, from_rotvec) +{ + // Rotation vector: 45 degrees around z-axis + double angle = simcoon::pi / 4.0; + vec::fixed<3> rotvec = {0, 0, angle}; + + Rotation r = Rotation::from_rotvec(rotvec, false); + + // Should match axis-angle construction + Rotation r_expected = Rotation::from_axis_angle(angle, 3); + + EXPECT_TRUE(r.equals(r_expected, 1.E-9)); +} + +// ============================================================================= +// Conversion Methods +// ============================================================================= + +TEST(TRotationClass, as_euler_roundtrip) +{ + // Test round-trip: euler -> rotation -> euler + double psi_in = 0.5; + double theta_in = 0.3; + double phi_in = 0.7; + + Rotation r = Rotation::from_euler(psi_in, theta_in, phi_in, "zxz", true, false); + vec::fixed<3> euler_out = r.as_euler("zxz", true, false); + + // Reconstruct rotation from output euler angles + Rotation r2 = Rotation::from_euler(euler_out(0), euler_out(1), euler_out(2), "zxz", true, false); + + // The two rotations should be equal + EXPECT_TRUE(r.equals(r2, 1.E-9)); +} + +TEST(TRotationClass, as_rotvec_roundtrip) +{ + // Create a rotation + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + + // Convert to rotvec and back + vec::fixed<3> rotvec = r.as_rotvec(); + Rotation r2 = Rotation::from_rotvec(rotvec); + + EXPECT_TRUE(r.equals(r2, 1.E-9)); +} + +TEST(TRotationClass, as_voigt_stress_rotation) +{ + // Test that voigt stress rotation matrix correctly transforms stress vectors + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + + mat::fixed<6,6> QS = r.as_voigt_stress_rotation(true); + vec::fixed<6> sigma = {100, 50, 0, 25, 0, 0}; + + // QS applied to stress should match apply_stress + vec::fixed<6> sigma_QS = QS * sigma; + vec::fixed<6> sigma_apply = r.apply_stress(sigma, true); + + EXPECT_LT(norm(sigma_QS - sigma_apply), 1.E-9); +} + +TEST(TRotationClass, as_voigt_strain_rotation) +{ + // Test that voigt strain rotation matrix correctly transforms strain vectors + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + + mat::fixed<6,6> QE = r.as_voigt_strain_rotation(true); + vec::fixed<6> epsilon = {0.01, -0.005, -0.005, 0.002, 0.001, 0.003}; + + // QE applied to strain should match apply_strain + vec::fixed<6> epsilon_QE = QE * epsilon; + vec::fixed<6> epsilon_apply = r.apply_strain(epsilon, true); + + EXPECT_LT(norm(epsilon_QE - epsilon_apply), 1.E-9); +} + +// ============================================================================= +// Apply Methods +// ============================================================================= + +TEST(TRotationClass, apply_vector) +{ + // 90 degree rotation around z-axis + Rotation r = Rotation::from_axis_angle(simcoon::pi / 2.0, 3); + + vec::fixed<3> v = {1, 0, 0}; + vec::fixed<3> v_rot = r.apply(v); + + // x-axis should rotate to y-axis + EXPECT_NEAR(v_rot(0), 0.0, 1.E-9); + EXPECT_NEAR(v_rot(1), 1.0, 1.E-9); + EXPECT_NEAR(v_rot(2), 0.0, 1.E-9); +} + +TEST(TRotationClass, apply_vector_inverse) +{ + Rotation r = Rotation::from_axis_angle(simcoon::pi / 2.0, 3); + + vec::fixed<3> v = {1, 0, 0}; + vec::fixed<3> v_rot = r.apply(v, false); + vec::fixed<3> v_back = r.apply(v_rot, true); + + EXPECT_LT(norm(v - v_back), 1.E-9); +} + +TEST(TRotationClass, apply_stress_consistency) +{ + // Test that apply_stress matches existing rotate_stress function + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + mat::fixed<3,3> R = r.as_matrix(); + + vec::fixed<6> sigma = {100, 50, 0, 25, 0, 0}; + + vec::fixed<6> sigma_class = r.apply_stress(sigma, true); + vec sigma_func = rotate_stress(vec(sigma), mat(R), true); + + EXPECT_LT(norm(vec(sigma_class) - sigma_func), 1.E-9); +} + +TEST(TRotationClass, apply_strain_consistency) +{ + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + mat::fixed<3,3> R = r.as_matrix(); + + vec::fixed<6> epsilon = {0.01, -0.005, -0.005, 0.002, 0.001, 0.003}; + + vec::fixed<6> epsilon_class = r.apply_strain(epsilon, true); + vec epsilon_func = rotate_strain(vec(epsilon), mat(R), true); + + EXPECT_LT(norm(vec(epsilon_class) - epsilon_func), 1.E-9); +} + +TEST(TRotationClass, apply_stiffness_consistency) +{ + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + mat::fixed<3,3> R = r.as_matrix(); + + // Create a simple stiffness matrix (isotropic) + mat L = eye(6, 6) * 100.0; + L(0,1) = L(1,0) = 30.0; + L(0,2) = L(2,0) = 30.0; + L(1,2) = L(2,1) = 30.0; + + mat::fixed<6,6> L_fixed; + L_fixed = L; + + mat::fixed<6,6> L_class = r.apply_stiffness(L_fixed, true); + mat L_func = rotate_stiffness(L, mat(R), true); + + EXPECT_LT(norm(mat(L_class) - L_func, "fro"), 1.E-9); +} + +// ============================================================================= +// Operations +// ============================================================================= + +TEST(TRotationClass, inverse) +{ + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + Rotation r_inv = r.inv(); + + // r * r_inv should be identity + Rotation product = r * r_inv; + EXPECT_TRUE(product.is_identity(1.E-9)); +} + +TEST(TRotationClass, composition) +{ + // Test that (r1 * r2).apply(v) == r1.apply(r2.apply(v)) + Rotation r1 = Rotation::from_axis_angle(0.5, 1); + Rotation r2 = Rotation::from_axis_angle(0.3, 2); + + vec::fixed<3> v = {1, 2, 3}; + + // Method 1: compose then apply + Rotation r_composed = r1 * r2; + vec::fixed<3> v1 = r_composed.apply(v); + + // Method 2: apply sequentially (r2 first, then r1) + vec::fixed<3> v2 = r1.apply(r2.apply(v)); + + EXPECT_LT(norm(v1 - v2), 1.E-9); +} + +TEST(TRotationClass, composition_inplace) +{ + Rotation r1 = Rotation::from_axis_angle(0.5, 1); + Rotation r2 = Rotation::from_axis_angle(0.3, 2); + + Rotation r_composed = r1 * r2; + + Rotation r1_copy = r1; + r1_copy *= r2; + + EXPECT_TRUE(r_composed.equals(r1_copy, 1.E-9)); +} + +TEST(TRotationClass, magnitude) +{ + // 45 degree rotation + double angle = simcoon::pi / 4.0; + Rotation r = Rotation::from_axis_angle(angle, 3); + + EXPECT_NEAR(r.magnitude(), angle, 1.E-9); + EXPECT_NEAR(r.magnitude(true), 45.0, 1.E-9); // degrees +} + +TEST(TRotationClass, slerp) +{ + Rotation r1 = Rotation::identity(); + Rotation r2 = Rotation::from_axis_angle(simcoon::pi / 2.0, 3); + + // t=0 should give r1 + Rotation r_t0 = r1.slerp(r2, 0.0); + EXPECT_TRUE(r_t0.equals(r1, 1.E-9)); + + // t=1 should give r2 + Rotation r_t1 = r1.slerp(r2, 1.0); + EXPECT_TRUE(r_t1.equals(r2, 1.E-9)); + + // t=0.5 should give 45 degrees rotation + Rotation r_t5 = r1.slerp(r2, 0.5); + EXPECT_NEAR(r_t5.magnitude(), simcoon::pi / 4.0, 1.E-9); +} + +TEST(TRotationClass, equals) +{ + Rotation r1 = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + Rotation r2 = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + Rotation r3 = Rotation::from_euler(0.6, 0.3, 0.7, "zxz"); + + EXPECT_TRUE(r1.equals(r2, 1.E-9)); + EXPECT_FALSE(r1.equals(r3, 1.E-9)); +} + +TEST(TRotationClass, apply_tensor) +{ + // Rotate a symmetric 3x3 tensor + Rotation r = Rotation::from_axis_angle(simcoon::pi / 2.0, 3); + mat::fixed<3,3> R = r.as_matrix(); + + // Diagonal tensor + mat::fixed<3,3> T = {{100, 10, 0}, {10, 50, 0}, {0, 0, 25}}; + + mat::fixed<3,3> T_rot = r.apply_tensor(T); + + // Verify R * T * R^T + mat::fixed<3,3> T_expected; + T_expected = R * T * R.t(); + EXPECT_LT(norm(T_rot - T_expected, "fro"), 1.E-9); + + // Inverse should recover original + mat::fixed<3,3> T_back = r.apply_tensor(T_rot, true); + EXPECT_LT(norm(T - T_back, "fro"), 1.E-9); +} + +TEST(TRotationClass, apply_compliance_consistency) +{ + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + mat::fixed<3,3> R = r.as_matrix(); + + // Create a compliance matrix + mat M = eye(6, 6) * 0.01; + M(0,1) = M(1,0) = -0.003; + M(0,2) = M(2,0) = -0.003; + M(1,2) = M(2,1) = -0.003; + + mat::fixed<6,6> M_fixed; + M_fixed = M; + + mat::fixed<6,6> M_class = r.apply_compliance(M_fixed, true); + mat M_func = rotate_compliance(M, mat(R), true); + + EXPECT_LT(norm(mat(M_class) - M_func, "fro"), 1.E-9); +} + +// ============================================================================= +// Apply Methods (Concentration Tensors) +// ============================================================================= + +TEST(TRotationClass, apply_strain_concentration_consistency) +{ + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + mat::fixed<3,3> R = r.as_matrix(); + + // Create a strain concentration tensor + mat A = eye(6, 6) + 0.1 * randu(6, 6); + mat::fixed<6,6> A_fixed; + A_fixed = A; + + mat::fixed<6,6> A_class = r.apply_strain_concentration(A_fixed, true); + mat A_func = rotate_strain_concentration(A, mat(R), true); + + EXPECT_LT(norm(mat(A_class) - A_func, "fro"), 1.E-9); +} + +TEST(TRotationClass, apply_stress_concentration_consistency) +{ + Rotation r = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + mat::fixed<3,3> R = r.as_matrix(); + + // Create a stress concentration tensor + mat B = eye(6, 6) + 0.1 * randu(6, 6); + mat::fixed<6,6> B_fixed; + B_fixed = B; + + mat::fixed<6,6> B_class = r.apply_stress_concentration(B_fixed, true); + mat B_func = rotate_stress_concentration(B, mat(R), true); + + EXPECT_LT(norm(mat(B_class) - B_func, "fro"), 1.E-9); +} + +// ============================================================================= +// Additional Factory Method Tests +// ============================================================================= + +TEST(TRotationClass, from_axis_angle_x) +{ + double angle = simcoon::pi / 3.0; // 60 degrees around x + Rotation r = Rotation::from_axis_angle(angle, 1); + mat::fixed<3,3> R_actual = r.as_matrix(); + + // Expected: rotation around x-axis + mat R_expected = {{1, 0, 0}, {0, cos(angle), -sin(angle)}, {0, sin(angle), cos(angle)}}; + + EXPECT_LT(norm(R_expected - R_actual, "fro"), 1.E-9); +} + +TEST(TRotationClass, from_axis_angle_y) +{ + double angle = simcoon::pi / 5.0; // 36 degrees around y + Rotation r = Rotation::from_axis_angle(angle, 2); + mat::fixed<3,3> R_actual = r.as_matrix(); + + // Expected: rotation around y-axis + mat R_expected = {{cos(angle), 0, sin(angle)}, {0, 1, 0}, {-sin(angle), 0, cos(angle)}}; + + EXPECT_LT(norm(R_expected - R_actual, "fro"), 1.E-9); +} + +// ============================================================================= +// Additional Euler Convention Tests +// ============================================================================= + +TEST(TRotationClass, from_euler_xyz) +{ + double psi = 0.4, theta = 0.6, phi = 0.8; + + // Build reference from individual axis rotations (extrinsic xyz = intrinsic zyx) + Rotation rx = Rotation::from_axis_angle(psi, 1); + Rotation ry = Rotation::from_axis_angle(theta, 2); + Rotation rz = Rotation::from_axis_angle(phi, 3); + + // Intrinsic xyz: apply x, then y', then z'' + // Quaternion order: rz * ry * rx + Rotation r_ref = rz * ry * rx; + Rotation r_xyz = Rotation::from_euler(psi, theta, phi, "xyz", true, false); + + EXPECT_TRUE(r_ref.equals(r_xyz, 1.E-9)); +} + +TEST(TRotationClass, from_euler_zyx) +{ + double psi = 0.3, theta = 0.5, phi = 0.2; + + Rotation rz = Rotation::from_axis_angle(psi, 3); + Rotation ry = Rotation::from_axis_angle(theta, 2); + Rotation rx = Rotation::from_axis_angle(phi, 1); + + // Intrinsic zyx: apply z, then y', then x'' + Rotation r_ref = rx * ry * rz; + Rotation r_zyx = Rotation::from_euler(psi, theta, phi, "zyx", true, false); + + EXPECT_TRUE(r_ref.equals(r_zyx, 1.E-9)); +} + +TEST(TRotationClass, from_euler_xyx) +{ + double psi = 0.4, theta = 0.7, phi = 0.9; + + Rotation rx1 = Rotation::from_axis_angle(psi, 1); + Rotation ry = Rotation::from_axis_angle(theta, 2); + Rotation rx2 = Rotation::from_axis_angle(phi, 1); + + // Intrinsic xyx: apply x, then y', then x'' + Rotation r_ref = rx2 * ry * rx1; + Rotation r_xyx = Rotation::from_euler(psi, theta, phi, "xyx", true, false); + + EXPECT_TRUE(r_ref.equals(r_xyx, 1.E-9)); +} + +TEST(TRotationClass, euler_roundtrip_xyz) +{ + double psi_in = 0.3, theta_in = 0.5, phi_in = 0.2; + Rotation r = Rotation::from_euler(psi_in, theta_in, phi_in, "xyz", true, false); + vec::fixed<3> angles = r.as_euler("xyz", true, false); + + Rotation r2 = Rotation::from_euler(angles(0), angles(1), angles(2), "xyz", true, false); + EXPECT_TRUE(r.equals(r2, 1.E-9)); +} + +TEST(TRotationClass, euler_extrinsic_vs_intrinsic) +{ + double a = 0.4, b = 0.6, c = 0.8; + + // Extrinsic XYZ = Intrinsic ZYX (reversed) + Rotation r_ext = Rotation::from_euler(a, b, c, "xyz", false, false); + Rotation r_int = Rotation::from_euler(c, b, a, "zyx", true, false); + + EXPECT_TRUE(r_ext.equals(r_int, 1.E-9)); +} + +// ============================================================================= +// Gimbal Lock Tests +// ============================================================================= + +TEST(TRotationClass, gimbal_lock_zxz_theta_zero) +{ + // theta=0 for zxz: psi and phi are not individually recoverable + double psi = 0.5, theta = 0.0, phi = 0.3; + Rotation r = Rotation::from_euler(psi, theta, phi, "zxz", true, false); + + // Round-trip through euler should still give the same rotation + vec::fixed<3> angles = r.as_euler("zxz", true, false); + Rotation r2 = Rotation::from_euler(angles(0), angles(1), angles(2), "zxz", true, false); + + EXPECT_TRUE(r.equals(r2, 1.E-9)); +} + +TEST(TRotationClass, gimbal_lock_zxz_theta_pi) +{ + // theta=pi for zxz: another singular configuration + double psi = 0.5, theta = simcoon::pi, phi = 0.3; + Rotation r = Rotation::from_euler(psi, theta, phi, "zxz", true, false); + + // The rotation itself must be valid + mat::fixed<3,3> R = r.as_matrix(); + EXPECT_LT(norm(R * R.t() - eye(3,3), "fro"), 1.E-9); + EXPECT_NEAR(det(R), 1.0, 1.E-9); + + // Round-trip should preserve the rotation + vec::fixed<3> angles = r.as_euler("zxz", true, false); + Rotation r2 = Rotation::from_euler(angles(0), angles(1), angles(2), "zxz", true, false); + EXPECT_TRUE(r.equals(r2, 1.E-9)); +} + +TEST(TRotationClass, gimbal_lock_xyz_theta_pi_half) +{ + // theta=pi/2 for xyz (Tait-Bryan): gimbal lock + double psi = 0.4, theta = simcoon::pi / 2.0, phi = 0.6; + Rotation r = Rotation::from_euler(psi, theta, phi, "xyz", true, false); + + mat::fixed<3,3> R = r.as_matrix(); + EXPECT_LT(norm(R * R.t() - eye(3,3), "fro"), 1.E-9); + + vec::fixed<3> angles = r.as_euler("xyz", true, false); + Rotation r2 = Rotation::from_euler(angles(0), angles(1), angles(2), "xyz", true, false); + EXPECT_TRUE(r.equals(r2, 1.E-9)); +} + +// ============================================================================= +// Edge Cases +// ============================================================================= + +TEST(TRotationClass, small_angle) +{ + // Very small rotation should not cause numerical issues + double small_angle = 1.E-10; + Rotation r = Rotation::from_axis_angle(small_angle, 3); + + vec::fixed<3> v = {1, 0, 0}; + vec::fixed<3> v_rot = r.apply(v); + + // Should be nearly unchanged + EXPECT_LT(norm(v - v_rot), 1.E-8); +} + +TEST(TRotationClass, rotation_180_degrees) +{ + // 180 degree rotation around z-axis + Rotation r = Rotation::from_axis_angle(simcoon::pi, 3); + + vec::fixed<3> v = {1, 0, 0}; + vec::fixed<3> v_rot = r.apply(v); + + // x should become -x + EXPECT_NEAR(v_rot(0), -1.0, 1.E-9); + EXPECT_NEAR(v_rot(1), 0.0, 1.E-9); + EXPECT_NEAR(v_rot(2), 0.0, 1.E-9); +} + +TEST(TRotationClass, random_rotation_properties) +{ + // Random rotation should produce a valid rotation matrix + Rotation r = Rotation::random(); + mat::fixed<3,3> R = r.as_matrix(); + + // R * R^T should be identity + EXPECT_LT(norm(R * R.t() - eye(3,3), "fro"), 1.E-9); + + // det(R) should be 1 + EXPECT_NEAR(det(R), 1.0, 1.E-9); +} + +// ============================================================================= +// Matrix Round-trip +// ============================================================================= + +TEST(TRotationClass, matrix_roundtrip) +{ + // Create rotation, convert to matrix, convert back + Rotation r1 = Rotation::from_euler(0.5, 0.3, 0.7, "zxz"); + mat::fixed<3,3> R = r1.as_matrix(); + Rotation r2 = Rotation::from_matrix(R); + + EXPECT_TRUE(r1.equals(r2, 1.E-9)); +} + +TEST(TRotationClass, quat_matrix_quat_roundtrip) +{ + // quat -> matrix -> quat round-trip + vec::fixed<4> q1 = {0.5, 0.5, 0.5, 0.5}; // This is already unit quaternion + q1 /= norm(q1); // Ensure normalization + + Rotation r1 = Rotation::from_quat(q1); + mat::fixed<3,3> R = r1.as_matrix(); + Rotation r2 = Rotation::from_matrix(R); + + // Note: q and -q represent the same rotation + EXPECT_TRUE(r1.equals(r2, 1.E-9)); +}