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