Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : // Defines the class BulgedCube.
5 :
6 : #pragma once
7 :
8 : #include <array>
9 : #include <cstddef>
10 : #include <limits>
11 : #include <optional>
12 :
13 : #include "DataStructures/Tensor/TypeAliases.hpp"
14 :
15 : /// \cond
16 : namespace PUP {
17 : class er;
18 : } // namespace PUP
19 : /// \endcond
20 :
21 : namespace domain {
22 : namespace CoordinateMaps {
23 :
24 : /*!
25 : * \ingroup CoordinateMapsGroup
26 : *
27 : * \brief Three dimensional map from the cube to a bulged cube.
28 : * The cube is shaped such that the surface is compatible
29 : * with the inner surface of Wedge<3>.
30 : * The shape of the object can be chosen to be cubical,
31 : * if the sphericity is set to 0, or to a sphere, if
32 : * the sphericity is set to 1. The sphericity can
33 : * be set to any number between 0 and 1 for a bulged cube.
34 : *
35 : * \details The volume map from the cube to a bulged cube is obtained by
36 : * interpolating between six surface maps, twelve bounding curves, and
37 : * eight corners. The surface map for the upper +z axis is obtained by
38 : * interpolating between a cubical surface and a spherical surface. The
39 : * two surfaces are chosen such that the latter circumscribes the former.
40 : *
41 : * We make a choice here as to whether we wish to use the logical coordinates
42 : * parameterizing these surface as they are, in which case we have the
43 : * equidistant choice of coordinates, or whether to apply a tangent map to them
44 : * which leads us to the equiangular choice of coordinates. In terms of the
45 : * logical coordinates, the equiangular coordinates are:
46 : *
47 : * \f[\textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4)\f]
48 : *
49 : * \f[\textrm{equiangular eta} : \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)\f]
50 : *
51 : * With derivatives:
52 : *
53 : * \f[\Xi'(\xi) = \frac{\pi}{4}(1+\Xi^2)\f]
54 : *
55 : * \f[\mathrm{H}'(\eta) = \frac{\pi}{4}(1+\mathrm{H}^2)\f]
56 : *
57 : * The equidistant coordinates are:
58 : *
59 : * \f[ \textrm{equidistant xi} : \Xi = \xi\f]
60 : *
61 : * \f[ \textrm{equidistant eta} : \mathrm{H} = \eta\f]
62 : *
63 : * with derivatives:
64 : *
65 : * <center>\f$\Xi'(\xi) = 1\f$, and \f$\mathrm{H}'(\eta) = 1\f$</center>
66 : *
67 : * We also define the variable \f$\rho\f$, given by:
68 : *
69 : * \f[\rho = \sqrt{1+\Xi^2+\mathrm{H}^2}\f]
70 : *
71 : * ### The Spherical Face Map
72 : * The surface map for the spherical face of radius \f$R\f$ lying in the
73 : * \f$+z\f$
74 : * direction in either choice of coordinates is then given by:
75 : *
76 : * \f[
77 : * \vec{\sigma}_{spherical}(\xi,\eta) =
78 : * \begin{bmatrix}
79 : * x(\xi,\eta)\\
80 : * y(\xi,\eta)\\
81 : * z(\xi,\eta)\\
82 : * \end{bmatrix} = \frac{R}{\rho}
83 : * \begin{bmatrix}
84 : * \Xi\\
85 : * \mathrm{H}\\
86 : * 1\\
87 : * \end{bmatrix}
88 : * \f]
89 : *
90 : * ### The Cubical Face Map
91 : * The surface map for the cubical face of side length \f$2L\f$ lying in the
92 : * \f$+z\f$ direction is given by:
93 : *
94 : * \f[
95 : * \vec{\sigma}_{cubical}(\xi,\eta) =
96 : * \begin{bmatrix}
97 : * x(\xi,\eta)\\
98 : * y(\xi,\eta)\\
99 : * L\\
100 : * \end{bmatrix} = L
101 : * \begin{bmatrix}
102 : * \Xi\\
103 : * \mathrm{H}\\
104 : * 1\\
105 : * \end{bmatrix}
106 : * \f]
107 : *
108 : * ### The Bulged Face Map
109 : * To construct the bulged map we interpolate between a cubical face map of
110 : * side length \f$2L\f$ and a spherical face map of radius \f$R\f$, with the
111 : * interpolation parameter being \f$s\f$, the `sphericity`.
112 : * The surface map for the bulged face lying in the \f$+z\f$ direction is then
113 : * given by:
114 : *
115 : * \f[
116 : * \vec{\sigma}_{+\zeta}(\xi,\eta) = \left\{(1-s)L + \frac{sR}{\rho}\right\}
117 : * \begin{bmatrix}
118 : * \Xi\\
119 : * \mathrm{H}\\
120 : * 1\\
121 : * \end{bmatrix}
122 : * \f]
123 : *
124 : * This equation defines the upper-z map \f$\vec{\sigma}_{+\zeta}\f$, and we
125 : * similarly define the other five surface maps \f$\vec{\sigma}_{+\eta}\f$,
126 : * \f$\vec{\sigma}_{+\xi}\f$, and so on by appropriate rotations.
127 : * We constrain L by demanding that the spherical face circumscribe the cube.
128 : * With this condition, we have \f$L = R/\sqrt3\f$.
129 : *
130 : * ### The General Formula for 3D Isoparametric Maps
131 : * The general formula is given by Eq. 1 in section 2.1 of Hesthaven's paper
132 : * "A Stable Penalty Method For The Compressible Navier-Stokes Equations III.
133 : * Multidimensional Domain Decomposition Schemes" available
134 : * <a href="
135 : * http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.699.1161&rep=rep1&type=pdf
136 : * "> here </a>.
137 : *
138 : * Hesthaven's formula is general in the degree of the shape functions used,
139 : * so for our purposes we take the special case where the shape functions are
140 : * linear in the interpolation variable, and define new variables accordingly.
141 : * However, our interpolation variables do not necessarily have to be the
142 : * logical coordinates themselves, though they often are. To make this
143 : * distinction clear, we will define the new interpolation variables
144 : * \f$\{\tilde{\xi},\tilde{\eta},\tilde{\zeta}\}\f$, which may either be the
145 : * logical coordinates themselves or a invertible transformation of them. For
146 : * the purposes of the bulged cube map, this transformation will be the same
147 : * transformation that takes the logical coordinates into the equiangular
148 : * coordinates. We will later see how this choice can lead to simplifications
149 : * in the final map.
150 : *
151 : * We define the following variables for
152 : * \f$\alpha, \beta, \gamma \in\{\tilde{\xi},\tilde{\eta},\tilde{\zeta}\}\f$:
153 : *
154 : * \f[
155 : * f^{\pm}_{\alpha} = \frac{1}{2}(1\pm\alpha)\\
156 : * f^{\pm\pm}_{\alpha \ \beta} = \frac{1}{4}(1\pm\alpha)(1\pm\beta)\\
157 : * f^{\pm\pm\pm}_{\alpha \ \beta \ \gamma} =
158 : * \frac{1}{8}(1\pm\alpha)(1\pm\beta)(1\pm\gamma)
159 : * \f]
160 : *
161 : * The formula involves six surfaces, which we will denote by
162 : * \f$\vec{\sigma}\f$, twelve curves, denoted by \f$\vec{\Gamma}\f$, and eight
163 : * vertices, denoted by \f$\vec{\pi}\f$, with the subscripts denoting which
164 : * face(s) these objects belong to. The full volume map is given by:
165 : *
166 : * \f{align*}
167 : * \vec{x}(\xi,\eta,\zeta) = &
168 : * f^{+}_{\tilde{\zeta}}\vec{\sigma}_{+\zeta}(\xi, \eta)+
169 : * f^{-}_{\tilde{\zeta}}\vec{\sigma}_{-\zeta}(\xi, \eta)\\
170 : * &+ f^{+}_{\tilde{\eta}}\vec{\sigma}_{+\eta}(\xi, \zeta)+
171 : * f^{-}_{\tilde{\eta}}\vec{\sigma}_{-\eta}(\xi, \zeta)+
172 : * f^{+}_{\tilde{\xi}}\vec{\sigma}_{+\xi}(\eta, \zeta)+
173 : * f^{-}_{\tilde{\xi}}\vec{\sigma}_{-\xi}(\eta, \zeta)\\
174 : * &- f^{++}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{+\xi +\eta}(\zeta)-
175 : * f^{-+}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{-\xi +\eta}(\zeta)-
176 : * f^{+-}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{+\xi -\eta}(\zeta)-
177 : * f^{--}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{-\xi -\eta}(\zeta)\\
178 : * &- f^{++}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{+\xi +\zeta}(\eta)-
179 : * f^{-+}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{-\xi +\zeta}(\eta)-
180 : * f^{+-}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{+\xi -\zeta}(\eta)-
181 : * f^{--}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{-\xi -\zeta}(\eta)\\
182 : * &- f^{++}_{\tilde{\eta} \ \tilde{\zeta}}\vec{\Gamma}_{+\eta +\zeta}(\xi)-
183 : * f^{-+}_{\tilde{\eta} \ \tilde{\zeta}}\vec{\Gamma}_{-\eta +\zeta}(\xi)-
184 : * f^{+-}_{\tilde{\eta} \ \tilde{\zeta}}\vec{\Gamma}_{+\eta -\zeta}(\xi)-
185 : * f^{--}_{\tilde{\eta} \tilde{\zeta}}\vec{\Gamma}_{-\eta -\zeta}(\xi)\\
186 : * &+ f^{+++}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{+\xi +\eta
187 : * +\zeta}+ f^{-++}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi
188 : * +\eta +\zeta}+ f^{+-+}_{\tilde{\xi} \ \tilde{\eta} \
189 : * \tilde{\zeta}}\vec{\pi}_{+\xi -\eta +\zeta}+
190 : * f^{--+}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi -\eta
191 : * +\zeta}\\
192 : * &+ f^{++-}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{+\xi +\eta
193 : * -\zeta}+ f^{-+-}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi
194 : * +\eta -\zeta}+ f^{+--}_{\tilde{\xi} \ \tilde{\eta} \
195 : * \tilde{\zeta}}\vec{\pi}_{+\xi -\eta -\zeta}+ f^{---}_{\tilde{\xi} \
196 : * \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi -\eta -\zeta} \f}
197 : *
198 : *
199 : * ### The Special Case for Octahedral Symmetry
200 : * The general formula is for the case in which there are six independently
201 : * specified bounding surfaces. In our case, the surfaces are obtained by
202 : * rotations and reflections of the upper-\f$\zeta\f$ face.
203 : *
204 : * We define the matrices corresponding to these transformations to be:
205 : *
206 : * \f[
207 : * S_{xy} =
208 : * \begin{bmatrix}
209 : * 0 & 1 & 0\\
210 : * 1 & 0 & 0\\
211 : * 0 & 0 & 1\\
212 : * \end{bmatrix},\
213 : *
214 : * S_{xz} =
215 : * \begin{bmatrix}
216 : * 0 & 0 & 1\\
217 : * 0 & 1 & 0\\
218 : * 1 & 0 & 0\\
219 : * \end{bmatrix},\
220 : *
221 : * S_{yz} =
222 : * \begin{bmatrix}
223 : * 1 & 0 & 0\\
224 : * 0 & 0 & 1\\
225 : * 0 & 1 & 0\\
226 : * \end{bmatrix}\f]
227 : *
228 : * \f[C_{zxy} =
229 : * \begin{bmatrix}
230 : * 0 & 0 & 1\\
231 : * 1 & 0 & 0\\
232 : * 0 & 1 & 0\\
233 : * \end{bmatrix},\
234 : *
235 : * C_{yzx} =
236 : * \begin{bmatrix}
237 : * 0 & 1 & 0\\
238 : * 0 & 0 & 1\\
239 : * 1 & 0 & 0\\
240 : * \end{bmatrix}\f]
241 : *
242 : * \f[N_{x} =
243 : * \begin{bmatrix}
244 : * -1 & 0 & 0\\
245 : * 0 & 1 & 0\\
246 : * 0 & 0 & 1\\
247 : * \end{bmatrix},\
248 : *
249 : * N_{y} =
250 : * \begin{bmatrix}
251 : * 1 & 0 & 0\\
252 : * 0 & -1 & 0\\
253 : * 0 & 0 & 1\\
254 : * \end{bmatrix},\
255 : *
256 : * N_{z} =
257 : * \begin{bmatrix}
258 : * 1 & 0 & 0\\
259 : * 0 & 1 & 0\\
260 : * 0 & 0 & -1\\
261 : * \end{bmatrix}
262 : * \f]
263 : *
264 : * The surface maps can now all be written in terms of
265 : * \f$\vec{\sigma}_{+\zeta}\f$ and these matrices:
266 : * <center>
267 : * \f$\vec{\sigma}_{-\zeta}(\xi, \eta) = N_z\vec{\sigma}_{+\zeta}(\xi, \eta)\\
268 : * \vec{\sigma}_{+\eta}(\xi, \zeta) = S_{yz}\vec{\sigma}_{+\zeta}(\xi, \zeta)\\
269 : * \vec{\sigma}_{-\eta}(\xi, \zeta) = N_yS_{yz}\vec{\sigma}_{+\zeta}(\xi,
270 : * \zeta)\\
271 : * \vec{\sigma}_{+\xi}(\eta, \zeta) = C_{zxy}\vec{\sigma}_{+\zeta}(\eta,
272 : * \zeta)\\
273 : * \vec{\sigma}_{-\xi}(\eta, \zeta) = N_xC_{zyx}\vec{\sigma}_{+\zeta}(\eta,
274 : * \zeta)\f$
275 : * </center>
276 : *
277 : * The four bounding curves \f$\vec{\Gamma}\f$ on the \f$+\zeta\f$ face are
278 : * given by:
279 : *
280 : * <center>
281 : * \f$\vec{\Gamma}_{+\xi,+\zeta}(\eta) = \vec{\sigma}_{+\zeta}(+1,\eta)\\
282 : * \vec{\Gamma}_{-\xi,+\zeta}(\eta) = \vec{\sigma}_{+\zeta}(-1,\eta)
283 : * = N_x\vec{\sigma}_{+\zeta}(+1, \eta)\\
284 : * \vec{\Gamma}_{+\eta,+\zeta}(\xi) = \vec{\sigma}_{+\zeta}(\xi,+1)
285 : * = S_{xy}\vec{\sigma}_{+\zeta}(+1, \xi)\\
286 : * \vec{\Gamma}_{-\eta,+\zeta}(\xi) = \vec{\sigma}_{+\zeta}(\xi,-1)
287 : * = N_yS_{xy}\vec{\sigma}_{+\zeta}(+1,\xi)\f$
288 : * </center>
289 : *
290 : * The bounding curves on the other surfaces can be obtained by transformations
291 : * on the \f$+\zeta\f$ face:
292 : *
293 : * <center>
294 : * \f$\vec{\Gamma}_{+\xi,-\zeta}(\eta) = N_z\vec{\sigma}_{+\zeta}(+1,\eta)\\
295 : * \vec{\Gamma}_{-\xi,-\zeta}(\eta) = N_z\vec{\sigma}_{+\zeta}(-1,\eta)
296 : * = N_zN_x\vec{\sigma}_{+\zeta}(+1,\eta)\\
297 : * \vec{\Gamma}_{+\eta,-\zeta}(\xi) = N_z\vec{\sigma}_{+\zeta}(\xi,+1)
298 : * = N_zS_{xy}\vec{\sigma}_{+\zeta}(+1, \xi)\\
299 : * \vec{\Gamma}_{-\eta,-\zeta}(\xi) = N_z\vec{\sigma}_{+\zeta}(\xi,-1)
300 : * = N_zN_yS_{xy}\vec{\sigma}_{+\zeta}(+1, \xi)\\
301 : * \vec{\Gamma}_{+\xi,+\eta}(\zeta) =
302 : * C_{zxy}\vec{\sigma}_{+\zeta}(+1,\zeta)\\
303 : * \vec{\Gamma}_{-\xi,+\eta}(\zeta) =
304 : * N_xC_{zxy}\vec{\sigma}_{+\zeta}(+1,\zeta)\\
305 : * \vec{\Gamma}_{+\xi,-\eta}(\zeta) = C_{zxy}\vec{\sigma}_{+\zeta}(-1,\zeta)
306 : * = C_{zxy}N_x\vec{\sigma}_{+\zeta}(+1,\zeta)\\
307 : * \vec{\Gamma}_{-\xi,-\eta}(\zeta) = N_xC_{zxy}\vec{\sigma}_{+\zeta}(-1,\zeta)
308 : * = N_xC_{zxy}N_x\vec{\sigma}_{+\zeta}(+1,\zeta)\f$
309 : * </center>
310 : *
311 : * Now we can write the volume map in terms of
312 : * \f$\vec{\sigma}_{+\zeta}\f$ only:
313 : * \f{align*}\vec{x}(\xi,\eta,\zeta) = &
314 : * (f^{+}_{\tilde{\zeta}} + f^{-}_{\tilde{\zeta}}N_z)
315 : * \vec{\sigma}_{+\zeta}(\xi, \eta)\\
316 : * &+ (f^{+}_{\tilde{\eta}} + f^{-}_{\tilde{\eta}}N_y)
317 : * S_{yz}\vec{\sigma}_{+\zeta}(\xi, \zeta)\\
318 : * &+ (f^{+}_{\tilde{\xi}} + f^{-}_{\tilde{\xi}}N_x)
319 : * C_{zxy}\vec{\sigma}_{+\zeta}(\eta, \zeta)\\
320 : * &- (f^{+}_{\tilde{\xi}}+f^{-}_{\tilde{\xi}}N_x)
321 : * (f^{+}_{\tilde{\eta}}+f^{-}_{\tilde{\eta}}N_y)
322 : * C_{zxy}\vec{\sigma}_{+\zeta}(+1, \zeta)\\
323 : * &- (f^{+}_{\tilde{\zeta}}+f^{-}_{\tilde{\zeta}}N_z)\left\{
324 : * (f^{+}_{\tilde{\xi}}+f^{-}_{\tilde{\xi}}N_x)\vec{\sigma}_{+\zeta}(+1, \eta)+
325 : * (f^{+}_{\tilde{\eta}}+f^{-}_{\tilde{\eta}}N_y)S_{xy}\vec{\sigma}_{+\zeta}(+1,
326 : * \xi)\right\}\\
327 : * &+ \frac{r}{\sqrt{3}}\vec{\tilde{\xi}}
328 : * \f}
329 : *
330 : * Note that we can now absorb all of the \f$f\f$s into the matrix prefactors
331 : * in the above equation and obtain a final set of matrices. We define the
332 : * following *blending matrices*:
333 : *
334 : * \f[
335 : * B_{\tilde{\xi}} =
336 : * \begin{bmatrix}
337 : * 0 & 0 & \tilde{\xi}\\
338 : * 1 & 0 & 0\\
339 : * 0 & 1 & 0\\
340 : * \end{bmatrix},\
341 : *
342 : * B_{\tilde{\eta}} =
343 : * \begin{bmatrix}
344 : * 1 & 0 & 0\\
345 : * 0 & 0 & \tilde{\eta}\\
346 : * 0 & 1 & 0\\
347 : * \end{bmatrix},\
348 : *
349 : * B_{\tilde{\zeta}} =
350 : * \begin{bmatrix}
351 : * 1 & 0 & 0\\
352 : * 0 & 1 & 0\\
353 : * 0 & 0 & \tilde{\zeta}\\
354 : * \end{bmatrix}\\
355 : *
356 : * B_{\tilde{\xi}\tilde{\eta}} =
357 : * \begin{bmatrix}
358 : * 0 & 0 & \tilde{\xi}\\
359 : * \tilde{\eta} & 0 & 0\\
360 : * 0 & 1 & 0\\
361 : * \end{bmatrix},\
362 : *
363 : * B_{\tilde{\xi}\tilde{\zeta}} =
364 : * \begin{bmatrix}
365 : * \tilde{\xi} & 0 & 0\\
366 : * 0 & 1 & 0\\
367 : * 0 & 0 & \tilde{\zeta}\\
368 : * \end{bmatrix},\
369 : *
370 : * B_{\tilde{\eta}\tilde{\zeta}} =
371 : * \begin{bmatrix}
372 : * 0 & 1 & 0\\
373 : * \tilde{\eta} & 0 & 0\\
374 : * 0 & 0 & \tilde{\zeta}\\
375 : * \end{bmatrix}\\
376 : *
377 : * B_{\tilde{\xi}\tilde{\eta}\tilde{\zeta}} =
378 : * \begin{bmatrix}
379 : * \tilde{\xi} & 0 & 0\\
380 : * 0 & \tilde{\eta} & 0\\
381 : * 0 & 0 & \tilde{\zeta}\\
382 : * \end{bmatrix}
383 : * \f]
384 : *
385 : * Now we can write the volume map in these terms:
386 : *
387 : * \f{align*}
388 : * \vec{x}(\xi,\eta,\zeta) = &
389 : * B_{\tilde{\zeta}}
390 : * \vec{\sigma}_{+\zeta}(\xi, \eta)\\& +
391 : * B_{\tilde{\eta}}
392 : * \vec{\sigma}_{+\zeta}(\xi, \zeta)+
393 : * B_{\tilde{\xi}}
394 : * \vec{\sigma}_{+\zeta}(\eta, \zeta)\\& -
395 : * B_{\tilde{\xi} \tilde{\eta}}
396 : * \vec{\sigma}_{+\zeta}(+1, \zeta)-
397 : * B_{\tilde{\xi} \tilde{\zeta}}
398 : * \vec{\sigma}_{+\zeta}(+1, \eta)+
399 : * B_{\tilde{\eta} \tilde{\zeta}}
400 : * \vec{\sigma}_{+\zeta}(+1, \xi)\\& +
401 : * B_{\tilde{\xi} \tilde{\eta} \tilde{\zeta}}
402 : * \vec{\sigma}_{+\zeta}(+1, +1)
403 : * \f}
404 : *
405 : * ### The Bulged Cube Map
406 : * We now use the result above to provide the mapping for the bulged cube.
407 : * First we will define the variables \f$\rho_A\f$ and \f$\rho_{AB}\f$, for
408 : * \f$A, B \in \{\Xi,\mathrm{H}, \mathrm{Z}\} \f$, where \f$\mathrm{Z}\f$
409 : * is \f$\tan(\zeta\pi/4)\f$ in the equiangular case and \f$\zeta\f$ in the
410 : * equidistant case:
411 : *
412 : * \f[
413 : * \rho_A = \sqrt{2 + A^2}\\
414 : * \rho_{AB} = \sqrt{1 + A^2 + B^2}
415 : * \f]
416 : * The final mapping is then:
417 : * \f[
418 : * \vec{x}(\xi,\eta,\zeta) = \frac{(1-s)R}{\sqrt{3}}
419 : * \begin{bmatrix}
420 : * \Xi\\
421 : * \mathrm{H}\\
422 : * \mathrm{Z}\\
423 : * \end{bmatrix} +
424 : * \frac{sR}{\sqrt{3}}
425 : * \begin{bmatrix}
426 : * \tilde{\xi}\\
427 : * \tilde{\eta}\\
428 : * \tilde{\zeta}\\
429 : * \end{bmatrix} + sR
430 : * \begin{bmatrix}
431 : * \tilde{\xi} & \Xi & \Xi\\
432 : * \mathrm{H} & \tilde{\eta} &\mathrm{H}\\
433 : * \mathrm{Z} & \mathrm{Z} & \tilde{\zeta}\\
434 : * \end{bmatrix}
435 : * \begin{bmatrix}
436 : * 1/\rho_{\mathrm{H}\mathrm{Z}}\\
437 : * 1/\rho_{\Xi\mathrm{Z}}\\
438 : * 1/\rho_{\Xi\mathrm{H}}\\
439 : * \end{bmatrix} - sR
440 : * \begin{bmatrix}
441 : * \Xi & \tilde{\xi} & \tilde{\xi}\\
442 : * \tilde{\eta} & \mathrm{H} &\tilde{\eta}\\
443 : * \tilde{\zeta} & \tilde{\zeta} & \mathrm{Z}\\
444 : * \end{bmatrix}
445 : * \begin{bmatrix}
446 : * 1/\rho_{\Xi}\\
447 : * 1/\rho_{\mathrm{H}}\\
448 : * 1/\rho_{\mathrm{Z}}\\
449 : * \end{bmatrix}
450 : * \f]
451 : *
452 : * Recall that the lower case Greek letters with tildes are the variables
453 : * used for the linear interpolation between the six bounding surfaces, and
454 : * that the upper case Greek letters are the coordinates along these surfaces -
455 : * both of which can be specified to be either
456 : * equidistant or equiangular. In the case where the
457 : * interpolation variable is chosen to match that of the
458 : * coordinates along the surface, we have \f$\tilde{\xi} = \Xi\f$, etc. In this
459 : * case, the formula reduces further. The reduced formula below is the one used
460 : * for this CoordinateMap. It is given by:
461 : *
462 : * \f[
463 : * \vec{x}(\xi,\eta,\zeta) =
464 : * \left\{
465 : * \frac{R}{\sqrt{3}}
466 : * + sR
467 : * \left(
468 : * 1/\rho_{\mathrm{H}\mathrm{Z}}+
469 : * 1/\rho_{\Xi\mathrm{Z}}+
470 : * 1/\rho_{\Xi\mathrm{H}}-
471 : * 1/\rho_{\Xi}-
472 : * 1/\rho_{\mathrm{H}}-
473 : * 1/\rho_{\mathrm{Z}}
474 : * \right)
475 : * \right\}
476 : * \begin{bmatrix}
477 : * \Xi\\
478 : * \mathrm{H}\\
479 : * \mathrm{Z}\\
480 : * \end{bmatrix}
481 : * \f]
482 : *
483 : * The inverse mapping is analytic in the angular directions. A root find
484 : * must be performed for the inverse mapping in the radial direction. This
485 : * one-dimensional formula is obtained by taking the magnitude of both sides
486 : * of the mapping, and changing variables from \f$\xi, \eta, \zeta\f$ to
487 : * \f$x, y, z\f$ and introducing \f$\rho^2 := \sqrt{\xi^2+\eta^2+\zeta^2}\f$.
488 : */
489 1 : class BulgedCube {
490 : public:
491 0 : static constexpr size_t dim = 3;
492 0 : BulgedCube(double radius, double sphericity, bool use_equiangular_map);
493 0 : BulgedCube() = default;
494 0 : ~BulgedCube() = default;
495 0 : BulgedCube(BulgedCube&&) = default;
496 0 : BulgedCube(const BulgedCube&) = default;
497 0 : BulgedCube& operator=(const BulgedCube&) = default;
498 0 : BulgedCube& operator=(BulgedCube&&) = default;
499 :
500 : template <typename T>
501 0 : std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
502 :
503 : /// The inverse function is only callable with doubles because the inverse
504 : /// might fail if called for a point out of range, and it is unclear
505 : /// what should happen if the inverse were to succeed for some points in a
506 : /// DataVector but fail for other points.
507 1 : std::optional<std::array<double, 3>> inverse(
508 : const std::array<double, 3>& target_coords) const;
509 :
510 : template <typename T>
511 0 : tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
512 : const std::array<T, 3>& source_coords) const;
513 :
514 : template <typename T>
515 0 : tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
516 : const std::array<T, 3>& source_coords) const;
517 :
518 : // NOLINTNEXTLINE(google-runtime-references)
519 0 : void pup(PUP::er& p);
520 :
521 0 : bool is_identity() const { return is_identity_; }
522 :
523 0 : static constexpr bool supports_hessian{true};
524 :
525 : private:
526 : /// Derivative of the physical coordinates with respect to the first of the
527 : /// three passed source coordinates. The Jacobian is assembled from calls
528 : /// with permuted arguments, exploiting the symmetry of the map.
529 : template <typename T>
530 1 : std::array<T, 3> xi_derivative(const T& xi, const T& eta,
531 : const T& zeta) const;
532 0 : friend bool operator==(const BulgedCube& lhs, const BulgedCube& rhs);
533 :
534 0 : double radius_{std::numeric_limits<double>::signaling_NaN()};
535 0 : double sphericity_{std::numeric_limits<double>::signaling_NaN()};
536 0 : bool use_equiangular_map_ = false;
537 0 : bool is_identity_ = false;
538 : };
539 :
540 0 : bool operator!=(const BulgedCube& lhs, const BulgedCube& rhs);
541 : } // namespace CoordinateMaps
542 : } // namespace domain
|
