SpECTRE  v2024.09.29
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)}\).

\begin{align*} n=1 \qquad& \mathrm{error}=3(\epsilon/30)^{1/2} &\qquad \sim \mathrm{5e-09}\\ n=2 \qquad& \mathrm{error}=3(\epsilon^2/840)^{1/3} &\qquad \sim \mathrm{7e-12}\\ n=3 \qquad& \mathrm{error}=3(\epsilon^3/45360)^{1/4} &\qquad \sim \mathrm{2e-13}\\ n=4 \qquad& \mathrm{error}=3(\epsilon^4/3991680)^{1/5} &\qquad \sim \mathrm{2e-14}\\ n=5 \qquad& \mathrm{error}=3(\epsilon^5/518918400)^{1/6} &\qquad \sim \mathrm{5e-15}\\ n=6 \qquad& \mathrm{error}=3(\epsilon^6/93405312000)^{1/7} &\qquad \sim \mathrm{1e-15} \end{align*}

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: