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             : 
      16             : /// \cond
      17             : namespace PUP {
      18             : class er;
      19             : }  // namespace PUP
      20             : /// \endcond
      21             : 
      22             : namespace domain {
      23             : namespace CoordinateMaps {
      24             : 
      25             : /*!
      26             :  * \ingroup CoordinateMapsGroup
      27             :  * \brief Non-linear map from \f$\xi \in [A, B]\rightarrow x \in [a, b]\f$.
      28             :  *
      29             :  * The formula for the mapping is:
      30             :  * \f{align}
      31             :  * x &= \frac{a}{2} \left(1-\mathrm{tan}\left(
      32             :  *      \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right) +
      33             :  *      \frac{b}{2} \left(1+\mathrm{tan}\left(
      34             :  *      \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right)\\
      35             :  * \xi &= \frac{A}{2} \left(1-\frac{4}{\pi}\mathrm{arctan}\left(
      36             :  *        \frac{2x-a-b}{b-a}\right)\right)+
      37             :  *        \frac{B}{2} \left(1+\frac{4}{\pi}\mathrm{arctan}\left(
      38             :  *        \frac{2x-a-b}{b-a}\right)\right)
      39             :  * \f}
      40             :  *
      41             :  * \note The intermediate step in which a tangent map is applied can be more
      42             :  * clearly understood if we define the coordinates:
      43             :  * \f{align}
      44             :  * \xi_{logical} &:= \frac{2\xi-B-A}{B-A} \in [-1, 1]\\
      45             :  * \Xi &:= \mathrm{tan}\left(\frac{\pi\xi_{logical}}{4}\right) \in [-1, 1]
      46             :  * \f}
      47             :  *
      48             :  * This map is intended to be used with the `Wedge` map when equiangular
      49             :  * coordinates are chosen for those maps. For more information on this choice
      50             :  * of coordinates, see the documentation for `Wedge`.
      51             :  */
      52           1 : class Equiangular {
      53             :  public:
      54           0 :   static constexpr size_t dim = 1;
      55             : 
      56           0 :   Equiangular(double A, double B, double a, double b);
      57             : 
      58           0 :   Equiangular() = default;
      59           0 :   ~Equiangular() = default;
      60           0 :   Equiangular(const Equiangular&) = default;
      61           0 :   Equiangular(Equiangular&&) = default;
      62           0 :   Equiangular& operator=(const Equiangular&) = default;
      63           0 :   Equiangular& operator=(Equiangular&&) = default;
      64             : 
      65             :   template <typename T>
      66           0 :   std::array<T, 1> operator()(const std::array<T, 1>& source_coords) const;
      67             : 
      68             :   /// The inverse function is only callable with doubles because the inverse
      69             :   /// might fail if called for a point out of range, and it is unclear
      70             :   /// what should happen if the inverse were to succeed for some points in a
      71             :   /// DataVector but fail for other points.
      72           1 :   std::optional<std::array<double, 1>> inverse(
      73             :       const std::array<double, 1>& target_coords) const;
      74             : 
      75             :   template <typename T>
      76           0 :   tnsr::Ij<T, 1, Frame::NoFrame> jacobian(
      77             :       const std::array<T, 1>& source_coords) const;
      78             : 
      79             :   template <typename T>
      80           0 :   tnsr::Ij<T, 1, Frame::NoFrame> inv_jacobian(
      81             :       const std::array<T, 1>& source_coords) const;
      82             : 
      83             :   // NOLINTNEXTLINE(google-runtime-references)
      84           0 :   void pup(PUP::er& p);
      85             : 
      86           0 :   static bool is_identity() { return false; }
      87             : 
      88           0 :   static constexpr bool supports_hessian{true};
      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