SpECTRE  v2022.10.04
domain::CoordinateMaps::FocallyLiftedInnerMaps::Endcap Class Reference

A FocallyLiftedInnerMap that maps a 3D unit right cylinder to a volume that connects portions of two spherical surfaces. More...

#include <FocallyLiftedEndcap.hpp>

Public Member Functions

Endcap (const std::array< double, 3 > &center, double radius, double z_plane)

Endcap (Endcap &&)=default

Endcap (const Endcap &)=default

Endcapoperator= (const Endcap &)=default

Endcapoperator= (Endcap &&)=default

template<typename T >
void forward_map (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > target_coords, const std::array< T, 3 > &source_coords) const

std::optional< std::array< double, 3 > > inverse (const std::array< double, 3 > &target_coords, double sigma_in) 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 >
void jacobian (const gsl::not_null< tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > * > jacobian_out, const std::array< T, 3 > &source_coords) const

template<typename T >
void inv_jacobian (const gsl::not_null< tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > * > inv_jacobian_out, const std::array< T, 3 > &source_coords) const

template<typename T >
void sigma (const gsl::not_null< tt::remove_cvref_wrap_t< T > * > sigma_out, const std::array< T, 3 > &source_coords) const

template<typename T >
void deriv_sigma (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > deriv_sigma_out, const std::array< T, 3 > &source_coords) const

template<typename T >
void dxbar_dsigma (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > dxbar_dsigma_out, const std::array< T, 3 > &source_coords) const

std::optional< double > lambda_tilde (const std::array< double, 3 > &parent_mapped_target_coords, const std::array< double, 3 > &projection_point, bool source_is_between_focus_and_target) const

template<typename T >
void deriv_lambda_tilde (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > deriv_lambda_tilde_out, const std::array< T, 3 > &target_coords, const T &lambda_tilde, const std::array< double, 3 > &projection_point) const

void pup (PUP::er &p)

Static Public Member Functions

static bool is_identity ()

Friends

bool operator== (const Endcap &lhs, const Endcap &rhs)

Detailed Description

A FocallyLiftedInnerMap that maps a 3D unit right cylinder to a volume that connects portions of two spherical surfaces.

Because FocallyLiftedEndcap is a FocallyLiftedInnerMap, it is meant to be a template parameter of FocallyLiftedMap, and its member functions are meant to be used by FocallyLiftedMap. See FocallyLiftedMap for further documentation.

Focally Lifted Endcap.

Details

The domain of the map is a 3D unit right cylinder with coordinates $$(\bar{x},\bar{y},\bar{z})$$ such that $$-1\leq\bar{z}\leq 1$$ and $$\bar{x}^2+\bar{y}^2 \leq 1$$. The range of the map has coordinates $$(x,y,z)$$.

Consider a sphere with center $$C^i$$ and radius $$R$$ that is intersected by a plane normal to the $$z$$ axis and located at $$z = z_\mathrm{P}$$. In the figure above, every point $$\bar{x}^i$$ in the blue region $$\sigma=0$$ maps to a point $$x_0^i$$ on a portion of the surface of the sphere.

Endcap provides the following functions:

forward_map

forward_map maps $$(\bar{x},\bar{y},\bar{z}=-1)$$ to the portion of the sphere with $$z \geq z_\mathrm{P}$$. The arguments to forward_map are $$(\bar{x},\bar{y},\bar{z})$$, but $$\bar{z}$$ is ignored. forward_map returns $$x_0^i$$, the 3D coordinates on that sphere, which are given by

\begin{align} x_0^0 &= R \frac{\sin(\bar{\rho} \theta_\mathrm{max}) \bar{x}}{\bar{\rho}} + C^0,\\ x_0^1 &= R \frac{\sin(\bar{\rho} \theta_\mathrm{max}) \bar{y}}{\bar{\rho}} + C^1,\\ x_0^2 &= R \cos(\bar{\rho} \theta_\mathrm{max}) + C^2. \end{align}

Here $$\bar{\rho}^2 \equiv (\bar{x}^2+\bar{y}^2)/\bar{R}^2$$, where $$\bar{R}$$ is the radius of the cylinder in barred coordinates, which is always unity, and where $$\theta_\mathrm{max}$$ is defined by $$\cos(\theta_\mathrm{max}) = (z_\mathrm{P}-C^2)/R$$. Note that when $$\bar{\rho}=0$$, we must evaluate $$\sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}$$ as $$\theta_\mathrm{max}$$.

sigma

$$\sigma$$ is a function that is zero at $$\bar{z}=-1$$ (which maps onto the sphere $$x^i=x_0^i$$) and unity at $$\bar{z}=+1$$ (corresponding to the upper surface of the FocallyLiftedMap). We define

\begin{align} \sigma &= \frac{\bar{z}+1}{2}. \end{align}

deriv_sigma

deriv_sigma returns

\begin{align} \frac{\partial \sigma}{\partial \bar{x}^j} &= (0,0,1/2). \end{align}

jacobian

jacobian returns $$\partial x_0^k/\partial \bar{x}^j$$. The arguments to jacobian are $$(\bar{x},\bar{y},\bar{z})$$, but $$\bar{z}$$ is ignored.

Differentiating Eqs.(1–3) above yields

\begin{align*} \frac{\partial x_0^2}{\partial \bar{x}} &= - R \theta_\mathrm{max} \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\bar{x}\\ \frac{\partial x_0^2}{\partial \bar{y}} &= - R \theta_\mathrm{max} \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\bar{y}\\ \frac{\partial x_0^0}{\partial \bar{x}} &= R \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}} + R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{x}^2\\ \frac{\partial x_0^0}{\partial \bar{y}} &= R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{x}\bar{y}\\ \frac{\partial x_0^1}{\partial \bar{x}} &= R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{x}\bar{y}\\ \frac{\partial x_0^1}{\partial \bar{y}} &= R \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}} + R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{y}^2\\ \frac{\partial x_0^i}{\partial \bar{z}} &=0. \end{align*}

Evaluating sinc functions

Note that $$\sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}$$ and its derivative appear in the above equations. We evaluate $$\sin(ax)/x$$ in a straightforward way, except we are careful to evaluate $$\sin(ax)/x = a$$ for the special case $$x=0$$.

The derivative of the sync function is more complicated to evaluate because of roundoff. Note that we can expand

\begin{align*} \frac{1}{x}\frac{d}{dx}\left(\frac{\sin(ax)}{x}\right) &= \frac{a}{x^2}\left(1 - 1 - \frac{2 (ax)^2}{3!} + \frac{4(ax)^4}{5!} - \frac{5(ax)^6}{7!} + \cdots \right), \end{align*}

where we kept the "1 - 1" above as a reminder that when evaluating this function directly as $$(a \cos(ax)/x^2 - \sin(ax)/x^3)$$, there can be significant roundoff because of the "1" in each of the two terms that are subtracted. The relative error in direct evaluation is $$3\epsilon/(ax)^2$$, where $$\epsilon$$ is machine precision (This expression comes from replacing "1 - 1" with $$\epsilon$$ and comparing the lowest-order contribution to the correct answer, i.e. the $$2(ax)^2/3!$$ term, with $$\epsilon$$, the error contribution). This means the error is 100% if $$(ax)^2=3\epsilon$$.

To avoid roundoff, we evaluate the series if $$ax$$ is small enough. Suppose we keep terms up to and including the $$(ax)^{2n}$$ term in the series. Then we evaluate the series if the next term, the $$(ax)^{2n+2}$$ term, is roundoff, i.e. if $$(2n+2)(ax)^{2n+2}/(2n+3)! < \epsilon$$. In this case, the direct evaluation has the maximum error if $$(2n+2)(ax)^{2n+2}/(2n+3)! = \epsilon$$. We showed above that the relative error in direct evaluation is $$3\epsilon/(ax)^2$$, which evaluates to $$(\epsilon^n (2n+2)/(2n+3)!)^{1/(n+1)}$$.

We gain less and less with each order, so we choose $$n=3$$. Then the series above can be rewritten to this order in the form

\begin{align*} \frac{1}{x}\frac{d}{dx}\left(\frac{\sin(ax)}{x}\right) &= -\frac{a^3}{3}\left(1 - \frac{3\cdot 4 (ax)^2}{5!} + \frac{3 \cdot 6(ax)^4}{7!}\right). \end{align*}

inverse

inverse takes $$x_0^i$$ and $$\sigma$$ as arguments, and returns $$(\bar{x},\bar{y},\bar{z})$$, or boost::none if $$x_0^i$$ or $$\sigma$$ are outside the range of the map. For example, if $$x_0^i$$ does not lie on the sphere, we return boost::none.

The easiest to compute is $$\bar{z}$$, which is given by inverting Eq. (4):

\begin{align} \bar{z} &= 2\sigma - 1. \end{align}

If $$\bar{z}$$ is outside the range $$[-1,1]$$ then we return boost::none.

To get $$\bar{x}$$ and $$\bar{y}$$, we invert Eqs (1–3). If $$x_0^0=x_0^1=0$$, then $$\bar{x}=\bar{y}=0$$. Otherwise, we compute

\begin{align} \bar{\rho} = \theta_\mathrm{max}^{-1} \tan^{-1}\left(\frac{\rho}{x_0^2-C^2}\right), \end{align}

where $$\rho^2 = (x_0^0-C^0)^2+(x_0^1-C^1)^2$$. Then

\begin{align} \bar{x} &= (x_0^0-C^0)\frac{\bar{\rho}}{\rho},\\ \bar{y} &= (x_0^1-C^1)\frac{\bar{\rho}}{\rho}. \end{align}

Note that if $$\bar{x}^2+\bar{y}^2 > 1$$, the original point is outside the range of the map so we return boost::none.

lambda_tilde

lambda_tilde takes as arguments a point $$x^i$$ and a projection point $$P^i$$, and computes $$\tilde{\lambda}$$, the solution to

\begin{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\end{align}

Since $$x_0^i$$ must lie on the sphere, $$\tilde{\lambda}$$ is the solution of the quadratic equation

\begin{align} |P^i + (x^i - P^i) \tilde{\lambda} - C^i |^2 - R^2 = 0. \end{align}

In solving the quadratic, we demand a root that is positive and less than or equal to unity, since $$x_0^i$$ is always between the projection point and $$x^i$$. If there are two suitable roots, this means that the entire sphere lies between $$P^i$$ and $$x^i$$; in this case if $$x^2 \geq z_\mathrm{P}$$ we choose the larger root, otherwise we choose the smaller one: this gives us the root with $$x_0^2 \geq z_\mathrm{P}$$, the portion of the sphere that is the range of Endcap. If there is no suitable root, this means that the point $$x^i$$ is not in the range of the map so we return a default-constructed std::optional.

deriv_lambda_tilde

deriv_lambda_tilde takes as arguments $$x_0^i$$, a projection point $$P^i$$, and $$\tilde{\lambda}$$, and returns $$\partial \tilde{\lambda}/\partial x^i$$. By differentiating Eq. (11), we find

\begin{align} \frac{\partial\tilde{\lambda}}{\partial x^j} &= \tilde{\lambda}^2 \frac{C^j - x_0^j}{|x_0^i - P^i|^2 + (x_0^i - P^i)(P_i - C_{i})}. \end{align}

inv_jacobian

inv_jacobian returns $$\partial \bar{x}^i/\partial x_0^k$$, where $$\sigma$$ is held fixed. The arguments to inv_jacobian are $$(\bar{x},\bar{y},\bar{z})$$, but $$\bar{z}$$ is ignored.

Note that $$\bar{x}$$ and $$\bar{y}$$ can be considered to depend only on $$x_0^0$$ and $$x_0^1$$ but not on $$x_0^2$$, because the point $$x_0^i$$ is constrained to lie on a sphere of radius $$R$$. Note that there is an alternative way to compute Eqs. (8) and (9) using only $$x_0^0$$ and $$x_0^1$$. To do this, define

\begin{align} \upsilon \equiv \sin(\bar{\rho}\theta_\mathrm{max}) &= \sqrt{\frac{(x_0^0-C^0)^2+(x_0^1-C^1)^2}{R^2}}. \end{align}

Then we can write

\begin{align} \frac{1}{\bar{\rho}}\sin(\bar{\rho}\theta_\mathrm{max}) &= \frac{\theta_\mathrm{max}\upsilon}{\arcsin(\upsilon)}, \end{align}

so that

\begin{align} \bar{x} &= \frac{x_0^0-C^0}{R}\left(\frac{1}{\bar{\rho}} \sin(\bar{\rho}\theta_\mathrm{max})\right)^{-1} \\ \bar{y} &= \frac{x_0^1-C^1}{R}\left(\frac{1}{\bar{\rho}} \sin(\bar{\rho}\theta_\mathrm{max})\right)^{-1}. \end{align}

We will compute $$\partial \bar{x}^i/\partial x_0^k$$ by differentiating Eqs. (15) and (16). Because those equations involve $$\bar{\rho}$$, we first establish some relations involving derivatives of $$\bar{\rho}$$. For ease of notation, we define

\begin{align} q \equiv \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}. \end{align}

First observe that

\begin{align} \frac{dq}{d\upsilon} = \frac{dq}{d\bar{\rho}} \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1}, \end{align}

where $$\upsilon$$ is the quantity defined by Eq. (13). Therefore

\begin{align} \frac{\partial q}{\partial x_0^0} &= \frac{\bar{x}}{\bar{\rho}R}\frac{dq}{d\bar{\rho}} \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1},\\ \frac{\partial q}{\partial x_0^1} &= \frac{\bar{y}}{\bar{\rho}R}\frac{dq}{d\bar{\rho}} \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1}, \end{align}

where we have differentiated Eq. (13), and where we have used Eqs. (15) and (16) to eliminate $$x_0^0$$ and $$x_0^1$$ in favor of $$\bar{x}$$ and $$\bar{y}$$ in the final result.

Let

\begin{align} \Sigma \equiv \frac{1}{\bar\rho} \frac{dq}{d\bar{\rho}}, \end{align}

since that combination will appear frequently in formulas below. Note that $$\Sigma$$ has a finite limit as $$\bar{\rho}\to 0$$, and it is evaluated according to the section on evaluating sinc functions above.

By differentiating Eqs. (15) and (16), and using Eqs. (19) and (20), we find

\begin{align} \frac{\partial \bar{x}}{\partial x_0^0} &= \frac{1}{R q} - \frac{\bar{x}^2 \Sigma}{R q} \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\ \frac{\partial \bar{x}}{\partial x_0^1} &= - \frac{\bar{x}\bar{y} \Sigma}{R q} \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\ \frac{\partial \bar{x}}{\partial x_0^2} &= 0,\\ \frac{\partial \bar{y}}{\partial x_0^0} &= \frac{\partial \bar{x}}{\partial x_0^1},\\ \frac{\partial \bar{y}}{\partial x_0^1} &= \frac{1}{R q} - \frac{\bar{y}^2 \Sigma}{R q} \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\ \frac{\partial \bar{y}}{\partial x_0^2} &= 0,\\ \frac{\partial \bar{z}}{\partial x_0^i} &= 0. \end{align}

Note that care must be taken to evaluate $$q = \sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}$$ and its derivative $$\Sigma$$ near $$\bar{\rho}=0$$; see the discussion above on evaluating sinc functions.

dxbar_dsigma

dxbar_dsigma returns $$\partial \bar{x}^i/\partial \sigma$$, where $$x_0^i$$ is held fixed.

From Eq. (6) we have

\begin{align} \frac{\partial \bar{x}^i}{\partial \sigma} &= (0,0,2). \end{align}

The documentation for this class was generated from the following file:
• src/Domain/CoordinateMaps/FocallyLiftedEndcap.hpp