SpECTRE
v2022.05.05

Map from a square or cube to a wedge. More...
#include <Wedge.hpp>
Public Types  
enum class  WedgeHalves { Both , UpperOnly , LowerOnly } 
Public Member Functions  
Wedge (double radius_inner, double radius_outer, double sphericity_inner, double sphericity_outer, OrientationMap< Dim > orientation_of_wedge, bool with_equiangular_map, WedgeHalves halves_to_use=WedgeHalves::Both, Distribution radial_distribution=Distribution::Linear)  
Wedge (Wedge &&)=default  
Wedge (const Wedge &)=default  
Wedge &  operator= (const Wedge &)=default 
Wedge &  operator= (Wedge &&)=default 
template<typename T >  
std::array< tt::remove_cvref_wrap_t< T >, Dim >  operator() (const std::array< T, Dim > &source_coords) const 
std::optional< std::array< double, Dim > >  inverse (const std::array< double, Dim > &target_coords) const 
For a \(+z\)oriented Wedge , returns invalid if \(z<=0\) or if \((x,y,z)\) is on or outside the cone defined by \((x^2/z^2 + y^2/z^2+1)^{1/2} = S/F\), where \(S = \frac{1}{2}(s_1 r_1  s_0 r_0)\) and \(F = \frac{1}{2\sqrt{3}}((1s_1) r_1  (1s_0) r_0)\). Here \(s_0,s_1\) and \(r_0,r_1\) are the specified sphericities and radii of the inner and outer \(z\) surfaces. The map is singular on the cone and on the xy plane. 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 >, Dim, Frame::NoFrame >  jacobian (const std::array< T, Dim > &source_coords) const 
template<typename T >  
tnsr::Ij< tt::remove_cvref_wrap_t< T >, Dim, Frame::NoFrame >  inv_jacobian (const std::array< T, Dim > &source_coords) const 
void  pup (PUP::er &p) 
Static Public Member Functions  
static constexpr bool  is_identity () 
Static Public Attributes  
static constexpr size_t  dim = Dim 
Friends  
template<size_t LocalDim>  
bool  operator== (const Wedge< LocalDim > &lhs, const Wedge< LocalDim > &rhs) 
Map from a square or cube to a wedge.
The mapping that goes from a reference cube (in 3D) or square (in 2D) to a wedge centered on a coordinate axis covering a volume between an inner surface and outer surface. Each surface can be given a curvature between flat (a sphericity of 0) or spherical (a sphericity of 1).
In 2D, the first logical coordinate corresponds to the radial coordinate, and the second logical coordinates correspond to the angular coordinate. In 3D, the first two logical coordinates correspond to the two angular coordinates, and the third to the radial coordinate. This difference originates from separate implementations for the 2D and 3D map that were merged. The 3D implementation can be changed to use the first logical coordinate as radial direction, but this requires propagating the change through the rest of the domain code (see issue https://github.com/sxscollaboration/spectre/issues/2988).
The following documentation is for the 3D map. The 2D map is obtained by setting either of the two angular coordinates to zero (and using \(\xi\) as radial coordinate).
The Wedge map is constructed by linearly interpolating between a bulged face of radius radius_of_inner_surface
to a bulged face of radius radius_of_outer_surface
, where the radius of each bulged face is defined to be the radius of the sphere circumscribing the bulge.
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:
\[\textrm{rho} : \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}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\]
Where
\[ \vec{x}(\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 bulged surface is itself constructed by linearly interpolating between a cubical face and a spherical face. The surface map for the cubical face of side length \(2L\) lying in the \(+z\) direction is given by:
\[\vec{\sigma}_{cubical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\]
Where
\[ \vec{x}(\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 this cubical face map and a spherical face map of radius \(R\), with the interpolation parameter being \(s\). The surface map for the bulged face lying in the \(+z\) direction is then given by:
\[\vec{\sigma}_{bulged}(\xi,\eta) = {(1s)L + \frac{sR}{\rho}} \begin{bmatrix} \Xi\\ \mathrm{H}\\ 1\\ \end{bmatrix}\]
We constrain L by demanding that the spherical face circumscribe the cube. With this condition, we have \(L = R/\sqrt3\).
The final map for the wedge which lies along the \(+z\) is obtained by interpolating between the two surfaces with the interpolation parameter being the logical coordinate \(\zeta\). This results in:
\[\vec{x}(\xi,\eta,\zeta) = \frac{1}{2}\left\{(1\zeta)\Big[(1s_{inner})\frac{R_{inner}}{\sqrt 3} + s_{inner}\frac{R_{inner}}{\rho}\Big] + (1+\zeta)\Big[(1s_{outer})\frac{R_{outer}}{\sqrt 3} +s_{outer} \frac{R_{outer}}{\rho}\Big] \right\}\begin{bmatrix} \Xi\\ \mathrm{H}\\ 1\\ \end{bmatrix}\]
We will define the variables \(F(\zeta)\) and \(S(\zeta)\), the frustum and sphere factors:
\[F(\zeta) = F_0 + F_1\zeta\]
\[S(\zeta) = S_0 + S_1\zeta\]
Where
\begin{align*}F_0 &= \frac{1}{2} \big\{ (1s_{outer})R_{outer} + (1s_{inner})R_{inner}\big\}\\ F_1 &= \partial_{\zeta} F = \frac{1}{2} \big\{ (1s_{outer})R_{outer}  (1s_{inner})R_{inner}\big\}\\ S_0 &= \frac{1}{2} \big\{ s_{outer}R_{outer} + s_{inner}R_{inner}\big\}\\ S_1 &= \partial_{\zeta} S = \frac{1}{2} \big\{ s_{outer}R_{outer}  s_{inner}R_{inner}\big\}\end{align*}
The map can then be rewritten as:
\[\vec{x}(\xi,\eta,\zeta) = \left\{\frac{F(\zeta)}{\sqrt 3} + \frac{S(\zeta)}{\rho}\right\}\begin{bmatrix} \Xi\\ \mathrm{H}\\ 1\\ \end{bmatrix}\]
We provide some common derivatives:
\[\partial_{\xi}z = \frac{S(\zeta)\Xi\Xi'}{\rho^3}\]
\[\partial_{\eta}z = \frac{S(\zeta)\mathrm{H}\mathrm{H}'}{\rho^3}\]
\[\partial_{\zeta}z = \frac{F'}{\sqrt 3} + \frac{S'}{\rho}\]
The Jacobian then is:
\[J = \begin{bmatrix} \Xi'z + \Xi\partial_{\xi}z & \Xi\partial_{\eta}z & \Xi\partial_{\zeta}z \\ \mathrm{H}\partial_{\xi}z & \mathrm{H}'z + \mathrm{H}\partial_{\eta}z & \mathrm{H}\partial_{\zeta}z\\ \partial_{\xi}z&\partial_{\eta}z &\partial_{\zeta}z \\ \end{bmatrix} \]
A common factor that shows up in the inverse jacobian is:
\[ T:= \frac{S(\zeta)}{(\partial_{\zeta}z)\rho^3}\]
The inverse Jacobian then is:
\[J^{1} = \frac{1}{z}\begin{bmatrix} \Xi'^{1} & 0 & \Xi\Xi'^{1}\\ 0 & \mathrm{H}'^{1} & \mathrm{H}\mathrm{H}'^{1}\\ T\Xi & T\mathrm{H} & T + F(\partial_{\zeta}z)^{1}/\sqrt 3\\ \end{bmatrix} \]
By default, Wedge linearly distributes its gridpoints in the radial direction. An exponential distribution of gridpoints can be obtained by linearly interpolating in the logarithm of the radius, in order to obtain a relatively higher resolution at smaller radii. Since this is a radial rescaling of Wedge, this option is only supported for fully spherical wedges with sphericity_inner
= sphericity_outer
= 1.
The linear interpolation done is:
\[ \log r = \frac{1\zeta}{2}\log R_{inner} + \frac{1+\zeta}{2}\log R_{outer} \]
The map then is:
\[\vec{x}(\xi,\eta,\zeta) = \frac{\sqrt{R_{inner}^{1\zeta}R_{outer}^{1+\zeta}}}{\rho}\begin{bmatrix} \Xi\\ \mathrm{H}\\ 1\\ \end{bmatrix}\]
The jacobian simplifies similarly.
Alternatively, an inverse radial distribution can be chosen where the linear interpolation is:
\[ \frac{1}{r} = \frac{R_\mathrm{inner} + R_\mathrm{outer}}{2 R_\mathrm{inner} R_\mathrm{outer}} + \frac{R_\mathrm{inner}  R_\mathrm{outer}}{2 R_\mathrm{inner} R_\mathrm{outer}} \zeta \]

strong 
domain::CoordinateMaps::Wedge< Dim >::Wedge  (  double  radius_inner, 
double  radius_outer,  
double  sphericity_inner,  
double  sphericity_outer,  
OrientationMap< Dim >  orientation_of_wedge,  
bool  with_equiangular_map,  
WedgeHalves  halves_to_use = WedgeHalves::Both , 

Distribution  radial_distribution = Distribution::Linear 

) 
Constructs a 3D wedge.
radius_inner  Distance from the origin to one of the corners which lie on the inner surface. 
radius_outer  Distance from the origin to one of the corners which lie on the outer surface. 
orientation_of_wedge  The orientation of the desired wedge relative to the orientation of the default wedge which is a wedge that has its curved surfaces pierced by the upperz axis. The logical xi and eta coordinates point in the cartesian x and y directions, respectively. 
sphericity_inner  Value between 0 and 1 which determines whether the inner surface is flat (value of 0), spherical (value of 1) or somewhere in between 
sphericity_outer  Value between 0 and 1 which determines whether the outer surface is flat (value of 0), spherical (value of 1) or somewhere in between 
with_equiangular_map  Determines whether to apply a tangent function mapping to the logical coordinates (for true ) or not (for false ). 
halves_to_use  Determines whether to construct a full wedge or only half a wedge. If constructing only half a wedge, the resulting shape has a face normal to the x direction (assuming default OrientationMap). If constructing half a wedge, an intermediate affine map is applied to the logical xi coordinate such that the interval [1,1] is mapped to the corresponding logical half of the wedge. For example, if UpperOnly is specified, [1,1] is mapped to [0,1], and if LowerOnly is specified, [1,1] is mapped to [1,0]. The case of Both means a full wedge, with no intermediate map applied. In all cases, the logical points returned by the inverse map will lie in the range [1,1] in each dimension. Half wedges are currently only useful in constructing domains for binary systems. 
radial_distribution  Determines how to distribute grid points along the radial direction. For wedges that are not exactly spherical, only Distribution::Linear is currently supported. 