Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <memory>
9 : #include <optional>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <unordered_set>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 :
16 : /// \cond
17 : namespace domain {
18 : namespace FunctionsOfTime {
19 : class FunctionOfTime;
20 : } // namespace FunctionsOfTime
21 : } // namespace domain
22 : namespace PUP {
23 : class er;
24 : } // namespace PUP
25 : /// \endcond
26 :
27 : namespace domain {
28 : namespace CoordinateMaps {
29 : namespace TimeDependent {
30 :
31 : /*!
32 : * \ingroup CoordMapsTimeDependentGroup
33 : * \brief Time-dependent spatial rotation in two or three dimensions.
34 : *
35 : * ### General Transformation
36 : *
37 : * Let the source coordinates \f$ \vec{\xi} \f$ be mapped to coordinates \f$
38 : * \vec{x} \f$ using the transformation
39 : *
40 : * \f[ \vec{x} = R(t)\vec{\xi}, \f]
41 : *
42 : * where \f$ R(t) \f$ is a rotation matrix of proper dimensionality (defined
43 : * below) and \f$ A\vec{v} \f$ is the standard matrix-vector multiplicaton. For
44 : * 2D rotation, \f$ \vec{\xi} = \left(\xi, \eta\right) \f$ and \f$ \vec{x} =
45 : * \left(x, y\right) \f$ while for 3D rotations \f$ \vec{\xi} = \left(\xi, \eta,
46 : * \zeta\right) \f$ and \f$ \vec{x} = \left(x, y, z\right) \f$.
47 : *
48 : * The inverse transformation is
49 : *
50 : * \f[ \vec{\xi} = R^T(t) \vec{x} \f]
51 : *
52 : * because the inverse of a rotation matrix is its transpose.
53 : *
54 : * The frame velocity \f$ \vec{v} = d\vec{x}/dt \f$ is
55 : *
56 : * \f[ \vec{v} = \frac{d}{dt}\big( R(t) \big) \vec{\xi} \f]
57 : *
58 : * where \f$ d(R(t))/dt \f$ is the time derivative of the rotation matrix.
59 : *
60 : * The components of the Jacobian \f$ \partial x^i/\partial\xi^j \f$ are
61 : * trivially related to the components of the rotation matrix by
62 : *
63 : * \f[ \partial x^i/\partial\xi^j = R_{ij}, \f]
64 : *
65 : * and similarly the components of the inverse Jacobian \f$ \partial
66 : * \xi^i/\partial x^j \f$ are
67 : *
68 : * \f[ \partial \xi^i/\partial x^j = R^{-1}_{ij} = R^T_{ij} = R_{ji}. \f]
69 : *
70 : * ### 2D Rotation Matrix
71 : *
72 : * The 2D rotaion matrix is defined in the usual way as
73 : *
74 : * \f[
75 : * R(t) =
76 : * \begin{bmatrix}
77 : * \cos(\theta(t)) & -\sin(\theta(t)) \\
78 : * \sin(\theta(t)) & \cos(\theta(t)) \\
79 : * \end{bmatrix}.
80 : * \f]
81 : *
82 : * We associate the polar coordinates \f$ \left( \mathrm{P}, \Phi\right) \f$
83 : * with the unmapped coordinates \f$ \left(\xi, \eta\right) \f$ and the polar
84 : * coordinates \f$ \left(r,\phi\right) \f$ with the mapped coordinates \f$
85 : * \left(x, y\right) \f$. We then have \f$ \phi = \Phi + \theta(t) \f$.
86 : *
87 : * The derivative of the rotation matrix is then
88 : *
89 : * \f[
90 : * R(t) =
91 : * \begin{bmatrix}
92 : * -\omega(t) \sin(\theta(t)) & -\omega(t)\cos(\theta(t)) \\
93 : * \omega(t) \cos(\theta(t)) & -\omega(t)\sin(\theta(t)) \\
94 : * \end{bmatrix}.
95 : * \f]
96 : *
97 : * where \f$ \omega(t) = d\theta(t)/dt \f$.
98 : *
99 : * \note This 2D rotation is assumed to be in the \f$ xy \f$-plane (about the
100 : * \f$ z \f$-axis).
101 : *
102 : * ### 3D Rotation Matrix
103 : *
104 : * For 3D rotations, we use quaternions to represent rotations about an
105 : * arbitrary axis. We define a unit quaternion as
106 : *
107 : * \f[
108 : * \mathbf{q}
109 : * = \left(q_0, q_1, q_2, q_3\right)
110 : * = \left(q_0, \vec{q}\right)
111 : * = \left(\cos(\frac{\theta(t)}{2}),
112 : * \hat{n}\sin(\frac{\theta(t)}{2})\right)
113 : * \f]
114 : *
115 : * where \f$ \hat{n} \f$ is our arbitrary rotation axis and \f$ \theta(t) \f$ is
116 : * the angle rotated about that axis. A rotation in 3D is then defined as
117 : *
118 : * \f[ \mathbf{x} = \mathbf{q}\mathbf{\xi}\mathbf{q}^* \f]
119 : *
120 : * where \f$ \mathbf{q}^* = \left(\cos(\theta(t)/2),
121 : * -\hat{n}\sin(\theta(t)/2)\right) \f$ and we promote the vectors to
122 : * quaternions as \f$ \mathbf{x} = \left(0, \vec{x}\right) \f$ and \f$
123 : * \mathbf{\xi} = \left(0, \vec{\xi}\right) \f$. This will rotate the vector \f$
124 : * \vec{\xi} \f$ about \f$ \hat{n} \f$ by an angle \f$ \theta(t) \f$,
125 : * transforming it into \f$ \vec{x} \f$.
126 : *
127 : * We can represent this rotation using quaternions as a rotation matrix of
128 : * the form
129 : *
130 : * \f[
131 : * R(t) =
132 : * \begin{bmatrix}
133 : * q_0^2 + q_1^2 - q_2^2 - q_3^2 & 2(q_1q_2 - q_0q_3) & 2(q_1q_3 + q_0q_2) \\
134 : * 2(q_1q_2 + q_0q_3) & q_0^2 + q_2^2 - q_1^2 - q_3^2 & 2(q_2q_3 - q_0q_1) \\
135 : * 2(q_1q_3 - q_0q_2) & 2(q_2q_3 + q_0q_1) & q_0^2 + q_3^2 - q_1^2 - q_2^2 \\
136 : * \end{bmatrix}.
137 : * \f]
138 : *
139 : * The derivative of this rotation matrix can expressed in a similar form
140 : *
141 : * \f[
142 : * R(t) =
143 : * \begin{bmatrix}
144 : * 2(q_0\dot{q_0} + q_1\dot{q_1} - q_2\dot{q_2} - q_3\dot{q_3})
145 : * & 2(\dot{q_1}q_2 + q_1\dot{q_2} - \dot{q_0}q_3 - q_0\dot{q_3})
146 : * & 2(\dot{q_1}q_3 + q_1\dot{q_3} + \dot{q_0}q_2 + q_0\dot{q_2}) \\
147 : * 2(\dot{q_1}q_2 + q_1\dot{q_2} + \dot{q_0}q_3 + q_0\dot{q_3})
148 : * & 2(q_0\dot{q_0} + q_2\dot{q_2} - q_1\dot{q_1} - q_3\dot{q_3})
149 : * & 2(\dot{q_2}q_3 + q_2\dot{q_3} - \dot{q_0}q_1 - q_0\dot{q_1}) \\
150 : * 2(\dot{q_1}q_3 + q_1\dot{q_3} - \dot{q_0}q_2 - q_0\dot{q_2})
151 : * & 2(\dot{q_2}q_3 + q_2\dot{q_3} + \dot{q_0}q_1 + q_0\dot{q_1})
152 : * & 2(q_0\dot{q_0} + q_3\dot{q_3} - q_1\dot{q_1} - q_2\dot{q_2}) \\
153 : * \end{bmatrix}.
154 : * \f]
155 : *
156 : * \note If you choose \f$ \hat{n} = (0, 0, 1) \f$, this rotation will be
157 : * equivalent to the 2D rotation.
158 : */
159 : template <size_t Dim>
160 1 : class Rotation {
161 : public:
162 : static_assert(Dim == 2 or Dim == 3,
163 : "Rotation map can only be constructed in 2 or 3 dimensions.");
164 0 : static constexpr size_t dim = Dim;
165 :
166 0 : explicit Rotation(std::string function_of_time_name);
167 0 : Rotation() = default;
168 :
169 : template <typename T>
170 0 : std::array<T, Dim> operator()(
171 : const std::array<T, Dim>& source_coords, double time,
172 : const std::unordered_map<
173 : std::string,
174 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
175 : functions_of_time) const;
176 :
177 : /// The inverse function is only callable with doubles because the inverse
178 : /// might fail if called for a point out of range, and it is unclear
179 : /// what should happen if the inverse were to succeed for some points in a
180 : /// DataVector but fail for other points.
181 1 : std::optional<std::array<double, Dim>> inverse(
182 : const std::array<double, Dim>& target_coords, double time,
183 : const std::unordered_map<
184 : std::string,
185 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
186 : functions_of_time) const;
187 :
188 : template <typename T>
189 0 : std::array<T, Dim> frame_velocity(
190 : const std::array<T, Dim>& source_coords, double time,
191 : const std::unordered_map<
192 : std::string,
193 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
194 : functions_of_time) const;
195 :
196 : template <typename T>
197 0 : tnsr::Ij<T, Dim, Frame::NoFrame> jacobian(
198 : const std::array<T, Dim>& source_coords, double time,
199 : const std::unordered_map<
200 : std::string,
201 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
202 : functions_of_time) const;
203 :
204 : template <typename T>
205 0 : tnsr::Ij<T, Dim, Frame::NoFrame> inv_jacobian(
206 : const std::array<T, Dim>& source_coords, double time,
207 : const std::unordered_map<
208 : std::string,
209 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
210 : functions_of_time) const;
211 :
212 : // NOLINTNEXTLINE(google-runtime-references)
213 0 : void pup(PUP::er& p);
214 :
215 0 : static bool is_identity() { return false; }
216 :
217 0 : static constexpr bool supports_hessian{true};
218 :
219 0 : const std::unordered_set<std::string>& function_of_time_names() const {
220 : return f_of_t_names_;
221 : }
222 :
223 : private:
224 : template <size_t LocalDim>
225 : // NOLINTNEXTLINE(readability-redundant-declaration)
226 0 : friend bool operator==(const Rotation<LocalDim>& lhs,
227 : const Rotation<LocalDim>& rhs);
228 0 : std::string f_of_t_name_;
229 0 : std::unordered_set<std::string> f_of_t_names_;
230 : };
231 :
232 : template <size_t Dim>
233 0 : bool operator!=(const Rotation<Dim>& lhs, const Rotation<Dim>& rhs);
234 :
235 : } // namespace TimeDependent
236 : } // namespace CoordinateMaps
237 : } // namespace domain
|
