Wedge2D.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 Wedge2D.
6 
7 #pragma once
8 
9 #include <array>
10 #include <boost/optional.hpp>
11 #include <cstddef>
12 #include <limits>
13 
15 #include "Domain/OrientationMap.hpp"
16 #include "Utilities/TypeTraits.hpp"
17 
18 namespace PUP {
19 class er;
20 } // namespace PUP
21 
22 namespace domain {
23 namespace CoordinateMaps {
24 
25 /*!
26  * \ingroup CoordinateMapsGroup
27  *
28  * \brief Two dimensional map from the logical square to a wedge, which is used
29  * by domains which use disks or annuli.
30  *
31  * \details The coordinate map is constructed by
32  * linearly interpolating between a bulged arc which is circumscribed by a
33  * circular arc of radius `radius_outer` and a bulged arc which is
34  * circumscribed by a circular arc of radius `radius_inner`. These arcs can be
35  * made to be straight or circular, based on the value of `circularity_inner`
36  * or `circularity_outer`, which can take on values between 0 (for straight)
37  * and 1 (for circular). These arcs extend \f$\pi/2\f$ in angle, and can be
38  * oriented along the +/- x or y axis. The choice of using either equiangular
39  * or equidistant coordinates along the arcs is specifiable with
40  * `with_equiangular_map`. The default logical coordinate that points in the
41  * angular direction is \f$\eta\f$. We introduce the auxiliary variable
42  * \f$\mathrm{H}\f$ which is a function of \f$\eta\f$. If we are using
43  * equiangular coordinates, we have:
44  *
45  * \f[\mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)\f]
46  *
47  * With derivative:
48  *
49  * \f[\mathrm{H}'(\eta) = \frac{\pi}{4}(1+\mathrm{H}^2)\f]
50  *
51  * If we are using equidistant coordinates, we have:
52  *
53  * \f[\mathrm{H}(\eta) = \eta\f]
54  *
55  * with derivative:
56  *
57  * <center>\f$\mathrm{H}'(\eta) = 1\f$</center>
58  *
59  * We also define the variable \f$\rho\f$, given by:
60  *
61  * \f[\rho = \sqrt{1+\mathrm{H}^2}\f]
62  *
63  * In terms of the the circularity \f$c\f$ and the radius \f$R\f$,
64  * the mapping is:
65  *
66  * \f[\vec{x}(\xi,\eta) =
67  * \frac{1}{2}\left\{(1-\xi)\Big[(1-c_{inner})\frac{R_{inner}}{\sqrt 2}
68  * + c_{inner}\frac{R_{inner}}{\rho}\Big] +
69  * (1+\xi)\Big[(1-c_{outer})\frac{R_{outer}}{\sqrt 2} +c_{outer}
70  * \frac{R_{outer}}{\rho}\Big] \right\}\begin{bmatrix}
71  * 1\\
72  * \mathrm{H}\\
73  * \end{bmatrix}\f]
74  *
75  * We will define the variables \f$T(\xi)\f$ and \f$A(\xi)\f$, the trapezoid
76  * and annulus factors: \f[T(\xi) = T_0 + T_1\xi\f] \f[A(\xi) = A_0 +
77  * A_1\xi\f]
78  * Where \f{align*}T_0 &= \frac{1}{2} \big\{ (1-c_{outer})R_{outer} +
79  * (1-c_{inner})R_{inner}\big\}\\
80  * T_1 &= \partial_{\xi} T = \frac{1}{2} \big\{ (1-c_{outer})R_{outer} -
81  * (1-c_{inner})R_{inner}\big\}\\
82  * A_0 &= \frac{1}{2} \big\{ c_{outer}R_{outer} + c_{inner}R_{inner}\big\}\\
83  * A_1 &= \partial_{\xi} A = \frac{1}{2} \big\{ c_{outer}R_{outer} -
84  * c_{inner}R_{inner}\big\}\f}
85  *
86  * The map can then be rewritten as:
87  * \f[\vec{x}(\xi,\eta) = \left\{\frac{T(\xi)}{\sqrt 2} +
88  * \frac{A(\xi)}{\rho}\right\}\begin{bmatrix}
89  * 1\\
90  * \mathrm{H}\\
91  * \end{bmatrix}\f]
92  *
93  *
94  * The Jacobian is: \f[J =
95  * \begin{bmatrix}
96  * \frac{T_1}{\sqrt 2} + \frac{A_1}{\rho} &
97  * \mathrm{H}\mathrm{H}'\frac{A(\xi)}{\rho^3} \\
98  * \mathrm{H}\partial_{\xi}x &
99  * \mathrm{H}\partial_{\eta}x + \mathrm{H}'x\\
100  * \end{bmatrix}
101  * \f]
102  *
103  * The inverse Jacobian is: \f[J^{-1} =
104  * \frac{1}{x}\begin{bmatrix}
105  * \frac{1}{\partial_{\xi}x}\Big\{
106  * \frac{T(\xi)}{\sqrt 2}+\frac{A(\xi)}{\rho^3}
107  * \Big\} & \mathrm{H}\frac{1}{\partial_{\xi}x}\frac{A(\xi)}{\rho^3}\\
108  * -\mathrm{H}\mathrm{H}'^{-1} & \mathrm{H}'^{-1}\\
109  * \end{bmatrix}
110  * \f]
111  *
112  * For a more detailed discussion, see the
113  * documentation for Wedge3D, where the default logical coordinate that points
114  * in the radial direction is \f$\zeta\f$ (In Wedge2D the logical coordinate
115  * that points in the radial direction is \f$\xi\f$), and one may set either
116  * of the two other logical coordinates to zero to obtain an equivalent
117  * Wedge2D map.
118  */
119 
120 class Wedge2D {
121  public:
122  static constexpr size_t dim = 2;
123 
124  Wedge2D(double radius_inner, double radius_outer, double circularity_inner,
125  double circularity_outer, OrientationMap<2> orientation_of_wedge,
126  bool with_equiangular_map) noexcept;
127 
128  Wedge2D() = default;
129  ~Wedge2D() = default;
130  Wedge2D(Wedge2D&&) = default;
131  Wedge2D& operator=(Wedge2D&&) = default;
132  Wedge2D(const Wedge2D&) = default;
133  Wedge2D& operator=(const Wedge2D&) = default;
134 
135  template <typename T>
137  const std::array<T, 2>& source_coords) const noexcept;
138 
139  /// Returns invalid if \f$x<=0\f$ (for a \f$+x\f$-oriented `Wedge2D`).
140  boost::optional<std::array<double, 2>> inverse(
141  const std::array<double, 2>& target_coords) const noexcept;
142 
143  template <typename T>
144  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 2, Frame::NoFrame> jacobian(
145  const std::array<T, 2>& source_coords) const noexcept;
146 
147  template <typename T>
148  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 2, Frame::NoFrame> inv_jacobian(
149  const std::array<T, 2>& source_coords) const noexcept;
150 
151  // clang-tidy: google runtime references
152  void pup(PUP::er& p); // NOLINT
153 
154  bool is_identity() const noexcept { return false; }
155 
156  private:
157  friend bool operator==(const Wedge2D& lhs, const Wedge2D& rhs) noexcept;
158 
159  double radius_inner_{};
160  double radius_outer_{};
161  double circularity_inner_{};
162  double circularity_outer_{};
163  OrientationMap<2> orientation_of_wedge_{};
164  bool with_equiangular_map_ = false;
165  double scaled_trapezoid_zero_{std::numeric_limits<double>::signaling_NaN()};
166  double annulus_zero_{std::numeric_limits<double>::signaling_NaN()};
167  double scaled_trapezoid_rate_{std::numeric_limits<double>::signaling_NaN()};
168  double annulus_rate_{std::numeric_limits<double>::signaling_NaN()};
169 };
170 bool operator!=(const Wedge2D& lhs, const Wedge2D& rhs) noexcept;
171 } // namespace CoordinateMaps
172 } // namespace domain
Definition: Strahlkorper.hpp:14
Definition: BlockId.hpp:16
T signaling_NaN(T... args)
Two dimensional map from the logical square to a wedge, which is used by domains which use disks or a...
Definition: Wedge2D.hpp:120
Defines classes for Tensor.
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
Defines type traits, some of which are future STL type_traits header.