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