|
SpECTRE
v2024.04.12
|
Map from \((\bar{x},\bar{y},\bar{z})\) to the volume contained between a 2D surface and the surface of a 2-sphere. 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 2-sphere.
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 co-linear. 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^i-P^i|^2,\\ b &= 2(x_0^i-P^i)(P^j-C^j)\delta_{ij},\\ c &= |P^i-C^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 2-surface and are defined such that \(P^i\), \(x_0^i\), and \(x^i\) are co-linear. 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 default-constructed std::optional if the appropriate \(\tilde{\lambda}\) cannot be found; this default-constructed 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 co-linear. 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^i-P^i|^2,\\ b &= 2(x^i-P^i)(P^j-C^j)\delta_{ij},\\ c &= |P^i-C^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 default-constructed 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 Newton-Raphson 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 2-surface) and on \(\sigma\) (encoding the distance away from the 2-surface).
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^j-P^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 InnerMaps function deriv_lambda_tilde.