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 <optional>
       9             : 
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Domain/CoordinateMaps/Distribution.hpp"
      12             : #include "Utilities/Serialization/PupStlCpp17.hpp"
      13             : 
      14             : /// \cond
      15             : namespace PUP {
      16             : class er;
      17             : }  // namespace PUP
      18             : /// \endcond
      19             : 
      20             : namespace domain::CoordinateMaps {
      21             : 
      22             : /*!
      23             :  * \ingroup CoordinateMapsGroup
      24             :  * \brief Maps \f$\xi\f$ in the 1D interval \f$[A, B]\f$ to \f$x\f$ in the
      25             :  * interval \f$[a, b]\f$ according to a `domain::CoordinateMaps::Distribution`.
      26             :  *
      27             :  * \details The mapping takes a `domain::CoordinateMaps::Distribution` and
      28             :  * distributes the grid points accordingly.
      29             :  *
      30             :  * The formula for the mapping is, in case of a `Linear` distribution
      31             :  * \f{align}
      32             :  * x &= \frac{b}{B-A} (\xi-A) + \frac{a}{B-A} (B-\xi)\\
      33             :  * \xi &=\frac{B}{b-a} (x-a) + \frac{A}{b-a} (b-x)
      34             :  * \f}
      35             :  *
      36             :  * For every other distribution we use this linear mapping to map onto an
      37             :  * interval \f$[-1, 1]\f$, so we define:
      38             :  * \f{align}
      39             :  * f(\xi) &:= \frac{A+B-2\xi}{A-B} \in [-1, 1]\\
      40             :  * g(x) &:= \frac{a+b-2x}{a-b} \in [-1, 1]
      41             :  * \f}
      42             :  *
      43             :  * With this an `Equiangular` distribution is described by
      44             :  *
      45             :  * \f{align}
      46             :  * x &= \frac{a}{2} \left(1-\mathrm{tan}\left(\frac{\pi}{4}f(\xi)\right)\right)
      47             :  * + \frac{b}{2} \left(1+\mathrm{tan}\left(\frac{\pi}{4}f(\xi)\right)\right)\\
      48             :  * \xi &=
      49             :  * \frac{A}{2} \left(1-\frac{4}{\pi}\mathrm{arctan}\left(g(x)\right)\right) +
      50             :  * \frac{B}{2} \left(1+\frac{4}{\pi}\mathrm{arctan}\left(g(x)\right)\right)
      51             :  * \f}
      52             :  *
      53             :  * \note The equiangular distribution is intended to be used with the `Wedge`
      54             :  * map when equiangular coordinates are chosen for those maps. For more
      55             :  * information on this choice of coordinates, see the documentation for `Wedge`.
      56             :  *
      57             :  *
      58             :  * For both the `Logarithmic` and `Inverse` distribution, we first specify a
      59             :  * position for the singularity \f$c:=\f$`singularity_pos` outside the target
      60             :  * interval \f$[a, b]\f$.
      61             :  *
      62             :  * The `Logarithmic` distribution further requires the introduction of a
      63             :  * variable \f$\sigma\f$ dependent on whether \f$c\f$ is left (\f$ \sigma =
      64             :  * 1\f$) or right (\f$ \sigma = -1\f$) of the target interval. With this, the
      65             :  * `Logarithmic` distribution is described by
      66             :  *
      67             :  * \f{align}
      68             :  * x &= \sigma\, {\rm exp}\left(\frac{\ln(b-c)+\ln(a-c)}{2} + f(\xi)
      69             :  * \frac{\ln(b-c)-\ln(a-c)}{2}\right) + c\\
      70             :  * \xi &= \frac{B-A}{2}\frac{2\ln(\sigma [x-c]) - \ln(b-c) - \ln(a-c)}{\ln(b-c)
      71             :  * - \ln(a-c)} + \frac{B+A}{2} \f}
      72             :  *
      73             :  * and the `Inverse` distribution by
      74             :  *
      75             :  * \f{align}
      76             :  * x &= \frac{2(a-c)(b-c)}{a+b-2c+(a-b)f(\xi)} + c\\
      77             :  * \xi &= -\frac{A-B}{2(a-b)} \left(\frac{2(a-c)(b-c)}{x-c} - (a-c) -
      78             :  * (b-c)\right) + \frac{A+B}{2}
      79             :  * \f}
      80             :  */
      81           1 : class Interval {
      82             :  public:
      83           0 :   static constexpr size_t dim = 1;
      84             : 
      85           0 :   Interval(double A, double B, double a, double b, Distribution distribution,
      86             :            std::optional<double> singularity_pos = std::nullopt);
      87             : 
      88           0 :   Interval() = default;
      89             : 
      90             :   template <typename T>
      91           0 :   std::array<T, 1> operator()(const std::array<T, 1>& source_coords) const;
      92             : 
      93           0 :   std::optional<std::array<double, 1>> inverse(
      94             :       const std::array<double, 1>& target_coords) const;
      95             : 
      96             :   template <typename T>
      97           0 :   tnsr::Ij<T, 1, Frame::NoFrame> jacobian(
      98             :       const std::array<T, 1>& source_coords) const;
      99             : 
     100             :   template <typename T>
     101           0 :   tnsr::Ij<T, 1, Frame::NoFrame> inv_jacobian(
     102             :       const std::array<T, 1>& source_coords) const;
     103             : 
     104             :   // NOLINTNEXTLINE(google-runtime-references)
     105           0 :   void pup(PUP::er& p);
     106             : 
     107           0 :   bool is_identity() const { return is_identity_; }
     108             : 
     109           0 :   static constexpr bool supports_hessian{true};
     110             : 
     111             :  private:
     112           0 :   friend bool operator==(const Interval& lhs, const Interval& rhs);
     113             : 
     114           0 :   double A_{std::numeric_limits<double>::signaling_NaN()};
     115           0 :   double B_{std::numeric_limits<double>::signaling_NaN()};
     116           0 :   double a_{std::numeric_limits<double>::signaling_NaN()};
     117           0 :   double b_{std::numeric_limits<double>::signaling_NaN()};
     118           0 :   Distribution distribution_{Distribution::Linear};
     119           0 :   std::optional<double> singularity_pos_{std::nullopt};
     120           0 :   bool is_identity_{false};
     121             : };
     122             : 
     123           0 : inline bool operator!=(const CoordinateMaps::Interval& lhs,
     124             :                        const CoordinateMaps::Interval& rhs) {
     125             :   return not(lhs == rhs);
     126             : }
     127             : 
     128             : }  // namespace domain::CoordinateMaps