SpECTRE  v2021.12.06
domain::CoordinateMaps::BulgedCube Class Reference

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
 
BulgedCubeoperator= (const BulgedCube &)=default
 
BulgedCubeoperator= (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::NoFramejacobian (const std::array< T, 3 > &source_coords) const
 
template<typename T >
tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrameinv_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)
 

Detailed Description

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.

Details

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:

\(\Xi'(\xi) = 1\), and \(\mathrm{H}'(\eta) = 1\)

We also define the variable \(\rho\), given by:

\[\rho = \sqrt{1+\Xi^2+\mathrm{H}^2}\]

The Spherical Face Map

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 Cubical Face Map

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} \]

The Bulged Face Map

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 for 3D Isoparametric Maps

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 Special Case for Octahedral Symmetry

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:

\(\vec{\sigma}_{-\zeta}(\xi, \eta) = N_z\vec{\sigma}_{+\zeta}(\xi, \eta)\\ \vec{\sigma}_{+\eta}(\xi, \zeta) = S_{yz}\vec{\sigma}_{+\zeta}(\xi, \zeta)\\ \vec{\sigma}_{-\eta}(\xi, \zeta) = N_yS_{yz}\vec{\sigma}_{+\zeta}(\xi, \zeta)\\ \vec{\sigma}_{+\xi}(\eta, \zeta) = C_{zxy}\vec{\sigma}_{+\zeta}(\eta, \zeta)\\ \vec{\sigma}_{-\xi}(\eta, \zeta) = N_xC_{zyx}\vec{\sigma}_{+\zeta}(\eta, \zeta)\)

The four bounding curves \(\vec{\Gamma}\) on the \(+\zeta\) face are given by:

\(\vec{\Gamma}_{+\xi,+\zeta}(\eta) = \vec{\sigma}_{+\zeta}(+1,\eta)\\ \vec{\Gamma}_{-\xi,+\zeta}(\eta) = \vec{\sigma}_{+\zeta}(-1,\eta) = N_x\vec{\sigma}_{+\zeta}(+1, \eta)\\ \vec{\Gamma}_{+\eta,+\zeta}(\xi) = \vec{\sigma}_{+\zeta}(\xi,+1) = S_{xy}\vec{\sigma}_{+\zeta}(+1, \xi)\\ \vec{\Gamma}_{-\eta,+\zeta}(\xi) = \vec{\sigma}_{+\zeta}(\xi,-1) = N_yS_{xy}\vec{\sigma}_{+\zeta}(+1,\xi)\)

The bounding curves on the other surfaces can be obtained by transformations on the \(+\zeta\) face:

\(\vec{\Gamma}_{+\xi,-\zeta}(\eta) = N_z\vec{\sigma}_{+\zeta}(+1,\eta)\\ \vec{\Gamma}_{-\xi,-\zeta}(\eta) = N_z\vec{\sigma}_{+\zeta}(-1,\eta) = N_zN_x\vec{\sigma}_{+\zeta}(+1,\eta)\\ \vec{\Gamma}_{+\eta,-\zeta}(\xi) = N_z\vec{\sigma}_{+\zeta}(\xi,+1) = N_zS_{xy}\vec{\sigma}_{+\zeta}(+1, \xi)\\ \vec{\Gamma}_{-\eta,-\zeta}(\xi) = N_z\vec{\sigma}_{+\zeta}(\xi,-1) = N_zN_yS_{xy}\vec{\sigma}_{+\zeta}(+1, \xi)\\ \vec{\Gamma}_{+\xi,+\eta}(\zeta) = C_{zxy}\vec{\sigma}_{+\zeta}(+1,\zeta)\\ \vec{\Gamma}_{-\xi,+\eta}(\zeta) = N_xC_{zxy}\vec{\sigma}_{+\zeta}(+1,\zeta)\\ \vec{\Gamma}_{+\xi,-\eta}(\zeta) = C_{zxy}\vec{\sigma}_{+\zeta}(-1,\zeta) = C_{zxy}N_x\vec{\sigma}_{+\zeta}(+1,\zeta)\\ \vec{\Gamma}_{-\xi,-\eta}(\zeta) = N_xC_{zxy}\vec{\sigma}_{+\zeta}(-1,\zeta) = N_xC_{zxy}N_x\vec{\sigma}_{+\zeta}(+1,\zeta)\)

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*}

The Bulged Cube Map

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}\).


The documentation for this class was generated from the following file: