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