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