EquatorialCompression.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 on the sphere.
25  * \image html EquatorialCompression.png "A sphere with an `aspect_ratio` of 3."
26  *
27  * \details A mapping from the sphere to itself which depends on a single
28  * parameter, the `aspect_ratio` \f$\alpha\f$, which is the ratio of the
29  * horizontal length to the vertical height for a given point. This parameter
30  * name was chosen because points with \f$\tan \theta = 1\f$ get mapped to
31  * points with \f$\tan \theta' = \alpha\f$. In general, gridpoints located
32  * at an angle \f$\theta\f$ from the pole are mapped to a new angle
33  * \f$\theta'\f$ satisfying \f$\tan \theta' = \alpha \tan \theta\f$.
34  *
35  * For an `aspect_ratio` greater than one, the gridpoints are mapped towards
36  * the equator, leading to an equatorially compressed grid. For an
37  * `aspect_ratio` less than one, the gridpoints are mapped towards the poles.
38  * Note that the aspect ratio must be positive.
39  *
40  * We define the auxiliary variables \f$ r := \sqrt{x^2 + y^2 +z^2}\f$
41  * and \f$ \rho := \sqrt{x^2 + y^2 + \alpha^{-2} z^2}\f$.
42  *
43  * The map corresponding to this transformation in cartesian coordinates
44  * is then given by:
45  *
46  * \f[\vec{x}'(x,y,z) =
47  * \frac{r}{\rho}\begin{bmatrix}
48  * x\\
49  * y\\
50  * \alpha^{-1} z\\
51  * \end{bmatrix}\f]
52  *
53  */
55  public:
56  static constexpr size_t dim = 3;
57  explicit EquatorialCompression(double aspect_ratio) noexcept;
58  EquatorialCompression() = default;
59  ~EquatorialCompression() = default;
62  EquatorialCompression& operator=(const EquatorialCompression&) = default;
63  EquatorialCompression& operator=(EquatorialCompression&&) = default;
64 
65  template <typename T>
67  const std::array<T, 3>& source_coords) const noexcept;
68 
69  boost::optional<std::array<double, 3>> inverse(
70  const std::array<double, 3>& target_coords) const noexcept;
71 
72  template <typename T>
73  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
74  const std::array<T, 3>& source_coords) const noexcept;
75 
76  template <typename T>
77  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
78  const std::array<T, 3>& source_coords) const noexcept;
79 
80  // clang-tidy: google runtime references
81  void pup(PUP::er& p) noexcept; // NOLINT
82 
83  bool is_identity() const noexcept { return is_identity_; }
84 
85  private:
86  template <typename T>
87  std::array<tt::remove_cvref_wrap_t<T>, 3> angular_distortion(
88  const std::array<T, 3>& coords, double inverse_alpha) const noexcept;
89  template <typename T>
90  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>
91  angular_distortion_jacobian(const std::array<T, 3>& coords,
92  double inverse_alpha) const noexcept;
93  friend bool operator==(const EquatorialCompression& lhs,
94  const EquatorialCompression& rhs) noexcept;
95 
96  double aspect_ratio_{std::numeric_limits<double>::signaling_NaN()};
97  double inverse_aspect_ratio_{std::numeric_limits<double>::signaling_NaN()};
98  bool is_identity_{false};
99 };
100 bool operator!=(const EquatorialCompression& lhs,
101  const EquatorialCompression& rhs) noexcept;
102 } // namespace CoordinateMaps
103 } // namespace domain
Redistributes gridpoints on the sphere.
Definition: EquatorialCompression.hpp:54
Definition: Strahlkorper.hpp:14
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.