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 :
16 : /// \cond
17 : namespace PUP {
18 : class er;
19 : } // namespace PUP
20 : /// \endcond
21 :
22 : namespace domain {
23 : namespace CoordinateMaps {
24 :
25 : /// \cond HIDDEN_SYMBOLS
26 : template <size_t Dim>
27 : class Rotation;
28 : /// \endcond
29 :
30 : /*!
31 : * \ingroup CoordinateMapsGroup
32 : * \brief Spatial rotation in two dimensions.
33 : *
34 : * Let \f$(R,\Phi)\f$ be the polar coordinates associated with
35 : * \f$(\xi,\eta)\f$.
36 : * Let \f$(r,\phi)\f$ be the polar coordinates associated with \f$(x,y)\f$.
37 : * Applies the spatial rotation \f$\phi = \Phi + \alpha\f$.
38 : *
39 : * The formula for the mapping is:
40 : *\f{eqnarray*}
41 : x &=& \xi \cos \alpha - \eta \sin \alpha \\
42 : y &=& \xi \sin \alpha + \eta \cos \alpha
43 : \f}.
44 : */
45 : template <>
46 1 : class Rotation<2> {
47 : public:
48 0 : static constexpr size_t dim = 2;
49 :
50 : /// Constructor.
51 : ///
52 : /// \param rotation_angle the angle \f$\alpha\f$ (in radians).
53 1 : explicit Rotation(double rotation_angle);
54 0 : Rotation() = default;
55 0 : ~Rotation() = default;
56 0 : Rotation(const Rotation&) = default;
57 0 : Rotation& operator=(const Rotation&) = default;
58 0 : Rotation(Rotation&&) = default; // NOLINT
59 0 : Rotation& operator=(Rotation&&) = default;
60 :
61 : template <typename T>
62 0 : std::array<T, 2> operator()(const std::array<T, 2>& source_coords) const;
63 :
64 : /// The inverse function is only callable with doubles because the inverse
65 : /// might fail if called for a point out of range, and it is unclear
66 : /// what should happen if the inverse were to succeed for some points in a
67 : /// DataVector but fail for other points.
68 1 : std::optional<std::array<double, 2>> inverse(
69 : const std::array<double, 2>& target_coords) const;
70 :
71 : template <typename T>
72 0 : tnsr::Ij<T, 2, Frame::NoFrame> jacobian(
73 : const std::array<T, 2>& source_coords) const;
74 :
75 : template <typename T>
76 0 : tnsr::Ij<T, 2, Frame::NoFrame> inv_jacobian(
77 : const std::array<T, 2>& source_coords) const;
78 :
79 : // NOLINTNEXTLINE(google-runtime-references)
80 0 : void pup(PUP::er& p);
81 :
82 0 : bool is_identity() const { return is_identity_; }
83 :
84 0 : static constexpr bool supports_hessian{true};
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<T, 3> operator()(const std::array<T, 3>& source_coords) const;
142 :
143 0 : std::optional<std::array<double, 3>> inverse(
144 : const std::array<double, 3>& target_coords) const;
145 :
146 : template <typename T>
147 0 : tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
148 : const std::array<T, 3>& source_coords) const;
149 :
150 : template <typename T>
151 0 : tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
152 : const std::array<T, 3>& source_coords) const;
153 :
154 : // NOLINTNEXTLINE(google-runtime-references)
155 0 : void pup(PUP::er& p);
156 :
157 0 : bool is_identity() const { return is_identity_; }
158 :
159 0 : static constexpr bool supports_hessian{true};
160 :
161 : private:
162 0 : friend bool operator==(const Rotation<3>& lhs, const Rotation<3>& rhs);
163 :
164 0 : double rotation_about_z_{std::numeric_limits<double>::signaling_NaN()};
165 0 : double rotation_about_rotated_y_{
166 : std::numeric_limits<double>::signaling_NaN()};
167 0 : double rotation_about_rotated_z_{
168 : std::numeric_limits<double>::signaling_NaN()};
169 0 : tnsr::ij<double, 3, Frame::Grid> rotation_matrix_{
170 : std::numeric_limits<double>::signaling_NaN()};
171 0 : bool is_identity_{false};
172 : };
173 :
174 0 : bool operator!=(const Rotation<3>& lhs, const Rotation<3>& rhs);
175 :
176 : } // namespace CoordinateMaps
177 : } // namespace domain
|
