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