Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines the class Equiangular.
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <array>
      10             : #include <cmath>
      11             : #include <cstddef>
      12             : #include <optional>
      13             : 
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
      16             : 
      17             : /// \cond
      18             : namespace PUP {
      19             : class er;
      20             : }  // namespace PUP
      21             : /// \endcond
      22             : 
      23             : namespace domain {
      24             : namespace CoordinateMaps {
      25             : 
      26             : /*!
      27             :  * \ingroup CoordinateMapsGroup
      28             :  * \brief Non-linear map from \f$\xi \in [A, B]\rightarrow x \in [a, b]\f$.
      29             :  *
      30             :  * The formula for the mapping is:
      31             :  * \f{align}
      32             :  * x &= \frac{a}{2} \left(1-\mathrm{tan}\left(
      33             :  *      \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right) +
      34             :  *      \frac{b}{2} \left(1+\mathrm{tan}\left(
      35             :  *      \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right)\\
      36             :  * \xi &= \frac{A}{2} \left(1-\frac{4}{\pi}\mathrm{arctan}\left(
      37             :  *        \frac{2x-a-b}{b-a}\right)\right)+
      38             :  *        \frac{B}{2} \left(1+\frac{4}{\pi}\mathrm{arctan}\left(
      39             :  *        \frac{2x-a-b}{b-a}\right)\right)
      40             :  * \f}
      41             :  *
      42             :  * \note The intermediate step in which a tangent map is applied can be more
      43             :  * clearly understood if we define the coordinates:
      44             :  * \f{align}
      45             :  * \xi_{logical} &:= \frac{2\xi-B-A}{B-A} \in [-1, 1]\\
      46             :  * \Xi &:= \mathrm{tan}\left(\frac{\pi\xi_{logical}}{4}\right) \in [-1, 1]
      47             :  * \f}
      48             :  *
      49             :  * This map is intended to be used with the `Wedge` map when equiangular
      50             :  * coordinates are chosen for those maps. For more information on this choice
      51             :  * of coordinates, see the documentation for `Wedge`.
      52             :  */
      53           1 : class Equiangular {
      54             :  public:
      55           0 :   static constexpr size_t dim = 1;
      56             : 
      57           0 :   Equiangular(double A, double B, double a, double b);
      58             : 
      59           0 :   Equiangular() = default;
      60           0 :   ~Equiangular() = default;
      61           0 :   Equiangular(const Equiangular&) = default;
      62           0 :   Equiangular(Equiangular&&) = default;
      63           0 :   Equiangular& operator=(const Equiangular&) = default;
      64           0 :   Equiangular& operator=(Equiangular&&) = default;
      65             : 
      66             :   template <typename T>
      67           0 :   std::array<tt::remove_cvref_wrap_t<T>, 1> operator()(
      68             :       const std::array<T, 1>& source_coords) const;
      69             : 
      70             :   /// The inverse function is only callable with doubles because the inverse
      71             :   /// might fail if called for a point out of range, and it is unclear
      72             :   /// what should happen if the inverse were to succeed for some points in a
      73             :   /// DataVector but fail for other points.
      74           1 :   std::optional<std::array<double, 1>> inverse(
      75             :       const std::array<double, 1>& target_coords) const;
      76             : 
      77             :   template <typename T>
      78           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 1, Frame::NoFrame> jacobian(
      79             :       const std::array<T, 1>& source_coords) const;
      80             : 
      81             :   template <typename T>
      82           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 1, Frame::NoFrame> inv_jacobian(
      83             :       const std::array<T, 1>& source_coords) const;
      84             : 
      85             :   // NOLINTNEXTLINE(google-runtime-references)
      86           0 :   void pup(PUP::er& p);
      87             : 
      88           0 :   static bool is_identity() { return false; }
      89             : 
      90             :  private:
      91           0 :   friend bool operator==(const Equiangular& lhs, const Equiangular& rhs);
      92             : 
      93           0 :   double A_{-1.0};
      94           0 :   double B_{1.0};
      95           0 :   double a_{-1.0};
      96           0 :   double b_{1.0};
      97           0 :   double length_of_domain_over_m_pi_4_{(B_ - A_) / M_PI_4};  // 4(B-A)/\pi
      98           0 :   double length_of_range_{2.0};                              // b-a
      99           0 :   double m_pi_4_over_length_of_domain_{M_PI_4 / (B_ - A_)};
     100           0 :   double one_over_length_of_range_{0.5};
     101             :   // The jacobian for the affine map with the same parameters.
     102           0 :   double linear_jacobian_times_m_pi_4_{length_of_range_ /
     103             :                                        length_of_domain_over_m_pi_4_};
     104             :   // The inverse jacobian for the affine map with the same parameters.
     105           0 :   double linear_inverse_jacobian_over_m_pi_4_{length_of_domain_over_m_pi_4_ /
     106             :                                               length_of_range_};
     107             : };
     108             : 
     109           0 : inline bool operator!=(const CoordinateMaps::Equiangular& lhs,
     110             :                        const CoordinateMaps::Equiangular& rhs) {
     111             :   return not(lhs == rhs);
     112             : }
     113             : 
     114             : }  // namespace CoordinateMaps
     115             : }  // namespace domain