|
SpECTRE
v2024.04.12
|
Three dimensional map from the cube to a bulged cube. The cube is shaped such that the surface is compatible with the inner surface of Wedge<3>. The shape of the object can be chosen to be cubical, if the sphericity is set to 0, or to a sphere, if the sphericity is set to 1. The sphericity can be set to any number between 0 and 1 for a bulged cube. More...
#include <BulgedCube.hpp>
Public Member Functions | |
| BulgedCube (double radius, double sphericity, bool use_equiangular_map) | |
| BulgedCube (BulgedCube &&)=default | |
| BulgedCube (const BulgedCube &)=default | |
| BulgedCube & | operator= (const BulgedCube &)=default |
| BulgedCube & | operator= (BulgedCube &&)=default |
| template<typename T > | |
| std::array< tt::remove_cvref_wrap_t< T >, 3 > | operator() (const std::array< T, 3 > &source_coords) const |
| std::optional< std::array< double, 3 > > | inverse (const std::array< double, 3 > &target_coords) const |
| The inverse function is only callable with doubles because the inverse might fail if called for a point out of range, and it is unclear what should happen if the inverse were to succeed for some points in a DataVector but fail for other points. | |
| template<typename T > | |
| tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > | jacobian (const std::array< T, 3 > &source_coords) const |
| template<typename T > | |
| tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > | inv_jacobian (const std::array< T, 3 > &source_coords) const |
| void | pup (PUP::er &p) |
| bool | is_identity () const |
Static Public Attributes | |
| static constexpr size_t | dim = 3 |
Friends | |
| bool | operator== (const BulgedCube &lhs, const BulgedCube &rhs) |
Three dimensional map from the cube to a bulged cube. The cube is shaped such that the surface is compatible with the inner surface of Wedge<3>. The shape of the object can be chosen to be cubical, if the sphericity is set to 0, or to a sphere, if the sphericity is set to 1. The sphericity can be set to any number between 0 and 1 for a bulged cube.
The volume map from the cube to a bulged cube is obtained by interpolating between six surface maps, twelve bounding curves, and eight corners. The surface map for the upper +z axis is obtained by interpolating between a cubical surface and a spherical surface. The two surfaces are chosen such that the latter circumscribes the former.
We make a choice here as to whether we wish to use the logical coordinates parameterizing these surface as they are, in which case we have the equidistant choice of coordinates, or whether to apply a tangent map to them which leads us to the equiangular choice of coordinates. In terms of the logical coordinates, the equiangular coordinates are:
\[\textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4)\]
\[\textrm{equiangular eta} : \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)\]
With derivatives:
\[\Xi'(\xi) = \frac{\pi}{4}(1+\Xi^2)\]
\[\mathrm{H}'(\eta) = \frac{\pi}{4}(1+\mathrm{H}^2)\]
The equidistant coordinates are:
\[ \textrm{equidistant xi} : \Xi = \xi\]
\[ \textrm{equidistant eta} : \mathrm{H} = \eta\]
with derivatives:
We also define the variable \(\rho\), given by:
\[\rho = \sqrt{1+\Xi^2+\mathrm{H}^2}\]
The surface map for the spherical face of radius \(R\) lying in the \(+z\) direction in either choice of coordinates is then given by:
\[ \vec{\sigma}_{spherical}(\xi,\eta) = \begin{bmatrix} x(\xi,\eta)\\ y(\xi,\eta)\\ z(\xi,\eta)\\ \end{bmatrix} = \frac{R}{\rho} \begin{bmatrix} \Xi\\ \mathrm{H}\\ 1\\ \end{bmatrix} \]
The surface map for the cubical face of side length \(2L\) lying in the \(+z\) direction is given by:
\[ \vec{\sigma}_{cubical}(\xi,\eta) = \begin{bmatrix} x(\xi,\eta)\\ y(\xi,\eta)\\ L\\ \end{bmatrix} = L \begin{bmatrix} \Xi\\ \mathrm{H}\\ 1\\ \end{bmatrix} \]
To construct the bulged map we interpolate between a cubical face map of side length \(2L\) and a spherical face map of radius \(R\), with the interpolation parameter being \(s\), the sphericity. The surface map for the bulged face lying in the \(+z\) direction is then given by:
\[ \vec{\sigma}_{+\zeta}(\xi,\eta) = \left\{(1-s)L + \frac{sR}{\rho}\right\} \begin{bmatrix} \Xi\\ \mathrm{H}\\ 1\\ \end{bmatrix} \]
This equation defines the upper-z map \(\vec{\sigma}_{+\zeta}\), and we similarly define the other five surface maps \(\vec{\sigma}_{+\eta}\), \(\vec{\sigma}_{+\xi}\), and so on by appropriate rotations. We constrain L by demanding that the spherical face circumscribe the cube. With this condition, we have \(L = R/\sqrt3\).
The general formula is given by Eq. 1 in section 2.1 of Hesthaven's paper "A Stable Penalty Method For The Compressible Navier-Stokes Equations III. Multidimensional Domain Decomposition Schemes" available here .
Hesthaven's formula is general in the degree of the shape functions used, so for our purposes we take the special case where the shape functions are linear in the interpolation variable, and define new variables accordingly. However, our interpolation variables do not necessarily have to be the logical coordinates themselves, though they often are. To make this distinction clear, we will define the new interpolation variables \(\{\tilde{\xi},\tilde{\eta},\tilde{\zeta}\}\), which may either be the logical coordinates themselves or a invertible transformation of them. For the purposes of the bulged cube map, this transformation will be the same transformation that takes the logical coordinates into the equiangular coordinates. We will later see how this choice can lead to simplifications in the final map.
We define the following variables for \(\alpha, \beta, \gamma \in\{\tilde{\xi},\tilde{\eta},\tilde{\zeta}\}\):
\[ f^{\pm}_{\alpha} = \frac{1}{2}(1\pm\alpha)\\ f^{\pm\pm}_{\alpha \ \beta} = \frac{1}{4}(1\pm\alpha)(1\pm\beta)\\ f^{\pm\pm\pm}_{\alpha \ \beta \ \gamma} = \frac{1}{8}(1\pm\alpha)(1\pm\beta)(1\pm\gamma) \]
The formula involves six surfaces, which we will denote by \(\vec{\sigma}\), twelve curves, denoted by \(\vec{\Gamma}\), and eight vertices, denoted by \(\vec{\pi}\), with the subscripts denoting which face(s) these objects belong to. The full volume map is given by:
\begin{align*} \vec{x}(\xi,\eta,\zeta) = & f^{+}_{\tilde{\zeta}}\vec{\sigma}_{+\zeta}(\xi, \eta)+ f^{-}_{\tilde{\zeta}}\vec{\sigma}_{-\zeta}(\xi, \eta)\\ &+ f^{+}_{\tilde{\eta}}\vec{\sigma}_{+\eta}(\xi, \zeta)+ f^{-}_{\tilde{\eta}}\vec{\sigma}_{-\eta}(\xi, \zeta)+ f^{+}_{\tilde{\xi}}\vec{\sigma}_{+\xi}(\eta, \zeta)+ f^{-}_{\tilde{\xi}}\vec{\sigma}_{-\xi}(\eta, \zeta)\\ &- f^{++}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{+\xi +\eta}(\zeta)- f^{-+}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{-\xi +\eta}(\zeta)- f^{+-}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{+\xi -\eta}(\zeta)- f^{--}_{\tilde{\xi} \ \tilde{\eta}}\vec{\Gamma}_{-\xi -\eta}(\zeta)\\ &- f^{++}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{+\xi +\zeta}(\eta)- f^{-+}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{-\xi +\zeta}(\eta)- f^{+-}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{+\xi -\zeta}(\eta)- f^{--}_{\tilde{\xi} \ \tilde{\zeta}}\vec{\Gamma}_{-\xi -\zeta}(\eta)\\ &- f^{++}_{\tilde{\eta} \ \tilde{\zeta}}\vec{\Gamma}_{+\eta +\zeta}(\xi)- f^{-+}_{\tilde{\eta} \ \tilde{\zeta}}\vec{\Gamma}_{-\eta +\zeta}(\xi)- f^{+-}_{\tilde{\eta} \ \tilde{\zeta}}\vec{\Gamma}_{+\eta -\zeta}(\xi)- f^{--}_{\tilde{\eta} \tilde{\zeta}}\vec{\Gamma}_{-\eta -\zeta}(\xi)\\ &+ f^{+++}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{+\xi +\eta +\zeta}+ f^{-++}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi +\eta +\zeta}+ f^{+-+}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{+\xi -\eta +\zeta}+ f^{--+}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi -\eta +\zeta}\\ &+ f^{++-}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{+\xi +\eta -\zeta}+ f^{-+-}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi +\eta -\zeta}+ f^{+--}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{+\xi -\eta -\zeta}+ f^{---}_{\tilde{\xi} \ \tilde{\eta} \ \tilde{\zeta}}\vec{\pi}_{-\xi -\eta -\zeta} \end{align*}
The general formula is for the case in which there are six independently specified bounding surfaces. In our case, the surfaces are obtained by rotations and reflections of the upper- \(\zeta\) face.
We define the matrices corresponding to these transformations to be:
\[ S_{xy} = \begin{bmatrix} 0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1\\ \end{bmatrix},\ S_{xz} = \begin{bmatrix} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0\\ \end{bmatrix},\ S_{yz} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0\\ \end{bmatrix}\]
\[C_{zxy} = \begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0\\ \end{bmatrix},\ C_{yzx} = \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ 1 & 0 & 0\\ \end{bmatrix}\]
\[N_{x} = \begin{bmatrix} -1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{bmatrix},\ N_{y} = \begin{bmatrix} 1 & 0 & 0\\ 0 & -1 & 0\\ 0 & 0 & 1\\ \end{bmatrix},\ N_{z} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1\\ \end{bmatrix} \]
The surface maps can now all be written in terms of \(\vec{\sigma}_{+\zeta}\) and these matrices:
The four bounding curves \(\vec{\Gamma}\) on the \(+\zeta\) face are given by:
The bounding curves on the other surfaces can be obtained by transformations on the \(+\zeta\) face:
Now we can write the volume map in terms of \(\vec{\sigma}_{+\zeta}\) only:
\begin{align*}\vec{x}(\xi,\eta,\zeta) = & (f^{+}_{\tilde{\zeta}} + f^{-}_{\tilde{\zeta}}N_z) \vec{\sigma}_{+\zeta}(\xi, \eta)\\ &+ (f^{+}_{\tilde{\eta}} + f^{-}_{\tilde{\eta}}N_y) S_{yz}\vec{\sigma}_{+\zeta}(\xi, \zeta)\\ &+ (f^{+}_{\tilde{\xi}} + f^{-}_{\tilde{\xi}}N_x) C_{zxy}\vec{\sigma}_{+\zeta}(\eta, \zeta)\\ &- (f^{+}_{\tilde{\xi}}+f^{-}_{\tilde{\xi}}N_x) (f^{+}_{\tilde{\eta}}+f^{-}_{\tilde{\eta}}N_y) C_{zxy}\vec{\sigma}_{+\zeta}(+1, \zeta)\\ &- (f^{+}_{\tilde{\zeta}}+f^{-}_{\tilde{\zeta}}N_z)\left\{ (f^{+}_{\tilde{\xi}}+f^{-}_{\tilde{\xi}}N_x)\vec{\sigma}_{+\zeta}(+1, \eta)+ (f^{+}_{\tilde{\eta}}+f^{-}_{\tilde{\eta}}N_y)S_{xy}\vec{\sigma}_{+\zeta}(+1, \xi)\right\}\\ &+ \frac{r}{\sqrt{3}}\vec{\tilde{\xi}} \end{align*}
Note that we can now absorb all of the \(f\)s into the matrix prefactors in the above equation and obtain a final set of matrices. We define the following blending matrices:
\[ B_{\tilde{\xi}} = \begin{bmatrix} 0 & 0 & \tilde{\xi}\\ 1 & 0 & 0\\ 0 & 1 & 0\\ \end{bmatrix},\ B_{\tilde{\eta}} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 0 & \tilde{\eta}\\ 0 & 1 & 0\\ \end{bmatrix},\ B_{\tilde{\zeta}} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & \tilde{\zeta}\\ \end{bmatrix}\\ B_{\tilde{\xi}\tilde{\eta}} = \begin{bmatrix} 0 & 0 & \tilde{\xi}\\ \tilde{\eta} & 0 & 0\\ 0 & 1 & 0\\ \end{bmatrix},\ B_{\tilde{\xi}\tilde{\zeta}} = \begin{bmatrix} \tilde{\xi} & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & \tilde{\zeta}\\ \end{bmatrix},\ B_{\tilde{\eta}\tilde{\zeta}} = \begin{bmatrix} 0 & 1 & 0\\ \tilde{\eta} & 0 & 0\\ 0 & 0 & \tilde{\zeta}\\ \end{bmatrix}\\ B_{\tilde{\xi}\tilde{\eta}\tilde{\zeta}} = \begin{bmatrix} \tilde{\xi} & 0 & 0\\ 0 & \tilde{\eta} & 0\\ 0 & 0 & \tilde{\zeta}\\ \end{bmatrix} \]
Now we can write the volume map in these terms:
\begin{align*} \vec{x}(\xi,\eta,\zeta) = & B_{\tilde{\zeta}} \vec{\sigma}_{+\zeta}(\xi, \eta)\\& + B_{\tilde{\eta}} \vec{\sigma}_{+\zeta}(\xi, \zeta)+ B_{\tilde{\xi}} \vec{\sigma}_{+\zeta}(\eta, \zeta)\\& - B_{\tilde{\xi} \tilde{\eta}} \vec{\sigma}_{+\zeta}(+1, \zeta)- B_{\tilde{\xi} \tilde{\zeta}} \vec{\sigma}_{+\zeta}(+1, \eta)+ B_{\tilde{\eta} \tilde{\zeta}} \vec{\sigma}_{+\zeta}(+1, \xi)\\& + B_{\tilde{\xi} \tilde{\eta} \tilde{\zeta}} \vec{\sigma}_{+\zeta}(+1, +1) \end{align*}
We now use the result above to provide the mapping for the bulged cube. First we will define the variables \(\rho_A\) and \(\rho_{AB}\), for \(A, B \in \{\Xi,\mathrm{H}, \mathrm{Z}\} \), where \(\mathrm{Z}\) is \(\tan(\zeta\pi/4)\) in the equiangular case and \(\zeta\) in the equidistant case:
\[ \rho_A = \sqrt{2 + A^2}\\ \rho_{AB} = \sqrt{1 + A^2 + B^2} \]
The final mapping is then:
\[ \vec{x}(\xi,\eta,\zeta) = \frac{(1-s)R}{\sqrt{3}} \begin{bmatrix} \Xi\\ \mathrm{H}\\ \mathrm{Z}\\ \end{bmatrix} + \frac{sR}{\sqrt{3}} \begin{bmatrix} \tilde{\xi}\\ \tilde{\eta}\\ \tilde{\zeta}\\ \end{bmatrix} + sR \begin{bmatrix} \tilde{\xi} & \Xi & \Xi\\ \mathrm{H} & \tilde{\eta} &\mathrm{H}\\ \mathrm{Z} & \mathrm{Z} & \tilde{\zeta}\\ \end{bmatrix} \begin{bmatrix} 1/\rho_{\mathrm{H}\mathrm{Z}}\\ 1/\rho_{\Xi\mathrm{Z}}\\ 1/\rho_{\Xi\mathrm{H}}\\ \end{bmatrix} - sR \begin{bmatrix} \Xi & \tilde{\xi} & \tilde{\xi}\\ \tilde{\eta} & \mathrm{H} &\tilde{\eta}\\ \tilde{\zeta} & \tilde{\zeta} & \mathrm{Z}\\ \end{bmatrix} \begin{bmatrix} 1/\rho_{\Xi}\\ 1/\rho_{\mathrm{H}}\\ 1/\rho_{\mathrm{Z}}\\ \end{bmatrix} \]
Recall that the lower case Greek letters with tildes are the variables used for the linear interpolation between the six bounding surfaces, and that the upper case Greek letters are the coordinates along these surfaces - both of which can be specified to be either equidistant or equiangular. In the case where the interpolation variable is chosen to match that of the coordinates along the surface, we have \(\tilde{\xi} = \Xi\), etc. In this case, the formula reduces further. The reduced formula below is the one used for this CoordinateMap. It is given by:
\[ \vec{x}(\xi,\eta,\zeta) = \left\{ \frac{R}{\sqrt{3}} + sR \left( 1/\rho_{\mathrm{H}\mathrm{Z}}+ 1/\rho_{\Xi\mathrm{Z}}+ 1/\rho_{\Xi\mathrm{H}}- 1/\rho_{\Xi}- 1/\rho_{\mathrm{H}}- 1/\rho_{\mathrm{Z}} \right) \right\} \begin{bmatrix} \Xi\\ \mathrm{H}\\ \mathrm{Z}\\ \end{bmatrix} \]
The inverse mapping is analytic in the angular directions. A root find must be performed for the inverse mapping in the radial direction. This one-dimensional formula is obtained by taking the magnitude of both sides of the mapping, and changing variables from \(\xi, \eta, \zeta\) to \(x, y, z\) and introducing \(\rho^2 := \sqrt{\xi^2+\eta^2+\zeta^2}\).