Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines the class FlatOffsetSphericalWedge.
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 four planes and
29 : * portions of two spherical surfaces.
30 : *
31 : * \image html FlatOffsetSphericalWedge.svg "Slices of FlatOffsetSphericalWedge"
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 lower-$z$ face of the volume is a portion of a spherical surface
40 : * of radius $R_1$ 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 upper-$z$ face of the volume is a portion of a spherical surface
43 : * of radius $R_2$ centered about the origin.
44 : *
45 : * The lower-$x$ face of the volume is a portion of the plane at constant
46 : * $x=0$, and the upper-$x$ face of the volume is a portion of the plane
47 : * at constant $x=D$. See the $x$ extents of the left panel of the
48 : * figure.
49 : *
50 : * Every constant-$x$ cross section of the volume looks like the right
51 : * panel of the figure: it is a two-dimensional wedge with 45 degree
52 : * opening angle. Both the inner and outer radii of this wedge vary
53 : * with $x$: The inner radius of the wedge is $R_1$ at $x=D$ and
54 : * $\sqrt{R_1^2-D^2}$ at $x=0$, and the outer radius of the wedge is
55 : * $R_2$ at $x=0$ and $\sqrt{R_2^2-D^2}$ at $x=D$.
56 : *
57 : * ### Relation to FlatOffsetWedge
58 : *
59 : * The lower-$z$ face of this FlatOffsetSphericalWedge is the same as
60 : * the upper-$z$ face of the similar map FlatOffsetWedge; Blocks using
61 : * these two maps are meant to abut at this surface. The parameter
62 : * $D$ means the same thing for the two maps, and the parameter $R$ in
63 : * FlatOffsetWedge is the same as what we call $R_1$ here.
64 : * For abutting blocks, the corresponding parameters of the two maps should
65 : * be equal, and both maps should have the same origin.
66 : *
67 : * The formulas below will be similar to those of the FlatOffsetWedge map,
68 : * but slightly more complicated.
69 : *
70 : * ### Restrictions
71 : *
72 : * We require $D<R_1$ or else the lower-$x$ face of the mapped
73 : * volume lies outside of the sphere of radius $R_1$.
74 : * We also require that $R_2^2 > R_1^2+D^2$ or else the outer sphere
75 : * and inner sphere will intersect at the upper-$x$ face and the map
76 : * will be singular.
77 : *
78 : * ## Equations for the map
79 : * Given our cube coordinates $\xi,\eta,\zeta$, each taking on values from
80 : * $-1$ to $+1$, we can derive the formulas for the map.
81 : *
82 : * Define
83 : * \begin{align}
84 : * q &\equiv \frac{D}{2R_1},\\
85 : * v &\equiv \frac{R_1}{R_2}.
86 : * \end{align}
87 : *
88 : * Notice that $q<1/2$, because of the restriction $D<R_1$. Also we
89 : * must have $v < \left(1+4q^2\right)^{-1/2}$ because of the
90 : * restriction $R_2^2 > R_1^2+D^2$.
91 : *
92 : * Note that $x$ is a function of $\xi$ only, for all points in the volume:
93 : * \begin{align}
94 : * x(\xi,\eta,\zeta) &= q R (\xi+1).
95 : * \end{align}
96 : *
97 : * ### Surface map for bottom and top surfaces
98 : * The $y$ and $z$ coordinates of the bottom 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_1}{\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 : * This is the same formula as for the upper face of the similar
113 : * FlatOffsetWedge map.
114 : *
115 : * The $y$ and $z$ coordinates for the mapped $\zeta=+1$ surface are
116 : * \begin{align}
117 : * \begin{bmatrix}
118 : * y(\xi,\eta,+1)\\
119 : * z(\xi,\eta,+1)
120 : * \end{bmatrix} =
121 : * \frac{R_2}{\sqrt{1+\eta^2}}
122 : * \sqrt{1-v^2q^2(\xi+1)^2}
123 : * \begin{bmatrix}
124 : * \eta \\
125 : * 1
126 : * \end{bmatrix}.
127 : * \end{align}
128 : * Because of the above restrictions on $v$ and $q$, and because $\xi$ and
129 : * $\eta$ range between $-1$ and $+1$, the arguments of the square
130 : * roots above are always positive.
131 : *
132 : * ### Full volume map
133 : * Adding the $\zeta$ dependence by linear interpolation gives us
134 : * the full volume map for the $y$ and $z$ coordinates:
135 : *
136 : * \begin{align}
137 : * \begin{bmatrix}
138 : * y(\xi,\eta,\zeta)\\
139 : * z(\xi,\eta,\zeta)
140 : * \end{bmatrix} =
141 : * \left[\frac{1+\zeta}{2}
142 : * \frac{R_2}{\sqrt{1+\eta^2}}
143 : * \sqrt{1-v^2q^2(\xi+1)^2}
144 : * + \frac{1-\zeta}{2}
145 : * \frac{R_1}{\sqrt{1+\eta^2}}
146 : * \sqrt{1-q^2(\xi-1)^2}\right]
147 : * \begin{bmatrix}
148 : * \eta \\
149 : * 1
150 : * \end{bmatrix}.
151 : * \end{align}
152 : *
153 : * ## Inverse
154 : *
155 : * The map can be inverted analytically:
156 : *
157 : * \begin{align}
158 : * \xi &= \frac{x}{qR}-1\\
159 : * \eta &= \frac{y}{z}\\
160 : * \zeta &= \frac{2z - W - P}{W-P},
161 : * \label{eq:zetaFromxyz}
162 : * \end{align}
163 : * where
164 : * \begin{align}
165 : * P &\equiv R_1\frac{\sqrt{1-(x-D)^2/R_1^2}}{\sqrt{1+y^2/z^2}}\\
166 : * &= R_1\frac{\sqrt{1-q^2(\xi-1)^2}}{\sqrt{1+\eta^2}},
167 : * \label{eq:Pdefinition} \\
168 : * W &\equiv R_2\frac{\sqrt{1-x^2/R_2^2}}{\sqrt{1+y^2/z^2}}\\
169 : * &= R_2\frac{\sqrt{1-v^2q^2(\xi+1)^2}}{\sqrt{1+\eta^2}}.
170 : * \label{eq:Wdefinition}
171 : * \end{align}
172 : *
173 : * It is easy to determine whether a point $(x,y,z)$ lies within the volume:
174 : * First compute $\xi$ and $\eta$ and check that they are both in $[-1,1]$.
175 : * If so, then the arguments of the square roots in
176 : * Eq. ($\ref{eq:Pdefinition}$) and Eq. ($\ref{eq:Wdefinition}$) are
177 : * guaranteed positive, and the denominator
178 : * of Eq. ($\ref{eq:zetaFromxyz}$) is guaranteed positive
179 : * by our condition $v < 1/\sqrt{1+4q^2}$.
180 : * Then it is straightforward to
181 : * compute $\zeta$ and then check if it is in $[-1,1]$.
182 : *
183 : * ## Jacobian
184 : *
185 : * Straightforward differentiation gives
186 : *
187 : * \begin{align}
188 : * \frac{\partial x}{\partial \xi} &= qR_1,\\
189 : * \frac{\partial x}{\partial \eta} &= 0,\\
190 : * \frac{\partial x}{\partial \zeta} &= 0,\\
191 : * \partial_\zeta \begin{bmatrix}
192 : * y\\
193 : * z
194 : * \end{bmatrix} &= \frac{W-P}{2}
195 : * \begin{bmatrix}
196 : * \eta \\
197 : * 1
198 : * \end{bmatrix}, \\
199 : * \partial_\xi \begin{bmatrix}
200 : * y\\
201 : * z
202 : * \end{bmatrix} &=
203 : * \frac{q^2}{2}\left[
204 : * P\frac{(1-\zeta)(1-\xi)}{1-q^2(1-\xi)^2}
205 : * -W\frac{v^2(1+\zeta)(1+\xi)}{1-v^2q^2(1+\xi)^2}
206 : * \right]
207 : * \begin{bmatrix}
208 : * \eta \\
209 : * 1
210 : * \end{bmatrix}, \\
211 : * \partial_\eta \begin{bmatrix}
212 : * y\\
213 : * z
214 : * \end{bmatrix} &=
215 : * -\frac{\eta}{2(1+\eta^2)}\left[
216 : * W(1+\zeta)
217 : * +P(1-\zeta)\right]
218 : * \begin{bmatrix}
219 : * \eta \\
220 : * 1
221 : * \end{bmatrix}+
222 : * \begin{bmatrix}
223 : * z \\
224 : * 0
225 : * \end{bmatrix}.
226 : * \end{align}
227 : *
228 : * Although it is straightforward to construct the inverse jacobian
229 : * analytically, for now we compute the inverse jacobian by numerically taking
230 : * the matrix inverse of the jacobian.
231 : *
232 : */
233 1 : class FlatOffsetSphericalWedge {
234 : public:
235 0 : static constexpr size_t dim = 3;
236 : /*!
237 : * \brief Constructs a FlatOffsetSphericalWedge.
238 : *
239 : * \param lower_face_x_width The width $D$ of the lower face in
240 : * the $x$ direction.
241 : * \param inner_radius The inner radius $R_1$.
242 : * \param outer_radius The outer radius $R_2$.
243 : */
244 1 : FlatOffsetSphericalWedge(double lower_face_x_width, double inner_radius,
245 : double outer_radius);
246 :
247 0 : FlatOffsetSphericalWedge() = default;
248 0 : ~FlatOffsetSphericalWedge() = default;
249 0 : FlatOffsetSphericalWedge(FlatOffsetSphericalWedge&&) = default;
250 0 : FlatOffsetSphericalWedge(const FlatOffsetSphericalWedge&) = default;
251 0 : FlatOffsetSphericalWedge& operator=(const FlatOffsetSphericalWedge&) =
252 : default;
253 0 : FlatOffsetSphericalWedge& operator=(FlatOffsetSphericalWedge&&) = default;
254 :
255 : template <typename T>
256 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
257 : const std::array<T, 3>& source_coords) const;
258 :
259 : /// The inverse function is only callable with doubles because the inverse
260 : /// might fail if called for a point out of range, and it is unclear
261 : /// what should happen if the inverse were to succeed for some points in a
262 : /// DataVector but fail for other points.
263 1 : std::optional<std::array<double, 3>> inverse(
264 : const std::array<double, 3>& target_coords) const;
265 :
266 : template <typename T>
267 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
268 : const std::array<T, 3>& source_coords) const;
269 :
270 : template <typename T>
271 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
272 : const std::array<T, 3>& source_coords) const;
273 :
274 : // clang-tidy: google runtime references
275 0 : void pup(PUP::er& p); // NOLINT
276 :
277 0 : static bool is_identity() { return false; }
278 :
279 : private:
280 0 : friend bool operator==(const FlatOffsetSphericalWedge& lhs,
281 : const FlatOffsetSphericalWedge& rhs);
282 0 : double lower_face_x_width_{std::numeric_limits<double>::signaling_NaN()};
283 0 : double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
284 0 : double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
285 : };
286 :
287 0 : bool operator!=(const FlatOffsetSphericalWedge& lhs,
288 : const FlatOffsetSphericalWedge& rhs);
289 : } // namespace domain::CoordinateMaps
|
