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