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