Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines the class FlattOffsetWedge.
6 :
7 : #pragma once
8 :
9 : #include <array>
10 : #include <cstddef>
11 : #include <limits>
12 : #include <optional>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
16 :
17 : /// \cond
18 : namespace PUP {
19 : class er;
20 : } // namespace PUP
21 : /// \endcond
22 :
23 : namespace domain::CoordinateMaps {
24 :
25 : /*!
26 : * \ingroup CoordinateMapsGroup
27 : *
28 : * \brief Map from a cube to a volume that connects three planes and
29 : * a portion of one spherical surface.
30 : *
31 : * \image html FlatOffsetWedge.svg "Two slices through a FlatOffsetWedge."
32 : *
33 : * \details A cube is mapped to the volume shown in the figure.
34 : *
35 : * A $y=0$ slice through the volume is shown in the left
36 : * panel of the figure, and a $x=$const slice through the volume is shown
37 : * in the right panel of the figure.
38 : *
39 : * The upper-$z$ face of the volume is a portion of a spherical surface
40 : * of radius $R$ centered about point $C$ in the figure, which is along
41 : * the $x$-axis a distance $D$ from the origin (point $O$ in the figure).
42 : * The lower-$z$ face of the volume is a two-dimensional rectangle of
43 : * constant $z$, located a distance $L$ above the origin.
44 : * This rectangle extends from
45 : * $0 \leq x \leq D$ in the $x$ direction, and from
46 : * $-L \leq y \leq+L$ in the $y$ direction.
47 : *
48 : * The lower-$x$ face of the volume is a portion of the plane at constant
49 : * $x=0$, and the upper-$x$ face of the volume is a portion of the plane
50 : * at constant $x=D$. See the $x$ extents of the left panel of the
51 : * figure.
52 : *
53 : * Every constant-$x$ cross section of the volume looks like the right panel
54 : * of the figure: it is a two-dimensional wedge with 45 degree
55 : * opening angle and a base of
56 : * width $2L$. The outer radius of this wedge varies with $x$ (this radius
57 : * is $R$ at $x=D$ and $\sqrt{R^2-D^2}$ at $x=0$).
58 : *
59 : * ### Restrictions
60 : *
61 : * We require $D^2 + L^2 < R^2$ or else the lower-$x$ face of the mapped
62 : * volume lies outside of the sphere of radius $R$. This implies $L<R$
63 : * and $D<R$.
64 : * However, the condition that the upper-$z$ and lower-$z$ surface don't touch
65 : * each other at the $\xi=-1$, $\eta=\pm 1$ corners is more stringent:
66 : * $D^2 + 2L^2 < R^2$. This implies that $L<R/\sqrt{2}$.
67 : *
68 : * ## Equations for the map
69 : * Given our cube coordinates $\xi,\eta,\zeta$, each taking on values from
70 : * $-1$ to $+1$, we can derive the formulas for the map.
71 : *
72 : * Define
73 : * \begin{equation}
74 : * q \equiv \frac{D}{2R}.
75 : * \end{equation}
76 : *
77 : * Notice that $q<1/2$, because of the restriction $D<R$.
78 : *
79 : * First, note that $x$ is a function of $\xi$ only, for all points in
80 : * the cube:
81 : * \begin{align}
82 : * x(\xi,\eta,\zeta) &= q R (\xi+1).
83 : * \end{align}
84 : *
85 : * ### Surface map for bottom and top surfaces
86 : * The $y$ and $z$ coordinates for the mapped $\zeta=-1$ surface are
87 : * \begin{align}
88 : * \begin{bmatrix}
89 : * y(\xi,\eta,-1)\\
90 : * z(\xi,\eta,-1)
91 : * \end{bmatrix} = L
92 : * \begin{bmatrix}
93 : * \eta\\
94 : * 1
95 : * \end{bmatrix}.
96 : * \end{align}
97 : *
98 : * The $y$ and $z$ coordinates of the top surface,
99 : * $\zeta=+1$, of the volume are given by
100 : * \begin{align}
101 : * \begin{bmatrix}
102 : * y(\xi,\eta,+1)\\
103 : * z(\xi,\eta,+1)
104 : * \end{bmatrix} =
105 : * \frac{R}{\sqrt{1+\eta^2}}
106 : * \sqrt{1-q^2(\xi-1)^2}
107 : * \begin{bmatrix}
108 : * \eta \\
109 : * 1
110 : * \end{bmatrix}.
111 : * \end{align}
112 : *
113 : * ### Full volume map
114 : * Adding the $\zeta$ dependence by linear interpolation gives us
115 : * the full volume map for the $y$ and $z$ coordinates:
116 : *
117 : * \begin{align}
118 : * \begin{bmatrix}
119 : * y(\xi,\eta,\zeta)\\
120 : * z(\xi,\eta,\zeta)
121 : * \end{bmatrix} =
122 : * \left[L \frac{1-\zeta}{2}
123 : * + \frac{1+\zeta}{2}
124 : * \frac{R}{\sqrt{1+\eta^2}}
125 : * \sqrt{1-q^2(\xi-1)^2}\right]
126 : * \begin{bmatrix}
127 : * \eta \\
128 : * 1
129 : * \end{bmatrix}.
130 : * \end{align}
131 : *
132 : * ## Inverse
133 : *
134 : * The map can be inverted analytically:
135 : *
136 : * \begin{align}
137 : * \xi &= \frac{x}{qR}-1\\
138 : * \eta &= \frac{y}{z}\\
139 : * \zeta &= \frac{2z - L - P}{P-L},
140 : * \label{eq:zetaFromxyz}
141 : * \end{align}
142 : * where
143 : * \begin{align}
144 : * P &\equiv R\frac{\sqrt{1-(x-D)^2/R^2}}{\sqrt{1+y^2/z^2}}\\
145 : * &= R\frac{\sqrt{1-q^2(\xi-1)^2}}{\sqrt{1+\eta^2}}.
146 : * \label{eq:Pdefinition}
147 : * \end{align}
148 : *
149 : * It is easy to determine whether a point $(x,y,z)$ lies within the volume:
150 : * First compute $\xi$ and $\eta$ and check that they are both in $[-1,1]$.
151 : * If so, then the arguments of the square roots in
152 : * Eq. ($\ref{eq:Pdefinition}$) are guaranteed positive, and the denominator
153 : * of Eq. ($\ref{eq:zetaFromxyz}$) is guaranteed positive
154 : * by our condition $R^2>2L^2+D^2$,
155 : * so it is straightforward to
156 : * compute $\zeta$ and then check if it is in $[-1,1]$.
157 : *
158 : * ## Jacobian
159 : *
160 : * Straightforward differentiation gives
161 : *
162 : * \begin{align}
163 : * \frac{\partial x}{\partial \xi} &= qR\\
164 : * \frac{\partial x}{\partial \eta} &= 0\\
165 : * \frac{\partial x}{\partial \zeta} &= 0,
166 : * \end{align}
167 : *
168 : * \begin{align}
169 : * \partial_\zeta \begin{bmatrix}
170 : * y\\
171 : * z
172 : * \end{bmatrix} &= \frac{1}{2}
173 : * \left[\frac{R}{\sqrt{1+\eta^2}}\sqrt{1-q^2(\xi-1)^2}
174 : * -L\right]
175 : * \begin{bmatrix}
176 : * \eta \\
177 : * 1
178 : * \end{bmatrix}, \\
179 : * \partial_\xi \begin{bmatrix}
180 : * y\\
181 : * z
182 : * \end{bmatrix} &=
183 : * \left[\frac{Rq^2 (1+\zeta)(1-\xi)}{2\sqrt{1+\eta^2}}
184 : * \left(1-q^2(\xi-1)^2\right)^{-1/2}\right]
185 : * \begin{bmatrix}
186 : * \eta \\
187 : * 1
188 : * \end{bmatrix}, \\
189 : * \partial_\eta \begin{bmatrix}
190 : * y\\
191 : * z
192 : * \end{bmatrix} &=
193 : * -\left[\frac{R\eta(1+\zeta)}{2(1+\eta^2)^{3/2}}
194 : * \sqrt{1-q^2(\xi-1)^2}\right]
195 : * \begin{bmatrix}
196 : * \eta \\
197 : * 1
198 : * \end{bmatrix}+
199 : * \begin{bmatrix}
200 : * z \\
201 : * 0
202 : * \end{bmatrix}.
203 : * \end{align}
204 : *
205 : * ## Inverse Jacobian
206 : *
207 : * Let
208 : * \begin{align}
209 : * \Xi &\equiv P\frac{1+\zeta}{P-L}.
210 : * \end{align}
211 : * Then
212 : * \begin{align}
213 : * \frac{\partial \zeta}{\partial x} &=
214 : * \Xi q\frac{\xi-1}{R\sqrt{1-q^2(\xi-1)^2}},\\
215 : * \frac{\partial \zeta}{\partial y} &=
216 : * \Xi \frac{\eta}{z\sqrt{1+\eta^2}},\\
217 : * \frac{\partial \zeta}{\partial z} &=
218 : * -\Xi \eta\frac{\eta}{z\sqrt{1+\eta^2}}
219 : * +\frac{2}{P-L},\\
220 : * \frac{\partial \xi}{\partial x} &= \frac{1}{qR}\\
221 : * \frac{\partial \xi}{\partial y} &= 0\\
222 : * \frac{\partial \xi}{\partial z} &= 0\\
223 : * \frac{\partial \eta}{\partial x} &= 0\\
224 : * \frac{\partial \eta}{\partial y} &= \frac{1}{z}\\
225 : * \frac{\partial \eta}{\partial z} &= -\frac{\eta}{z}.
226 : * \end{align}
227 : *
228 : */
229 1 : class FlatOffsetWedge {
230 : public:
231 0 : static constexpr size_t dim = 3;
232 : /*!
233 : * \brief Constructs a FlatOffsetWedge.
234 : *
235 : * \param lower_face_y_half_width The half-width $L$ of the lower face in
236 : * the $y$ direction.
237 : * \param lower_face_x_width The width $D$ of the lower face in
238 : * the $x$ direction.
239 : * \param outer_radius The outer radius $R$.
240 : */
241 1 : FlatOffsetWedge(double lower_face_y_half_width, double lower_face_x_width,
242 : double outer_radius);
243 :
244 0 : FlatOffsetWedge() = default;
245 0 : ~FlatOffsetWedge() = default;
246 0 : FlatOffsetWedge(FlatOffsetWedge&&) = default;
247 0 : FlatOffsetWedge(const FlatOffsetWedge&) = default;
248 0 : FlatOffsetWedge& operator=(const FlatOffsetWedge&) = default;
249 0 : FlatOffsetWedge& operator=(FlatOffsetWedge&&) = default;
250 :
251 : template <typename T>
252 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
253 : const std::array<T, 3>& source_coords) const;
254 :
255 : /// The inverse function is only callable with doubles because the inverse
256 : /// might fail if called for a point out of range, and it is unclear
257 : /// what should happen if the inverse were to succeed for some points in a
258 : /// DataVector but fail for other points.
259 1 : std::optional<std::array<double, 3>> inverse(
260 : const std::array<double, 3>& target_coords) const;
261 :
262 : template <typename T>
263 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
264 : const std::array<T, 3>& source_coords) const;
265 :
266 : template <typename T>
267 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
268 : const std::array<T, 3>& source_coords) const;
269 :
270 : // clang-tidy: google runtime references
271 0 : void pup(PUP::er& p); // NOLINT
272 :
273 0 : static bool is_identity() { return false; }
274 :
275 : private:
276 0 : friend bool operator==(const FlatOffsetWedge& lhs,
277 : const FlatOffsetWedge& rhs);
278 0 : double lower_face_y_half_width_{std::numeric_limits<double>::signaling_NaN()};
279 0 : double lower_face_x_width_{std::numeric_limits<double>::signaling_NaN()};
280 0 : double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
281 : };
282 :
283 0 : bool operator!=(const FlatOffsetWedge& lhs, const FlatOffsetWedge& rhs);
284 : } // namespace domain::CoordinateMaps
|
