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 <limits>
       9             : #include <optional>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
      13             : 
      14             : /// \cond
      15             : namespace PUP {
      16             : class er;
      17             : }  // namespace PUP
      18             : /// \endcond
      19             : 
      20             : namespace domain {
      21             : namespace CoordinateMaps {
      22             : 
      23             : /*!
      24             :  * \ingroup CoordinateMapsGroup
      25             :  *
      26             :  * \brief Redistributes gridpoints within the unit sphere.
      27             :  * \image html SpecialMobius.png "A sphere with a `mu` of 0.25."
      28             :  *
      29             :  * \details A special case of the conformal Mobius transformation that
      30             :  * maps the unit ball to itself. This map depends on a single
      31             :  * parameter, `mu` \f$ = \mu\f$, which is the x-coordinate of the preimage
      32             :  * of the origin under this map. This map has the fixed points \f$x=1\f$ and
      33             :  * \f$x=-1\f$. The map is singular for \f$\mu=1\f$ but we have found that this
      34             :  * map is accurate up to 12 decimal places for values of \f$\mu\f$ up to 0.96.
      35             :  *
      36             :  * We define the auxiliary variables
      37             :  * \f[ r := \sqrt{x^2 + y^2 +z^2}\f]
      38             :  * and
      39             :  * \f[ \lambda := \frac{1}{1 - 2 x \mu + \mu^2 r^2}\f]
      40             :  *
      41             :  * The map corresponding to this transformation in cartesian coordinates
      42             :  * is then given by:
      43             :  *
      44             :  * \f[\vec{x}'(x,y,z) =
      45             :  * \lambda\begin{bmatrix}
      46             :  * x(1+\mu^2) - \mu(1+r^2)\\
      47             :  * y(1-\mu^2)\\
      48             :  * z(1-\mu^2)\\
      49             :  * \end{bmatrix}\f]
      50             :  *
      51             :  * The inverse map is the same as the forward map with \f$\mu\f$
      52             :  * replaced by \f$-\mu\f$.
      53             :  *
      54             :  * This map is intended to be used only inside the unit sphere.  A
      55             :  * point inside the unit sphere maps to another point inside the unit
      56             :  * sphere. The map can have undesirable behavior at certain points
      57             :  * outside the unit sphere: The map is singular at
      58             :  * \f$(x,y,z) = (1/\mu, 0, 0)\f$ (which is outside the unit sphere
      59             :  * since \f$|\mu| < 1\f$). Moreover, a point on the \f$x\f$-axis
      60             :  * arbitrarily close to the singularity maps to an arbitrarily large
      61             :  * value on the \f$\pm x\f$-axis, where the sign depends on which side
      62             :  * of the singularity the point is on.
      63             :  *
      64             :  * A general Mobius transformation is a function on the complex plane, and
      65             :  * takes the form \f$ f(z) = \frac{az+b}{cz+d}\f$, where
      66             :  * \f$z, a, b, c, d \in \mathbb{C}\f$, and \f$ad-bc\neq 0\f$.
      67             :  *
      68             :  * The special case used in this map is the function
      69             :  * \f$ f(z) = \frac{z - \mu}{1 - z\mu}\f$. This has the desired properties:
      70             :  * - The unit disk in the complex plane is mapped to itself.
      71             :  *
      72             :  * - The x-axis is mapped to itself.
      73             :  *
      74             :  * - \f$f(\mu) = 0\f$.
      75             :  *
      76             :  * The three-dimensional version of this map is obtained by rotating the disk
      77             :  * in the plane about the x-axis.
      78             :  *
      79             :  * This map is useful for performing transformations along the x-axis
      80             :  * that preserve the unit disk. A concrete example of this is in the BBH
      81             :  * domain, where two BBHs with a center-of-mass at x=\f$\mu\f$ can be shifted
      82             :  * such that the new center of mass is now located at x=0. Additionally,
      83             :  * the spherical shape of the outer wave-zone is preserved and, as a mobius
      84             :  * map, the spherical coordinate shapes of the black holes is also preserved.
      85             :  */
      86           1 : class SpecialMobius {
      87             :  public:
      88           0 :   static constexpr size_t dim = 3;
      89           0 :   explicit SpecialMobius(double mu);
      90           0 :   SpecialMobius() = default;
      91           0 :   ~SpecialMobius() = default;
      92           0 :   SpecialMobius(SpecialMobius&&) = default;
      93           0 :   SpecialMobius(const SpecialMobius&) = default;
      94           0 :   SpecialMobius& operator=(const SpecialMobius&) = default;
      95           0 :   SpecialMobius& operator=(SpecialMobius&&) = default;
      96             : 
      97             :   template <typename T>
      98           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
      99             :       const std::array<T, 3>& source_coords) const;
     100             : 
     101             :   /// Returns std::nullopt for target_coords outside the unit sphere.
     102             :   /// The inverse function is only callable with doubles because the inverse
     103             :   /// might fail if called for a point out of range, and it is unclear
     104             :   /// what should happen if the inverse were to succeed for some points in a
     105             :   /// DataVector but fail for other points.
     106           1 :   std::optional<std::array<double, 3>> inverse(
     107             :       const std::array<double, 3>& target_coords) const;
     108             : 
     109             :   template <typename T>
     110           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
     111             :       const std::array<T, 3>& source_coords) const;
     112             : 
     113             :   template <typename T>
     114           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
     115             :       const std::array<T, 3>& source_coords) const;
     116             : 
     117             :   // NOLINTNEXTLINE(google-runtime-references)
     118           0 :   void pup(PUP::er& p);
     119             : 
     120           0 :   bool is_identity() const { return is_identity_; }
     121             : 
     122             :  private:
     123             :   template <typename T>
     124           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> mobius_distortion(
     125             :       const std::array<T, 3>& coords, double mu) const;
     126             :   template <typename T>
     127             :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>
     128           0 :   mobius_distortion_jacobian(const std::array<T, 3>& coords, double mu) const;
     129           0 :   friend bool operator==(const SpecialMobius& lhs, const SpecialMobius& rhs);
     130             : 
     131           0 :   double mu_{std::numeric_limits<double>::signaling_NaN()};
     132           0 :   bool is_identity_{false};
     133             : };
     134           0 : bool operator!=(const SpecialMobius& lhs, const SpecialMobius& rhs);
     135             : }  // namespace CoordinateMaps
     136             : }  // namespace domain