Public Member Functions |
Static Public Member Functions |
Static Public Attributes |
Friends |
List of all members

domain::CoordinateMaps::FocallyLiftedMap< InnerMap > Class Template Reference

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) noexcept | |

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 noexcept |

std::optional< std::array< double, 3 > > | inverse (const std::array< double, 3 > &target_coords) const noexcept |

template<typename T > | |

tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > | jacobian (const std::array< T, 3 > &source_coords) const noexcept |

template<typename T > | |

tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > | inv_jacobian (const std::array< T, 3 > &source_coords) const noexcept |

void | pup (PUP::er &p) noexcept |

## Static Public Member Functions | |

static bool | is_identity () noexcept |

## Static Public Attributes | |

static constexpr size_t | dim = 3 |

## Friends | |

bool | operator== (const FocallyLiftedMap< InnerMap > &lhs, const FocallyLiftedMap< InnerMap > &rhs) noexcept |

class domain::CoordinateMaps::FocallyLiftedMap< InnerMap >

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 `InnerMap`

s function `deriv_lambda_tilde`

.

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

- src/Domain/CoordinateMaps/FocallyLiftedMap.hpp
- src/Domain/CoordinateMaps/FocallyLiftedMap.cpp

© Copyright 2017 - 2021 SXS Collaboration, Distributed under the MIT License