SpECTRE
v2024.09.29

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, const std::array< double, Dim  1 > &opening_angles=make_array< Dim  1 >(M_PI_2), bool with_adapted_equiangular_map=true)  
Constructs a centered wedge (one with no focal offset) More...  
Wedge (double radius_inner, std::optional< double > radius_outer, double cube_half_length, std::array< double, Dim > focal_offset, OrientationMap< Dim > orientation_of_wedge, bool with_equiangular_map, WedgeHalves halves_to_use=WedgeHalves::Both, Distribution radial_distribution=Distribution::Linear)  
Constructs a wedge with a focal offset. More...  
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 coordinate corresponds 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 the 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 centered 3D map, as we will defer the dicussion of Wedge
s with a focal_offset_
to a later section. The 2D map is obtained by setting either of the two angular coordinates to zero (and using \(\xi\) as the radial coordinate). Note that there is also a normalization factor of \(\sqrt{3}\) that appears in multiple expressions in the 3D case that becomes \(\sqrt{2}\) in the 2D case.
The Wedge map is constructed by linearly interpolating between a bulged face of radius radius_inner_
to a bulged face of radius radius_outer_
, 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. Wedge
s have variable opening_angles_
which, for centered Wedge
s, are the angular sizes of the wedge in the \(\xi\) and \(\eta\) directions (for the 3D case) in the target frame. By default, Wedge
s have opening angles of \(\pi/2\), so we will discuss that case here and defer both the discussion of generalized opening angles and the interaction between opening angles and nonzero focal offsets for later sections.
For a Wedge with \(\xi\) and \(\eta\) opening angles of \(\pi/2\), the equiangular coordinates in terms of the logical coordinates are:
\begin{align} \textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4) \label{eq:equiangular_xi_pi_over_2} \end{align}
\begin{align} \textrm{equiangular eta} : \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4) \label{eq:equiangular_eta_pi_over_2} \end{align}
With derivatives:
\begin{align} \Xi'(\xi) &= \frac{\pi}{4}(1+\Xi^2) \\ \mathrm{H}'(\eta) &= \frac{\pi}{4}(1+\mathrm{H}^2) \end{align}
The equidistant coordinates are:
\begin{align} \textrm{equidistant xi} : \Xi = \xi \\ \textrm{equidistant eta} : \mathrm{H} = \eta \end{align}
with derivatives:
\begin{align} \Xi'(\xi) &= 1 \\ \mathrm{H}'(\eta) &= 1 \end{align}
We also define the variable \(\rho\), given by:
\begin{align} \textrm{rho} : \rho = \sqrt{1+\Xi^2+\mathrm{H}^2} \end{align}
The surface map for the spherical face of radius \(R\) lying in the \(+z\) direction in either choice of coordinates is then given by:
\begin{align} \vec{\sigma}_{spherical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi}) \end{align}
Where
\begin{align} \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} \end{align}
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:
\begin{align} \vec{\sigma}_{cubical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi}) \end{align}
Where
\begin{align} \vec{x}(\xi,\eta) = \begin{bmatrix} x(\xi,\eta) \\ y(\xi,\eta) \\ L \\ \end{bmatrix} = L\begin{bmatrix} \Xi \\ \mathrm{H} \\ 1 \\ \end{bmatrix} \end{align}
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\), called the sphericity and which ranges from 0 to 1, with 0 corresponding to a flat surface and 1 corresponding to a spherical surface. The surface map for the bulged face lying in the \(+z\) direction is then given by:
\begin{align} \vec{\sigma}_{bulged}(\xi,\eta) = \left\{(1s)L + \frac{sR}{\rho}\right\} \begin{bmatrix} \Xi \\ \mathrm{H} \\ 1 \\ \end{bmatrix} \end{align}
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\) axis is obtained by interpolating between the two surfaces with the interpolation parameter being the logical coordinate \(\zeta\). For a wedge whose gridpoints are linearly distributed in the radial direction (radial_distribution_
is domain::CoordinateMaps::Distribution::Linear), this interpolation results in the following map:
\begin{align} \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} \end{align}
We will define the variables \(F(\zeta)\) and \(S(\zeta)\), the frustum and sphere factors (in the linear case):
\begin{align} F(\zeta) &= F_0 + F_1\zeta \label{eq:frustum_factor} \\ S(\zeta) &= S_0 + S_1\zeta \label{eq:sphere_factor} \end{align}
Where
\begin{align} F_0 &= \frac{1}{2} \big\{ (1s_{outer})R_{outer} + (1s_{inner})R_{inner} \big\} \label{eq:frustum_zero_linear} \\ F_1 &= \partial_{\zeta}F = \frac{1}{2} \big\{ (1s_{outer})R_{outer}  (1s_{inner})R_{inner} \big\} \label{eq:frustum_rate_linear} \\ S_0 &= \frac{1}{2} \big\{ s_{outer}R_{outer} + s_{inner}R_{inner} \big\} \label{eq:sphere_zero_linear} \\ S_1 &= \partial_{\zeta}S = \frac{1}{2} \big\{ s_{outer}R_{outer}  s_{inner}R_{inner}\big\} \label{eq:sphere_rate_linear} \end{align}
The map can then be rewritten as:
\begin{align} \vec{x}(\xi,\eta,\zeta) = \left\{ \frac{F(\zeta)}{\sqrt 3} + \frac{S(\zeta)}{\rho} \right\} \begin{bmatrix} \Xi \\ \mathrm{H} \\ 1 \\ \end{bmatrix} \end{align}
The inverse map is given by:
\begin{align} \xi &= \frac{x}{z} \\ \eta &= \frac{y}{z} \\ \zeta &= \frac{z  \left(\frac{F_0}{\sqrt{3}} + \frac{S_0}{\rho}\right)} {\left(\frac{F_1}{\sqrt{3}} + \frac{S_1}{\rho}\right)} \end{align}
We provide some common derivatives:
\begin{align} \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'(\zeta)}{\rho} \end{align}
The Jacobian then is:
\begin{align} 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} \label{eq:jacobian_centered_wedge} \end{align}
A common factor that shows up in the inverse Jacobian is:
\begin{align} T:= \frac{S(\zeta)}{(\partial_{\zeta}z)\rho^3} \end{align}
The inverse Jacobian then is:
\begin{align} 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} \end{align}
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 for a logarithmic radial distribution (radial_distribution_
is domain::CoordinateMaps::Distribution::Logarithmic) is:
\begin{align} \ln r = \frac{1\zeta}{2}\ln R_{inner} + \frac{1+\zeta}{2}\ln R_{outer} \end{align}
The map then is:
\begin{align} \vec{x}(\xi,\eta,\zeta) = \frac{\sqrt{R_{inner}^{1\zeta}R_{outer}^{1+\zeta}}}{\rho} \begin{bmatrix} \Xi \\ \mathrm{H} \\ 1 \\ \end{bmatrix} \end{align}
We can rewrite this map to take on the same form as the map for the linear radial distribution, where we set
\begin{align} F(\zeta) &= 0 \\ S(\zeta) &= \sqrt{R_{inner}^{1\zeta}R_{outer}^{1+\zeta}} \\ \end{align}
Which gives us
\begin{align} \vec{x}(\xi,\eta,\zeta) = \frac{S(\zeta)}{\rho} \begin{bmatrix} \Xi \\ \mathrm{H} \\ 1 \\ \end{bmatrix} \end{align}
The Jacobian then is still Eq. ( \(\ref{eq:jacobian_centered_wedge}\)) but where \(F(\zeta)\) and \(S(\zeta)\) are the quantities defined here for the logarithmic distribution.
Alternatively, an inverse radial distribution (radial_distribution_
is domain::CoordinateMaps::Distribution::Inverse) can be chosen where the linear interpolation is:
\begin{align} \frac{1}{r} = \frac{R_\mathrm{inner} + R_\mathrm{outer}} {2 R_\mathrm{inner}R_\mathrm{outer}} + \frac{R_\mathrm{inner}  R_\mathrm{outer}} {2R_\mathrm{inner} R_\mathrm{outer}} \zeta \end{align}
Which can be rewritten as:
\begin{align} \frac{1}{r} = \frac{1\zeta}{2R_{inner}} + \frac{1+\zeta}{2R_{outer}} \end{align}
The map likewise takes the form:
\begin{align} \vec{x}(\xi,\eta,\zeta) = \frac{S(\zeta)}{\rho} \begin{bmatrix} \Xi \\ \mathrm{H} \\ 1 \\ \end{bmatrix} \end{align}
Where
\begin{align} F(\zeta) &= 0 \\ S(\zeta) &= \frac{2R_{inner}R_{outer}} {(1 + \zeta)R_{inner} + (1  \zeta)R_{outer}} \end{align}
Again, the Jacobian is still Eq. ( \(\ref{eq:jacobian_centered_wedge}\)) but where \(F(\zeta)\) and \(S(\zeta)\) are the quantities defined here for the inverse distribution.
Consider the following map on \(\xi \in [1,1]\), which maps this interval onto a parameterized curve that extends one fourth of a circle.
\begin{align} \vec{\Gamma}(\xi) = \frac{R}{\sqrt{1+\xi^2}} \begin{bmatrix} 1 \\ \xi \\ \end{bmatrix}. \label{eq:quarter_circle} \end{align}
It is convenient to compute the polar coordinate \(\theta\) of the mapped point as a function of \(\xi\):
\begin{align} \theta(\xi) = \tan^{1}\left(\frac{\Gamma_y(\xi)}{\Gamma_x(\xi)}\right). \label{eq:polar_coord} \end{align}
The opening angle of the map is defined to be:
\begin{align} \Delta \theta = \theta(1)  \theta(1), \label{eq:define_opening_angle} \end{align}
We can see that with \(\xi=\pm 1\), we have \(\Gamma_x = R/\sqrt{2}\) and \(\Gamma_y=\pm R/\sqrt{2}\), giving us \(\theta(1) = \pi/4\) and \(\theta(1) = \pi/4\). This wedge has an opening angle \(\pi/2\) radians, as expected.
On the other hand, the following map has an opening angle of \(\theta_O\):
\begin{align} \vec{\Gamma}(\xi) = \frac{R}{\sqrt{1+\tan^2{(\theta_O/2)}\xi^2}} \begin{bmatrix} 1 \\ \tan{(\theta_O/2)}\xi \\ \end{bmatrix}. \end{align}
Let us also consider the generalized map
\begin{align} \vec{\Gamma}(\xi) = \frac{R}{\sqrt{1+\Xi^2}} \begin{bmatrix} 1 \\ \Xi \\ \end{bmatrix}, \end{align}
where \(\Xi(\xi)\) is a function of \(\xi\). \(\theta(\xi)\) can then be written as
\begin{align} \theta(\xi) = \tan^{1}(\Xi). \label{eq:theta} \end{align}
For the map \(\Xi(\xi) = \tan(\pi\xi/4)\), Eq. ( \(\ref{eq:theta}\)) yields \(\theta(\xi) = \pi\xi/4\) and \(\Delta\theta = \pi/2\). Note that this choice of \(\Xi(\xi)\) is equivalent to a reparameterization of the previous map given in Eq. ( \(\ref{eq:quarter_circle}\)). The reparameterization of the curve \(\vec{\Gamma}(\xi)\) via the tangent map yields an empirically superior gridpoint distribution in practice. That this reparameterization should have this property can be motivated by an observation of the following:
\begin{align} \frac{\mathrm{d}\tan^{1}\Xi}{\mathrm{d}\xi} = \frac{1}{1+\Xi^2}\frac{\mathrm{d}\Xi}{\mathrm{d}\xi} = \frac{\pi}{4}. \end{align}
In other words, this parameterization has the property that the logical coordinate \(\xi\) subtends the angle \(\theta\) at a constant rate. In general, we say that a curve \(\vec{\Gamma}(\xi)\) is parameterized equiangularly if
\begin{align} \frac{\mathrm{d}\theta}{\mathrm{d}\xi} = \text{const}. \end{align}
As for the map
\begin{align} \Xi(\xi) = \tan{(\theta_O/2)}\frac{\tan{(\theta_D \xi/2)}}{\tan{(\theta_D/2)}}, \end{align}
this choice of \(\Xi(\xi)\) results in a \(\vec{\Gamma}(\xi)\) with opening angle \(\theta_O\), which is equiangularly distributed if \(\theta_O = \theta_D\). In the Wedge map, the argument with_adapted_equiangular_map
controls whether to set \(\theta_O = \theta_D\) (the true
case) or to set \(\theta_D = \pi/2\) (the false
case). When working with a 3D Wedge, the opening angles for the Wedge can be separately controlled for both the \(\xi\) and \(\eta\) directions, but with_adapted_equiangular_map
will apply to both directions. Additionally in the 3D case, it is not possible to set with_equiangular_map_
to true
for all of the six wedges of a sphere unless every opening angle is \(\pi/2\). In the BinaryCompactObject domain, the outer \(+y\), \(y\), \(+z\), and \(z\) Wedge
s are allowed to have a userspecified opening angle in the \(\xi\)direction, with a corresponding \(\theta_D\) equal to this opening angle, while in the \(\eta\)direction the opening angle is set to \(\pi/2\). The two end cap Wedge
s in the \(+x\) and \(x\) directions have angular dimensions and gridpoint distributions determined by the other four Wedge
s, as the six Wedge
s must conforming have gridpoint distributions at the \(\xi = \pm1\), \(\eta = \pm 1\) boundaries.
In the case of the rectangular BinaryCompactObject domain, it becomes desirable to offset the center of the spherical excision surface relative to the center of the cubical surface surrounding it. To enable the offsetting of the central excision, the Wedge map must be generalized according to the focal lifting method, which we will now discuss.
We consider the problem of creating parameterized volumes from parameterized surfaces. Consider a parameterized surface \(\vec{\sigma}_{parent}(\xi,\eta)\), also referred to as the parent surface. We define focal lifting as the projection of this parent surface into a threedimensional parameterized volume \(\vec{x}(\xi,\eta, \zeta)\) with respect to some focus \(\vec{x}_0\) and lifting scale factor \(\Lambda(\xi,\eta,\zeta)\). The resulting volume is then said to be a focally lifted volume. These volume maps can be cast into the following form:
\begin{align} \vec{x}  \vec{x}_0 = \Lambda(\vec{\sigma}_{parent}\vec{x}_0), \label{eq:focal_lifting} \end{align}
which makes apparent how the mapped point \(\vec{x}(\xi,\eta,\zeta)\) is obtained. The parametric equations for the generalized 3D Wedge maps can all be written in the above form, which we will refer to as focally lifted form. In the case of the 3D Wedge map with no focal offset, we have:
\begin{align} \vec{x}_0 &= 0 \\ \Lambda &= \left\{\frac{F(\zeta)}{\sqrt{3}} + \frac{S(\zeta)}{\rho} \right\} \\ \vec{\sigma}_{parent} &= \begin{bmatrix} \Xi, \mathrm{H}, 1 \end{bmatrix}^T \end{align}
The above map can be thought of as constructing a wedge from a biunit cube centered at the origin. Points on the parent surface are scaled by a factor of \(\Lambda(\xi,\eta,\zeta)\) to obtain the corresponding point in the volume. When generalizing the map to have a focus shifted from the origin (obtained by setting focal_offset_
to be nonzero), we scale the original parent surface \(\vec{\sigma}_{parent} = [\Xi, \mathrm{H},1]^T\) by a factor \(L\), and let the focus \(\vec{x_0}\) shift away from the origin. The generalized wedge map is then given by:
\begin{align} \vec{x}  \vec{x}_0 = \left\{\frac{F(\zeta)}{L\sqrt 3} + \frac{S(\zeta)}{L\rho}\right\} \begin{bmatrix} L\Xi  x_0 \\ L\mathrm{H}  y_0 \\ Lz_0 \\ \end{bmatrix} \end{align}
where we are now defining \(\rho\) to be
\begin{align} \rho = \sqrt{(\Xi  x_0/L)^2 + (\mathrm{H}  y_0/L)^2 + (1  z_0/L)^2}. \label{eq:generalized_rho} \end{align}
This map is often written as:
\begin{align} \vec{x}  \vec{x}_0 = \left\{\frac{F(\zeta)}{\sqrt{3}} + \frac{S(\zeta)}{\rho}\right\}(\vec{\sigma}_0  \vec{x}_0/L), \label{eq:focally_lifted_map_with_s_and_f_factors} \end{align}
where \(\vec{\sigma}_0 = [\Xi, \mathrm{H},1]^T\), as the parent surface \(\vec{\sigma}_{parent}\) is now \(L\vec{\sigma}_0\). We give the quantity in braces the name \(z_{\Lambda} = L\Lambda\), generalized z. With this definition, we can rewrite Eq. ( \(\ref{eq:focally_lifted_map_with_s_and_f_factors}\)) in the simpler form,
\begin{align} \vec{x}  \vec{x}_0 = z_{\Lambda}(\vec{\sigma}_0  \vec{x}_0/L). \label{eq:focally_lifted_map_with_generalized_z_coef} \end{align}
The map can be inverted by first solving for \(z_{\Lambda}\) in terms of the target coordinates. We make use of the fact that the parent surface \(\vec{\sigma}_{parent}\) has a constant normal vector \(\hat{n} = \hat{z}\).
\begin{align} z_{\Lambda} = \frac{(\vec{x}  \vec{x}_0)\cdot\hat{n}} {(\vec{\sigma}_0\vec{x}_0/L)\cdot\hat{n}}. \end{align}
In other words, when \(\hat{n} = \hat{z}\),
\begin{align} z_{\Lambda} = \left\{\frac{F(\zeta)}{\sqrt{3}} + \frac{S(\zeta)}{\rho}\right\} = \frac{z  z_0}{1  z_0/L} \end{align}
Moving all the known quantities in Eq. ( \(\ref{eq:focally_lifted_map_with_generalized_z_coef}\)) to the left hand side results in the following expression that solves for the source coordinates \(\xi\) and \(\eta\) in terms of the target coordinates:
\begin{align} \frac{\vec{x}  \vec{x}_0}{z_{\Lambda}} + \frac{\vec{x}_0}{L} = \vec{\sigma}_0(\xi,\eta) = \begin{bmatrix} \Xi \\ \mathrm{H} \\ 1 \\ \end{bmatrix}, \end{align}
Note that \(\vec{\sigma}_0  \vec{x}_0/L = \sqrt{(\Xi  x_0/L)^2 + (\mathrm{H}  y_0/L)^2 + (1  z_0/L)^2} = \rho\), indicating that an expression for \(\rho\) in terms of the target coordinates can be computed via taking the magnitude of both sides of Eq. ( \(\ref{eq:focally_lifted_map_with_generalized_z_coef}\)):
\begin{align} \vec{x}  \vec{x}_0 = z_{\Lambda}\vec{\sigma}_0  \vec{x}_0/L = z_{\Lambda}\rho. \end{align}
The quantity \(\rho\) is then given by:
\begin{align} \rho = \frac{\vec{x}  \vec{x}_0}{z_{\Lambda}}. \end{align}
With \(\rho\) computed, the radial source coordinate \(\zeta\) can be computed from
\begin{align} z_{\Lambda} = \left\{\frac{F(\zeta)}{\sqrt{3}} + \frac{S(\zeta)}{\rho} \right\} = \left\{\frac{F_0}{\sqrt{3}} + \frac{S_0}{\rho} + \frac{F_1\zeta}{\sqrt{3}} + \frac{S_1\zeta}{\rho}\right\}, \end{align}
which gives
\begin{align} \zeta = \frac{z_{\Lambda}  \left(\frac{F_0}{\sqrt{3}} + \frac{S_0}{\rho}\right)} {\left(\frac{F_1}{\sqrt{3}} + \frac{S_1}{\rho}\right)}. \end{align}
To compute the Jacobian, it is useful to first note that \(\rho\) (Eq. ( \(\ref{eq:generalized_rho}\))) is the magnitude of the vector
\begin{align} \vec{\rho} = \vec{\sigma}_0  \vec{x}_0/L = \begin{bmatrix} \Xi  x_0/L \\ \mathrm{H}  y_0/L \\ 1  z_0/L \end{bmatrix} \end{align}
and that we can express the target coordinates in Eq. ( \(\ref{eq:focally_lifted_map_with_generalized_z_coef}\)) in terms of the components of \(\vec{\rho}\):
\begin{align} x &= z_{\Lambda}\rho_x + x_0 \\ y &= z_{\Lambda}\rho_y + y_0 \\ z &= z_{\Lambda}\rho_z + z_0 \end{align}
Some common terms used in the Jacobian are the derivatives of \(z_{\Lambda}\) with respect to the source coordinates:
\begin{align} \partial_{\xi}z_{\Lambda} &= \frac{S(\zeta)\Xi'\rho_x}{\rho^3} \\ \partial_{\eta}z_{\Lambda} &= \frac{S(\zeta)\mathrm{H}'\rho_y}{\rho^3} \\ \partial_{\zeta}z_{\Lambda} &= \frac{F'(\zeta)}{\sqrt{3}} + \frac{S'(\zeta)}{\rho} \end{align}
The Jacobian then is:
\begin{align} J = \begin{bmatrix} \Xi'z_{\Lambda} + \rho_x\partial_{\xi}z_{\Lambda} & \rho_x\partial_{\eta}z_{\Lambda} & \rho_x\partial_{\zeta}z_{\Lambda} \\ \rho_y\partial_{\xi}z_{\Lambda} & \mathrm{H}'z_{\Lambda} + \rho_y\partial_{\eta}z_{\Lambda} & \rho_y\partial_{\zeta}z_{\Lambda} \\ \rho_z\partial_{\xi}z_{\Lambda} & \rho_z\partial_{\eta}z_{\Lambda} & \rho_z\partial_{\zeta}z_{\Lambda} \\ \end{bmatrix} \end{align}
A common factor that shows up in this inverse Jacobian is:
\begin{align} T:= \frac{S(\zeta)}{(\partial_{\zeta}z_{\Lambda})\rho^3} \end{align}
And the inverse Jacobian is then:
\begin{align} J^{1} = \frac{1}{z_{\Lambda}}\begin{bmatrix} \Xi'^{1} & 0 & \rho_x(\Xi'\rho_z)^{1} \\ 0 & \mathrm{H}'^{1} & \rho_y(\mathrm{H}'\rho_z)^{1} \\ T\rho_x & T\rho_y & T\rho_z + F(\partial_{\zeta}z_{\Lambda}\rho_z)^{1}/\sqrt{3} \end{bmatrix} \end{align}
The default Wedge map is oriented in the \(+z\) direction, so the construction of a Wedge oriented along a different direction requires an additional OrientationMap \(R\) to be passed to orientation_of_wedge
. When offsetting a rotated Wedge, the coordinates passed as parameters to focal_offset
are in the coordinate frame in which the Wedge is rotated (the target frame). However, the focal lifting procedure (shown in Eq. ( \(\ref{eq:focal_lifting}\))) is done in the default frame in which the Wedge is facing the \(+z\) direction, so the focal offset \(\vec{x}_0\) is first hit by the inverse rotation \(R^{1}\) and then the rotated focus \(R^{1}\vec{x}_0\) is used internally as the focus for the \(+z\) Wedge. When the focal lifting calculation has completed, the rotation of the \(+z\) Wedge into the desired orientation by \(R\) also rotates the focus into the desired location. When performing the inverse operation, the focus is similarly rotated into the default frame, where the inversion is performed.
When a Wedge is created with a nonzero focal offset, the resulting shape can take on a variety of possible angular sizes, depending on where the focus is placed relative to the default centered location. The reader might note that the angular size of a Wedge can also be controlled by passing an argument to the opening_angles
parameter in the Wedge constructor. While both of these methods allow the angular size of a Wedge to be changed, the user is prevented from employing both of them at the same time. In particular, when the the offset is set to some nonzero value, the opening_angles_
member variable is set to \(\pi/2\). Note that the opening_angles_
member being set to \(\pi/2\) does not imply the resulting Wedge will have an angular size of \(\pi/2\). On the contrary, the Wedge will have the angular size that is determined by the application of the focal lifting method on the parent surface, which is the upper \(+z\) face of a cube that is centered at the origin.
Because opening_angles_
is set to \(\pi/2\) when there is a nonzero focal offset, when there is a nonzero focal offset and with_equiangular_map_
is true
, \(\Xi\) is given by Eq. ( \(\ref{eq:equiangular_xi_pi_over_2}\)) and \(\mathrm{H}\) by Eq. ( \(\ref{eq:equiangular_eta_pi_over_2}\)), just as it is for the case of a centered Wedge with opening_angles_
of \(\pi/2\).

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 , 

const std::array< double, Dim  1 > &  opening_angles = make_array< Dim  1 >(M_PI_2) , 

bool  with_adapted_equiangular_map = true 

) 
Constructs a centered wedge (one with no focal offset)
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 gridpoints along the radial direction. For wedges that are not exactly spherical, only Distribution::Linear is currently supported. 
opening_angles  Determines the angular size of the wedge. The default value is \(\pi/2\), which corresponds to a wedge size of \(\pi/2\). For this setting, four Wedges can be put together to cover \(2\pi\) in angle along a great circle. This option is meant to be used with the equiangular map option turned on. 
with_adapted_equiangular_map  Determines whether to adapt the point distribution in the wedge to match its physical angular size. When true , angular distances are proportional to logical distances. Note that it is not possible to use adapted maps in every Wedge of a Sphere unless each Wedge has the same size along both angular directions. 
domain::CoordinateMaps::Wedge< Dim >::Wedge  (  double  radius_inner, 
std::optional< double >  radius_outer,  
double  cube_half_length,  
std::array< double, Dim >  focal_offset,  
OrientationMap< Dim >  orientation_of_wedge,  
bool  with_equiangular_map,  
WedgeHalves  halves_to_use = WedgeHalves::Both , 

Distribution  radial_distribution = Distribution::Linear 

) 
Constructs a wedge with a focal offset.
Can construct an offset Wedge with a spherical inner surface and either a spherical or a flat outer surface. If radius_outer
has a value, a spherical Wedge will be constructed, and if not, a flat one will be constructed.
Note that because the focal offset is what determines the angular size of the Wedge, opening angles cannot be used with offset Wedges.
In the event that focal_offset
happens to be zero, the Wedge's member variables and behavior will be set up to be equivalent to that of a centered Wedge:
cube_half_length
will be ignoredradius_outer
is std::nullopt
, the outer radius of the Wedge will be set to \(\sqrt{\mathrm{Dim}}L\), where \(L\) is the cube_half_length
radius_inner  Distance from the origin to one of the corners which lie on the inner surface. 
radius_outer  If this has a value, it creates a spherical Wedge (inner and outer sphericity are 1) where this is the distance from the origin to one of the corners that lie on the inner surface. If this is std::nullopt , it creates a Wedge with a flat outer surface (inner sphericity is 1 and outer sphericity is 0). In the event that radius_outer == std::nullopt and focal_offset is zero, the outer radius will instead be set to \(\sqrt{\mathrm{Dim}}L\), where \(L\) is the cube_half_length . The outer radius is given a value in this circumstance so that it can be handled as a centered Wedge (one with no offset). 
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. 
cube_half_length  Half the length of the parent surface (see Wedge documentation for more details). If focal_offset is zero, this parameter has no effect and is ignored so that the Wedge can be handled as a centered Wedge (one with no offset). 
focal_offset  The target frame coordinates of the focus from which the Wedge is focally lifted. 
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 gridpoints along the radial direction. For wedges that are not exactly spherical, only Distribution::Linear is currently supported. 