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