Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines the class Affine.
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <array>
      10             : #include <cstddef>
      11             : #include <optional>
      12             : 
      13             : #include "DataStructures/Tensor/TypeAliases.hpp"
      14             : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
      15             : 
      16             : /// \cond
      17             : namespace PUP {
      18             : class er;
      19             : }  // namespace PUP
      20             : /// \endcond
      21             : 
      22             : namespace domain {
      23           1 : namespace CoordinateMaps {
      24             : 
      25             : /*!
      26             :  * \ingroup CoordinateMapsGroup
      27             :  * \brief Affine map from \f$\xi \in [A, B]\rightarrow x \in [a, b]\f$.
      28             :  *
      29             :  * The formula for the mapping is...
      30             :  * \f[
      31             :  * x = \frac{b}{B-A} (\xi-A) +\frac{a}{B-A}(B-\xi)
      32             :  * \f]
      33             :  * \f[
      34             :  * \xi =\frac{B}{b-a} (x-a) +\frac{A}{b-a}(b-x)
      35             :  * \f]
      36             :  */
      37           1 : class Affine {
      38             :  public:
      39           0 :   static constexpr size_t dim = 1;
      40             : 
      41           0 :   Affine(double A, double B, double a, double b);
      42             : 
      43           0 :   Affine() = default;
      44           0 :   ~Affine() = default;
      45           0 :   Affine(const Affine&) = default;
      46           0 :   Affine(Affine&&) = default;  // NOLINT
      47           0 :   Affine& operator=(const Affine&) = default;
      48           0 :   Affine& operator=(Affine&&) = default;
      49             : 
      50             :   template <typename T>
      51           0 :   std::array<tt::remove_cvref_wrap_t<T>, 1> operator()(
      52             :       const std::array<T, 1>& source_coords) const;
      53             : 
      54             :   /// The inverse function is only callable with doubles because the inverse
      55             :   /// might fail if called for a point out of range, and it is unclear
      56             :   /// what should happen if the inverse were to succeed for some points in a
      57             :   /// DataVector but fail for other points.
      58           1 :   std::optional<std::array<double, 1>> inverse(
      59             :       const std::array<double, 1>& target_coords) const;
      60             : 
      61             :   template <typename T>
      62           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 1, Frame::NoFrame> jacobian(
      63             :       const std::array<T, 1>& source_coords) const;
      64             : 
      65             :   template <typename T>
      66           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 1, Frame::NoFrame> inv_jacobian(
      67             :       const std::array<T, 1>& source_coords) const;
      68             : 
      69             :   // NOLINTNEXTLINE(google-runtime-references)
      70           0 :   void pup(PUP::er& p);
      71             : 
      72           0 :   bool is_identity() const { return is_identity_; }
      73             : 
      74             :  private:
      75           0 :   friend bool operator==(const Affine& lhs, const Affine& rhs);
      76             : 
      77           0 :   double A_{-1.0};
      78           0 :   double B_{1.0};
      79           0 :   double a_{-1.0};
      80           0 :   double b_{1.0};
      81           0 :   double length_of_domain_{2.0};  // B-A
      82           0 :   double length_of_range_{2.0};   // b-a
      83           0 :   double jacobian_{length_of_range_ / length_of_domain_};
      84           0 :   double inverse_jacobian_{length_of_domain_ / length_of_range_};
      85           0 :   bool is_identity_{false};
      86             : };
      87             : 
      88           0 : inline bool operator!=(const CoordinateMaps::Affine& lhs,
      89             :                        const CoordinateMaps::Affine& rhs) {
      90             :   return not(lhs == rhs);
      91             : }
      92             : 
      93             : }  // namespace CoordinateMaps
      94             : }  // namespace domain