Line data Source code
1 1 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : /// \file 5 : /// Defines the class Rotation. 6 : 7 : #pragma once 8 : 9 : #include <array> 10 : #include <cstddef> 11 : #include <limits> 12 : #include <optional> 13 : 14 : #include "DataStructures/Tensor/Tensor.hpp" 15 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp" 16 : 17 : /// \cond 18 : namespace PUP { 19 : class er; 20 : } // namespace PUP 21 : /// \endcond 22 : 23 : namespace domain { 24 : namespace CoordinateMaps { 25 : 26 : /// \cond HIDDEN_SYMBOLS 27 : template <size_t Dim> 28 : class Rotation; 29 : /// \endcond 30 : 31 : /*! 32 : * \ingroup CoordinateMapsGroup 33 : * \brief Spatial rotation in two dimensions. 34 : * 35 : * Let \f$(R,\Phi)\f$ be the polar coordinates associated with 36 : * \f$(\xi,\eta)\f$. 37 : * Let \f$(r,\phi)\f$ be the polar coordinates associated with \f$(x,y)\f$. 38 : * Applies the spatial rotation \f$\phi = \Phi + \alpha\f$. 39 : * 40 : * The formula for the mapping is: 41 : *\f{eqnarray*} 42 : x &=& \xi \cos \alpha - \eta \sin \alpha \\ 43 : y &=& \xi \sin \alpha + \eta \cos \alpha 44 : \f}. 45 : */ 46 : template <> 47 1 : class Rotation<2> { 48 : public: 49 0 : static constexpr size_t dim = 2; 50 : 51 : /// Constructor. 52 : /// 53 : /// \param rotation_angle the angle \f$\alpha\f$ (in radians). 54 1 : explicit Rotation(double rotation_angle); 55 0 : Rotation() = default; 56 0 : ~Rotation() = default; 57 0 : Rotation(const Rotation&) = default; 58 0 : Rotation& operator=(const Rotation&) = default; 59 0 : Rotation(Rotation&&) = default; // NOLINT 60 0 : Rotation& operator=(Rotation&&) = default; 61 : 62 : template <typename T> 63 0 : std::array<tt::remove_cvref_wrap_t<T>, 2> operator()( 64 : const std::array<T, 2>& source_coords) const; 65 : 66 : /// The inverse function is only callable with doubles because the inverse 67 : /// might fail if called for a point out of range, and it is unclear 68 : /// what should happen if the inverse were to succeed for some points in a 69 : /// DataVector but fail for other points. 70 1 : std::optional<std::array<double, 2>> inverse( 71 : const std::array<double, 2>& target_coords) const; 72 : 73 : template <typename T> 74 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 2, Frame::NoFrame> jacobian( 75 : const std::array<T, 2>& source_coords) const; 76 : 77 : template <typename T> 78 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 2, Frame::NoFrame> inv_jacobian( 79 : const std::array<T, 2>& source_coords) const; 80 : 81 : // NOLINTNEXTLINE(google-runtime-references) 82 0 : void pup(PUP::er& p); 83 : 84 0 : bool is_identity() const { return is_identity_; } 85 : 86 : private: 87 0 : friend bool operator==(const Rotation<2>& lhs, const Rotation<2>& rhs); 88 : 89 0 : double rotation_angle_{std::numeric_limits<double>::signaling_NaN()}; 90 0 : tnsr::ij<double, 2, Frame::Grid> rotation_matrix_{ 91 : std::numeric_limits<double>::signaling_NaN()}; 92 0 : bool is_identity_{false}; 93 : }; 94 : 95 0 : bool operator!=(const Rotation<2>& lhs, const Rotation<2>& rhs); 96 : 97 : /*! 98 : * \ingroup CoordinateMapsGroup 99 : * \brief Spatial rotation in three dimensions using Euler angles 100 : * 101 : * Rotation angles should be specified in degrees. 102 : * First rotation \f$\alpha\f$ is about z axis. 103 : * Second rotation \f$\beta\f$ is about rotated y axis. 104 : * Third rotation \f$\gamma\f$ is about rotated z axis. 105 : * These rotations are of the \f$(\xi,\eta,\zeta)\f$ coordinate system with 106 : * respect 107 : * to the grid coordinates \f$(x,y,z)\f$. 108 : * 109 : * The formula for the mapping is: 110 : * \f{eqnarray*} 111 : * x &=& \xi (\cos\gamma \cos\beta \cos\alpha - \sin\gamma \sin\alpha) 112 : * + \eta (-\sin\gamma \cos\beta \cos\alpha - \cos\gamma \sin\alpha) 113 : * + \zeta \sin\beta \cos\alpha \\ 114 : * y &=& \xi (\cos\gamma \cos\beta \sin\alpha + \sin\gamma \cos\alpha) 115 : * + \eta (-\sin\gamma \cos\beta \sin\alpha + \cos\gamma \cos\alpha) 116 : * + \zeta \sin\beta \sin\alpha \\ 117 : * z &=& -\xi \cos\gamma \sin\beta + \eta \sin\gamma \sin\beta 118 : * + \zeta \cos\beta 119 : * \f} 120 : */ 121 : template <> 122 1 : class Rotation<3> { 123 : public: 124 0 : static constexpr size_t dim = 3; 125 : 126 : /// Constructor. 127 : /// 128 : /// \param rotation_about_z the angle \f$\alpha\f$ (in radians). 129 : /// \param rotation_about_rotated_y the angle \f$\beta\f$ (in radians). 130 : /// \param rotation_about_rotated_z the angle \f$\gamma\f$ (in radians). 131 1 : Rotation(double rotation_about_z, double rotation_about_rotated_y, 132 : double rotation_about_rotated_z); 133 0 : Rotation() = default; 134 0 : ~Rotation() = default; 135 0 : Rotation(const Rotation&) = default; 136 0 : Rotation& operator=(const Rotation&) = default; 137 0 : Rotation(Rotation&&) = default; // NOLINT 138 0 : Rotation& operator=(Rotation&&) = default; 139 : 140 : template <typename T> 141 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()( 142 : const std::array<T, 3>& source_coords) const; 143 : 144 0 : std::optional<std::array<double, 3>> inverse( 145 : const std::array<double, 3>& target_coords) const; 146 : 147 : template <typename T> 148 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian( 149 : const std::array<T, 3>& source_coords) const; 150 : 151 : template <typename T> 152 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian( 153 : const std::array<T, 3>& source_coords) const; 154 : 155 : // NOLINTNEXTLINE(google-runtime-references) 156 0 : void pup(PUP::er& p); 157 : 158 0 : bool is_identity() const { return is_identity_; } 159 : 160 : private: 161 0 : friend bool operator==(const Rotation<3>& lhs, const Rotation<3>& rhs); 162 : 163 0 : double rotation_about_z_{std::numeric_limits<double>::signaling_NaN()}; 164 0 : double rotation_about_rotated_y_{ 165 : std::numeric_limits<double>::signaling_NaN()}; 166 0 : double rotation_about_rotated_z_{ 167 : std::numeric_limits<double>::signaling_NaN()}; 168 0 : tnsr::ij<double, 3, Frame::Grid> rotation_matrix_{ 169 : std::numeric_limits<double>::signaling_NaN()}; 170 0 : bool is_identity_{false}; 171 : }; 172 : 173 0 : bool operator!=(const Rotation<3>& lhs, const Rotation<3>& rhs); 174 : 175 : } // namespace CoordinateMaps 176 : } // namespace domain