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