Equiangular.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines the class Equiangular.
6 
7 #pragma once
8 
9 #include <array>
10 #include <boost/optional.hpp>
11 #include <cmath>
12 #include <cstddef>
13 
15 #include "Utilities/TypeTraits.hpp"
16 
17 namespace PUP {
18 class er;
19 } // namespace PUP
20 
21 namespace CoordinateMaps {
22 
23 /*!
24  * \ingroup CoordinateMapsGroup
25  * \brief Non-linear map from \f$\xi \in [A, B]\rightarrow x \in [a, b]\f$.
26  *
27  * The formula for the mapping is:
28  * \f{align}
29  * x &= \frac{a}{2} \left(1-\mathrm{tan}\left(
30  * \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right) +
31  * \frac{b}{2} \left(1+\mathrm{tan}\left(
32  * \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right)\\
33  * \xi &= \frac{A}{2} \left(1-\frac{4}{\pi}\mathrm{arctan}\left(
34  * \frac{2x-a-b}{b-a}\right)\right)+
35  * \frac{B}{2} \left(1+\frac{4}{\pi}\mathrm{arctan}\left(
36  * \frac{2x-a-b}{b-a}\right)\right)
37  * \f}
38  *
39  * \note The intermediate step in which a tangent map is applied can be more
40  * clearly understood if we define the coordinates:
41  * \f{align}
42  * \xi_{logical} &:= \frac{2\xi-B-A}{B-A} \in [-1, 1]\\
43  * \Xi &:= \mathrm{tan}\left(\frac{\pi\xi_{logical}}{4}\right) \in [-1, 1]
44  * \f}
45  *
46  * This map is intended to be used with `Wedge2D` and `Wedge3D` when equiangular
47  * coordinates are chosen for those maps. For more information on this choice
48  * of coordinates, see the documentation for Wedge3D.
49  */
50 class Equiangular {
51  public:
52  static constexpr size_t dim = 1;
53 
54  Equiangular(double A, double B, double a, double b) noexcept;
55 
56  Equiangular() = default;
57  ~Equiangular() = default;
58  Equiangular(const Equiangular&) = default;
59  Equiangular(Equiangular&&) noexcept = default;
60  Equiangular& operator=(const Equiangular&) = default;
61  Equiangular& operator=(Equiangular&&) = default;
62 
63  template <typename T>
65  const std::array<T, 1>& source_coords) const noexcept;
66 
67  boost::optional<std::array<double, 1>> inverse(
68  const std::array<double, 1>& target_coords) const noexcept;
69 
70  template <typename T>
71  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 1, Frame::NoFrame> jacobian(
72  const std::array<T, 1>& source_coords) const noexcept;
73 
74  template <typename T>
75  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 1, Frame::NoFrame> inv_jacobian(
76  const std::array<T, 1>& source_coords) const noexcept;
77 
78  // clang-tidy: google-runtime-references
79  void pup(PUP::er& p) noexcept; // NOLINT
80 
81  bool is_identity() const noexcept { return false; }
82 
83  private:
84  friend bool operator==(const Equiangular& lhs,
85  const Equiangular& rhs) noexcept;
86 
87  double A_{-1.0};
88  double B_{1.0};
89  double a_{-1.0};
90  double b_{1.0};
91  double length_of_domain_over_m_pi_4_{(B_ - A_) / M_PI_4}; // 4(B-A)/\pi
92  double length_of_range_{2.0}; // b-a
93  double m_pi_4_over_length_of_domain_{M_PI_4 / (B_ - A_)};
94  double one_over_length_of_range_{0.5};
95  // The jacobian for the affine map with the same parameters.
96  double linear_jacobian_times_m_pi_4_{length_of_range_ /
97  length_of_domain_over_m_pi_4_};
98  // The inverse jacobian for the affine map with the same parameters.
99  double linear_inverse_jacobian_over_m_pi_4_{length_of_domain_over_m_pi_4_ /
100  length_of_range_};
101 };
102 
103 inline bool operator!=(const CoordinateMaps::Equiangular& lhs,
104  const CoordinateMaps::Equiangular& rhs) noexcept {
105  return not(lhs == rhs);
106 }
107 
108 } // namespace CoordinateMaps
Definition: Strahlkorper.hpp:14
Contains all coordinate maps.
Definition: Affine.cpp:14
Defines a list of useful type aliases for tensors.
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
Non-linear map from .
Definition: Equiangular.hpp:50
Defines type traits, some of which are future STL type_traits header.