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