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