SpECTRE
v2023.05.16

Map from \((\bar{x},\bar{y},\bar{z})\) to the volume contained between a 2D surface and the surface of a 2sphere. More...
#include <FocallyLiftedMap.hpp>
Public Member Functions  
FocallyLiftedMap (const std::array< double, 3 > ¢er, const std::array< double, 3 > &proj_center, double radius, bool source_is_between_focus_and_target, InnerMap inner_map)  
FocallyLiftedMap (FocallyLiftedMap &&)=default  
FocallyLiftedMap (const FocallyLiftedMap &)=default  
FocallyLiftedMap &  operator= (const FocallyLiftedMap &)=default 
FocallyLiftedMap &  operator= (FocallyLiftedMap &&)=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) 
Static Public Member Functions  
static bool  is_identity () 
Static Public Attributes  
static constexpr size_t  dim = 3 
Friends  
bool  operator== (const FocallyLiftedMap< InnerMap > &lhs, const FocallyLiftedMap< InnerMap > &rhs) 
Map from \((\bar{x},\bar{y},\bar{z})\) to the volume contained between a 2D surface and the surface of a 2sphere.
We are given the radius \(R\) and center \(C^i\) of a sphere, and a projection point \(P^i\). Also, through the class defined by the InnerMap
template parameter, we are given the functions \(f^i(\bar{x},\bar{y},\bar{z})\) and \(\sigma(\bar{x},\bar{y},\bar{z})\) defined below; these functions define the mapping from \((\bar{x},\bar{y},\bar{z})\) to the 2D surface.
The above figure is a 2D representation of the focally lifted map, where the shaded region in \(\bar{x}^i\) coordinates on the left is mapped to the shaded region in the \(x^i\) coordinates on the right. Shown is an arbitrary point \(\bar{x}^i\) that is mapped to \(x^i\); for that point, the corresponding \(x_0^i\) and \(x_1^i\) (defined below) are shown on the right, as is the projection point \(P^i\). Also shown are the level surfaces \(\sigma=0\) and \(\sigma=1\).
The input coordinates are labeled \((\bar{x},\bar{y},\bar{z})\). Let \(x_0^i = f^i(\bar{x},\bar{y},\bar{z})\) be the three coordinates of the 2D surface in 3D space.
Now let \(x_1^i\) be a point on the surface of the sphere, constructed so that \(P^i\), \(x_0^i\), and \(x_1^i\) are colinear. In particular, \(x_1^i\) is determined by the equation
\begin{align} x_1^i = P^i + (x_0^i  P^i) \lambda,\end{align}
where \(\lambda\) is a scalar factor that depends on \(x_0^i\) and that can be computed by solving a quadratic equation. This quadratic equation is derived by demanding that \(x_1^i\) lies on the sphere:
\begin{align} P^i + (x_0^i  P^i) \lambda  C^i ^2  R^2 = 0, \end{align}
where \(A^i^2\) means \(\delta_{ij} A^i A^j\).
The quadratic equation, Eq. (2), takes the usual form \(a\lambda^2+b\lambda+c=0\), with
\begin{align*} a &= x_0^iP^i^2,\\ b &= 2(x_0^iP^i)(P^jC^j)\delta_{ij},\\ c &= P^iC^i^2  R^2. \end{align*}
Now assume that \(x_0^i\) lies on a level surface defined by some scalar function \(\sigma(\bar{x}^i)=0\), where \(\sigma\) is normalized so that \(\sigma=1\) on the sphere. Then, once \(\lambda\) has been computed and \(x_1^i\) has been determined, the map is given by
\begin{align}x^i = x_0^i + (x_1^i  x_0^i) \sigma(\bar{x}^i).\end{align}
Note that classes using FocallyLiftedMap will place restrictions on \(P^i\), \(C^i\), \(x_0^i\), and \(R\). For example, we demand that \(P^i\) does not lie on either the \(\sigma=0\) or \(\sigma=1\) surfaces depicted in the right panel of the figure, and we demand that the \(\sigma=0\) and \(\sigma=1\) surfaces do not intersect; otherwise the map is singular.
Also note that the quadratic Eq. (2) typically has more than one root, corresponding to two intersections of the sphere. The boolean parameter source_is_between_focus_and_target
that is passed into the constructor of FocallyLiftedMap
is used to choose the appropriate root, or to error if a suitable root is not found. source_is_between_focus_and_target
should be true if the source point lies between the projection point \(P^i\) and the sphere. source_is_between_focus_and_target
is known only by each particular CoordinateMap
that uses FocallyLiftedMap
.
Differentiating Eq. (1) above yields
\begin{align} \frac{\partial x_1^i}{\partial x_0^j} = \lambda \delta^i_j + (x_0^i  P^i) \frac{\partial \lambda}{\partial x_0^j}. \end{align}
and differentiating Eq. (2) and then inserting Eq. (1) yields
\begin{align} \frac{\partial\lambda}{\partial x_0^j} &= \lambda^2 \frac{C_j  x_1^j}{x_1^i  P^i^2 + (x_1^i  P^i)(P_i  C_{i})}. \end{align}
The Jacobian can be found by differentiating Eq. (3) above and using the chain rule, recognizing that \(x_1^i\) depends on \(\bar{x}^i\) only through its dependence on \(x_0^i\); this is because there is no explicit dependence on \(\bar{x}^i\) in Eq. (1) (which determines \(x_1^i\)) or Eq. (2) (which determines \(\lambda\)).
\begin{align} \frac{\partial x^i}{\partial \bar{x}^j} &= \sigma \frac{\partial x_1^i}{\partial x_0^k} \frac{\partial x_0^k}{\partial \bar{x}^j} + (1\sigma)\frac{\partial x_0^i}{\partial \bar{x}^j} + \frac{\partial \sigma}{\partial \bar{x}^j} (x_1^i  x_0^i), \nonumber \\ &= (1\sigma+\lambda\sigma) \frac{\partial x_0^i}{\partial \bar{x}^j} + \sigma (x_0^i  P^i) \frac{\partial \lambda}{\partial x_0^k} \frac{\partial x_0^k}{\partial \bar{x}^j} + \frac{\partial \sigma}{\partial \bar{x}^j} (x_1^i  x_0^i), \end{align}
where in the last line we have substituted Eq. (4).
The class defined by the InnerMap
template parameter provides the function deriv_sigma
, which returns \(\partial \sigma/\partial \bar{x}^j\), and the function jacobian
, which returns \(\partial x_0^k/\partial \bar{x}^j\).
Given \(x^i\), we wish to compute \(\bar{x}^i\).
We first find the coordinates \(x_0^i\) that lie on the 2surface and are defined such that \(P^i\), \(x_0^i\), and \(x^i\) are colinear. See the right panel of the above figure. \(x_0^i\) is determined by the equation
\begin{align} x_0^i = P^i + (x^i  P^i) \tilde{\lambda},\end{align}
where \(\tilde{\lambda}\) is a scalar factor that depends on \(x^i\) and is determined by the class defined by the InnerMap
template parameter. InnerMap
provides a function lambda_tilde
that takes \(x^i\) and \(P^i\) as arguments and returns \(\tilde{\lambda}\) (or a defaultconstructed std::optional
if the appropriate \(\tilde{\lambda}\) cannot be found; this defaultconstructed value indicates that the point \(x^i\) is outside the range of the map).
Now consider the coordinates \(x_1^i\) that lie on the sphere and are defined such that \(P^i\), \(x_1^i\), and \(x^i\) are colinear. See the right panel of the figure. \(x_1^i\) is determined by the equation
\begin{align} x_1^i = P^i + (x^i  P^i) \bar{\lambda},\end{align}
where \(\bar{\lambda}\) is a scalar factor that depends on \(x^i\) and is the solution of a quadratic that is derived by demanding that \(x_1^i\) lies on the sphere:
\begin{align} P^i + (x^i  P^i) \bar{\lambda}  C^i ^2  R^2 = 0. \end{align}
Eq. (9) is a quadratic equation that takes the usual form \(a\bar{\lambda}^2+b\bar{\lambda}+c=0\), with
\begin{align*} a &= x^iP^i^2,\\ b &= 2(x^iP^i)(P^jC^j)\delta_{ij},\\ c &= P^iC^i^2  R^2. \end{align*}
Note that we don't actually need to compute \(x_1^i\). Instead, we can determine \(\sigma\) by the relation (obtained by solving Eq. (3) for \(\sigma\) and then inserting Eqs. (7) and (8) to eliminate \(x_1^i\) and \(x_0^i\))
\begin{align} \sigma = \frac{\tilde{\lambda}1}{\tilde{\lambda}\bar{\lambda}}. \end{align}
The denominator of Eq. (10) is nonzero for nonsingular maps: From Eqs. (7) and (8), \(\bar{\lambda}=\tilde{\lambda}\) means that \(x_1^i=x_0^i\), which means that the \(\sigma=0\) and \(\sigma=1\) surfaces intersect, i.e. the map is singular.
Once we have \(x_0^i\) and \(\sigma\), the point \((\bar{x},\bar{y},\bar{z})\) is uniquely determined by InnerMap
. The inverse
function of InnerMap
takes \(x_0^i\) and \(\sigma\) as arguments, and returns \((\bar{x},\bar{y},\bar{z})\), or a defaultconstructed std::optional
if \(x_0^i\) or \(\sigma\) are outside the range of the map.
The inverse function described above will sometimes have errors that are noticeably larger than roundoff. Therefore we apply a single NewtonRaphson iteration to refine the result of the inverse map: Suppose we are given \(x^i\), and we have computed \(\bar{x}^i\) by the above procedure. We then correct \(\bar{x}^i\) by adding
\begin{align} \delta \bar{x}^i = \left(x^j  x^j(\bar{x})\right) \frac{\partial \bar{x}^i}{\partial x^j}, \end{align}
where \(x^j(\bar{x})\) is the result of applying the forward map to \(\bar{x}^i\) and \(\partial \bar{x}^i/\partial x^j\) is the inverse jacobian.
We write the inverse Jacobian as
\begin{align} \frac{\partial \bar{x}^i}{\partial x^j} = \frac{\partial \bar{x}^i}{\partial x_0^k} \frac{\partial x_0^k}{\partial x^j} + \frac{\partial \bar{x}^i}{\partial \sigma} \frac{\partial \sigma}{\partial x^j}, \end{align}
where we have recognized that \(\bar{x}^i\) depends both on \(x_0^k\) (the corresponding point on the 2surface) and on \(\sigma\) (encoding the distance away from the 2surface).
We now evaluate Eq. (12). The InnerMap
class provides a function inv_jacobian
that returns \(\partial \bar{x}^i/\partial x_0^k\) (where \(\sigma\) is held fixed), and a function dxbar_dsigma
that returns \(\partial \bar{x}^i/\partial \sigma\) (where \(x_0^i\) is held fixed). The factor \(\partial x_0^j/\partial x^i\) can be computed by differentiating Eq. (7), which yields
\begin{align} \frac{\partial x_0^j}{\partial x^i} &= \tilde{\lambda} \delta_i^j + \frac{x_0^jP^j}{\tilde{\lambda}} \frac{\partial\tilde{\lambda}}{\partial x^i}, \end{align}
where \(\partial \tilde{\lambda}/\partial x^i\) is provided by the deriv_lambda_tilde
function of InnerMap
. Note that for nonsingular maps there is no worry that \(\tilde{\lambda}\) is zero in the denominator of the second term of Eq. (13); if \(\tilde{\lambda}=0\) then \(x_0^i=P^i\) by Eq. (7), and therefore the map is singular.
To evaluate the remaining unknown factor in Eq. (12), \(\partial \sigma/\partial x^j\), note that [combining Eqs. (1), (7), and (8)] \(\bar{\lambda}=\tilde{\lambda}\lambda\). Therefore Eq. (10) is equivalent to
\begin{align} \sigma &= \frac{\tilde{\lambda}1}{\tilde{\lambda}(1\lambda)}. \end{align}
Differentiating this expression yields
\begin{align} \frac{\partial \sigma}{\partial x^i} &= \frac{\partial \sigma}{\partial \lambda} \frac{\partial \lambda}{\partial x_0^j} \frac{\partial x_0^j}{\partial x^i} + \frac{\partial \sigma}{\partial \tilde\lambda} \frac{\partial \tilde\lambda}{\partial x^i}\\ &= \frac{\sigma}{1\lambda} \frac{\partial \lambda}{\partial x_0^j} \frac{\partial x_0^j}{\partial x^i} + \frac{1}{\tilde{\lambda}^2(1\lambda)} \frac{\partial \tilde\lambda}{\partial x^i}, \end{align}
where the second factor in the first term can be evaluated using Eq. (5), the third factor in the first term can be evaluated using Eq. (13), and the second factor in the second term is provided by InnerMap
s function deriv_lambda_tilde
.