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 on the sphere.
      27             :  * \image html EquatorialCompression.png "A sphere with an `aspect_ratio` of 3."
      28             :  *
      29             :  * \details A mapping from the sphere to itself which redistributes points
      30             :  * towards (or away from) a user-specifed axis, indicated by `index_pole_axis_`.
      31             :  * Once the axis is selected, the map is determined by a single parameter,
      32             :  * the `aspect_ratio` \f$\alpha\f$, which is the ratio of the distance
      33             :  * perpendicular to the polar axis to the distance along the polar axis for
      34             :  * a given point. This parameter name was chosen because points with
      35             :  * \f$\tan \theta = 1\f$ get mapped to points with \f$\tan \theta' = \alpha\f$.
      36             :  * In general, gridpoints located at an angle \f$\theta\f$ from the pole are
      37             :  * mapped to a new angle
      38             :  * \f$\theta'\f$ satisfying \f$\tan \theta' = \alpha \tan \theta\f$.
      39             :  *
      40             :  * For an `aspect_ratio` greater than one, the gridpoints are mapped towards
      41             :  * the equator, leading to an equatorially compressed grid. For an
      42             :  * `aspect_ratio` less than one, the gridpoints are mapped towards the poles.
      43             :  * Note that the aspect ratio must be positive.
      44             :  *
      45             :  * Suppose the polar axis were the z-axis, given by `index_pole_axis_ == 2`.
      46             :  * We can then define the auxiliary variables \f$ r := \sqrt{x^2 + y^2 +z^2}\f$
      47             :  * and \f$ \rho := \sqrt{x^2 + y^2 + \alpha^{-2} z^2}\f$.
      48             :  *
      49             :  * The map corresponding to this transformation in cartesian coordinates
      50             :  * is then given by:
      51             :  *
      52             :  * \f[\vec{x}'(x,y,z) =
      53             :  * \frac{r}{\rho}\begin{bmatrix}
      54             :  * x\\
      55             :  * y\\
      56             :  * \alpha^{-1} z\\
      57             :  * \end{bmatrix}.\f]
      58             :  *
      59             :  * The mappings for polar axes along the x and y axes are similarly obtained.
      60             :  */
      61           1 : class EquatorialCompression {
      62             :  public:
      63           0 :   static constexpr size_t dim = 3;
      64           0 :   explicit EquatorialCompression(double aspect_ratio,
      65             :                                  size_t index_pole_axis = 2);
      66           0 :   EquatorialCompression() = default;
      67           0 :   ~EquatorialCompression() = default;
      68           0 :   EquatorialCompression(EquatorialCompression&&) = default;
      69           0 :   EquatorialCompression(const EquatorialCompression&) = default;
      70           0 :   EquatorialCompression& operator=(const EquatorialCompression&) = default;
      71           0 :   EquatorialCompression& operator=(EquatorialCompression&&) = default;
      72             : 
      73             :   template <typename T>
      74           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
      75             :       const std::array<T, 3>& source_coords) const;
      76             : 
      77             :   /// The inverse function is only callable with doubles because the inverse
      78             :   /// might fail if called for a point out of range, and it is unclear
      79             :   /// what should happen if the inverse were to succeed for some points in a
      80             :   /// DataVector but fail for other points.
      81           1 :   std::optional<std::array<double, 3>> inverse(
      82             :       const std::array<double, 3>& target_coords) const;
      83             : 
      84             :   template <typename T>
      85           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
      86             :       const std::array<T, 3>& source_coords) const;
      87             : 
      88             :   template <typename T>
      89           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
      90             :       const std::array<T, 3>& source_coords) const;
      91             : 
      92             :   // NOLINTNEXTLINE(google-runtime-references)
      93           0 :   void pup(PUP::er& p);
      94             : 
      95           0 :   bool is_identity() const { return is_identity_; }
      96             : 
      97             :  private:
      98             :   template <typename T>
      99           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> angular_distortion(
     100             :       const std::array<T, 3>& coords, double inverse_alpha) const;
     101             :   template <typename T>
     102             :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>
     103           0 :   angular_distortion_jacobian(const std::array<T, 3>& coords,
     104             :                               double inverse_alpha) const;
     105           0 :   friend bool operator==(const EquatorialCompression& lhs,
     106             :                          const EquatorialCompression& rhs);
     107             : 
     108           0 :   double aspect_ratio_{std::numeric_limits<double>::signaling_NaN()};
     109           0 :   double inverse_aspect_ratio_{std::numeric_limits<double>::signaling_NaN()};
     110           0 :   bool is_identity_{false};
     111           0 :   size_t index_pole_axis_{};
     112             : };
     113           0 : bool operator!=(const EquatorialCompression& lhs,
     114             :                 const EquatorialCompression& rhs);
     115             : }  // namespace CoordinateMaps
     116             : }  // namespace domain