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