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             : 
      15             : /// \cond
      16             : namespace PUP {
      17             : class er;
      18             : }  // namespace PUP
      19             : /// \endcond
      20             : 
      21             : namespace domain {
      22           1 : namespace CoordinateMaps {
      23             : 
      24             : /*!
      25             :  * \ingroup CoordinateMapsGroup
      26             :  * \brief Affine map from \f$\xi \in [A, B]\rightarrow x \in [a, b]\f$.
      27             :  *
      28             :  * The formula for the mapping is...
      29             :  * \f[
      30             :  * x = \frac{b}{B-A} (\xi-A) +\frac{a}{B-A}(B-\xi)
      31             :  * \f]
      32             :  * \f[
      33             :  * \xi =\frac{B}{b-a} (x-a) +\frac{A}{b-a}(b-x)
      34             :  * \f]
      35             :  */
      36           1 : class Affine {
      37             :  public:
      38           0 :   static constexpr size_t dim = 1;
      39             : 
      40           0 :   Affine(double A, double B, double a, double b);
      41             : 
      42           0 :   Affine() = default;
      43           0 :   ~Affine() = default;
      44           0 :   Affine(const Affine&) = default;
      45           0 :   Affine(Affine&&) = default;  // NOLINT
      46           0 :   Affine& operator=(const Affine&) = default;
      47           0 :   Affine& operator=(Affine&&) = default;
      48             : 
      49             :   template <typename T>
      50           0 :   std::array<T, 1> operator()(const std::array<T, 1>& source_coords) const;
      51             : 
      52             :   /// The inverse function is only callable with doubles because the inverse
      53             :   /// might fail if called for a point out of range, and it is unclear
      54             :   /// what should happen if the inverse were to succeed for some points in a
      55             :   /// DataVector but fail for other points.
      56           1 :   std::optional<std::array<double, 1>> inverse(
      57             :       const std::array<double, 1>& target_coords) const;
      58             : 
      59             :   template <typename T>
      60           0 :   tnsr::Ij<T, 1, Frame::NoFrame> jacobian(
      61             :       const std::array<T, 1>& source_coords) const;
      62             : 
      63             :   template <typename T>
      64           0 :   tnsr::Ij<T, 1, Frame::NoFrame> inv_jacobian(
      65             :       const std::array<T, 1>& source_coords) const;
      66             : 
      67             :   // NOLINTNEXTLINE(google-runtime-references)
      68           0 :   void pup(PUP::er& p);
      69             : 
      70           0 :   bool is_identity() const { return is_identity_; }
      71             : 
      72           0 :   static constexpr bool supports_hessian{true};
      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