|
SpECTRE
v2024.04.12
|
Map from 3D unit right cylinder to a 3D volume that connects a portion of a spherical surface with a disk. More...
#include <UniformCylindricalFlatEndcap.hpp>
Public Member Functions | |
| UniformCylindricalFlatEndcap (const std::array< double, 3 > ¢er_one, const std::array< double, 3 > ¢er_two, double radius_one, double radius_two, double z_plane_one) | |
| UniformCylindricalFlatEndcap (UniformCylindricalFlatEndcap &&)=default | |
| UniformCylindricalFlatEndcap (const UniformCylindricalFlatEndcap &)=default | |
| UniformCylindricalFlatEndcap & | operator= (const UniformCylindricalFlatEndcap &)=default |
| UniformCylindricalFlatEndcap & | operator= (UniformCylindricalFlatEndcap &&)=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 UniformCylindricalFlatEndcap &lhs, const UniformCylindricalFlatEndcap &rhs) |
Map from 3D unit right cylinder to a 3D volume that connects a portion of a spherical surface with a disk.
Consider a sphere with center \(C_1\) and radius \(R_1\), and a disk in the \(xy\) plane with center \(C_2\) and radius \(R_2\). Let sphere 1 be intersected by a plane normal to the \(z\) axis and located at \(z = z_{\mathrm{P}1}\),
UniformCylindricalFlatEndcap maps 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\)) to the shaded area in the figure above (with coordinates \((x,y,z)\)). The "bottom" of the cylinder \(\bar{z}=-1\) is mapped to the portion of sphere 1 that has \(z \geq z_{\mathrm{P}1}\), and on this portion of the sphere the angular coordinate \(\theta_1 = \acos((z-C_1^z)/R_1)\) is uniform in \(\bar{\rho} = \sqrt{\bar{x}^2+\bar{y}^2}\) and the angular coordinate \(\phi_1 = \atan((y-C_1^y)/(x-C_1^x))\) is the same as \(\phi = \atan(\bar{y}/\bar{x})\). The "top" of the cylinder \(\bar{z}=+1\) is mapped to the disk, and the radial polar coordinate coordinate \(\sqrt{x^2+y^2}\) on this disk is equal to \(R_2\bar\rho\) and the angular coordinate \(\phi_2 = \atan((y-C_2^y)/(x-C_2^x))\) is the same as \(\phi\).
UniformCylindricalFlatEndcap is intended to be composed with Wedge<2> maps to construct a portion of a cylindrical domain for a binary system.
UniformCylindricalFlatEndcap can be used to construct a domain that is similar to, but not identical to, the one described briefly in the Appendix of [32]. UniformCylindricalFlatEndcap is used to construct the Blocks analogous to those labeled 'MA wedge' and 'MB wedge' in Figure 20 of that paper.
UniformCylindricalFlatEndcap provides the following functions:
operator() maps \((\bar{x},\bar{y},\bar{z})\) to \((x,y,z)\) according to
\begin{align} x &= C_1^x+\lambda(C_2^x-C_1^x) + \cos\phi\left(R_1\sin\theta_1 + \lambda(R_2\bar\rho-R_1\sin\theta_1)\right), \label{eq:x0} \\ y &= C_1^y+\lambda(C_2^y-C_1^y) + \sin\phi\left(R_1\sin\theta_1 + \lambda(R_2\bar\rho-R_1\sin\theta_1)\right), \label{eq:x1} \\ z &= C_1^z+\lambda(C_2^z-C_1^z) + (1-\lambda)R_1\cos\theta_1 \label{eq:x2}. \end{align}
Here
\begin{align} \lambda &= \frac{\bar{z}+1}{2},\label{eq:lambdafromzbar}\\ \theta_1 &= \bar{\rho} \theta_{1 \mathrm{max}},\label{eq:deftheta1}\\ \phi &= \atan(\bar{y}/\bar{x})\label{eq:defphi}, \end{align}
where \(\theta_{1 \mathrm{max}}\) is defined by
\begin{align} \cos(\theta_{1\mathrm{max}}) &= (z_{\mathrm{P}1}-C_1^z)/R_1,\\ \end{align}
and
\begin{align} \bar{\rho} &= \sqrt{\bar{x}^2+\bar{y}^2}/\bar{R} \label{eq:defrhobar}, \end{align}
where \(\bar{R}\) is the radius of the cylinder in barred coordinates, which is always unity.
Given \((x,y,z)\) we want to find \((\bar{x},\bar{y},\bar{z})\). From Eq. ( \(\ref{eq:x2}\)) we can write \(\lambda\) as a function of \(\bar\rho\):
\begin{align} \lambda &= \frac{z - C_1^z - R_1\cos\theta_1} {C_2^z-C_1^z - R_1\cos\theta_1} \label{eq:lambda_from_rho}. \end{align}
Then by eliminating \(\phi\) from Eqs. ( \(\ref{eq:x0}\)) and ( \(\ref{eq:x1}\)) we find that \(\bar{\rho}\) is the solution of \(Q(\bar{\rho})=0\), where
\begin{align} Q(\bar{\rho}) &= \left(x-C_1^x-\lambda(C_2^x-C_1^x)\right)^2+ \left(y-C_1^y-\lambda(C_2^y-C_1^y)\right)^2- \left((1-\lambda)R_1\sin\theta_1 + \lambda \bar\rho R_2\right)^2.\label{eq:defQ} \end{align}
Here \(\lambda\) and \(\theta_1\), are functions of \(\bar{\rho}\) through Eqs. ( \(\ref{eq:lambda_from_rho}\)) and ( \(\ref{eq:deftheta1}\)).
We solve \(Q(\bar{\rho})=0\) numerically; it is a one-dimensional root-finding problem.
Once we have determined \(\bar{\rho}\), we then obtain \(\lambda\) from Eq. ( \(\ref{eq:lambda_from_rho}\)), and we obtain \(\phi\) from
\begin{align} \tan\phi &= \frac{y-C_1^y-\lambda(C_2^y-C_1^y)}{x-C_1^x-\lambda(C_2^x-C_1^x)}. \end{align}
Then \(\bar{z}\) is obtained from Eq. ( \(\ref{eq:lambdafromzbar}\)) and \(\bar{x}\) and \(\bar{y}\) are obtained from
\begin{align} \bar{x} &= \bar{\rho}\bar{R}\cos\phi,\\ \bar{y} &= \bar{\rho}\bar{R}\sin\phi. \end{align}
We solve \(Q(\bar{\rho})=0\) numerically for \(\bar{\rho}\), where \(Q(\bar{\rho})\) is given by Eq. ( \(\ref{eq:defQ}\)).
Note that the root we care about must have \(0\leq\lambda\leq 1\); therefore from Eq. ( \(\ref{eq:lambda_from_rho}\)) we have
\begin{align} \bar{\rho}_{\mathrm{min}} &= \left\{\begin{array}{ll} 0 & \text{for } z-C_1^z \geq R_1, \\ \displaystyle \frac{1}{\theta_{1 \mathrm{max}}} \cos^{-1}\left(\frac{z-C_1^z}{R_1}\right) & \text{otherwise} \end{array}\right.\label{eq:rhobarmin}\\ \bar{\rho}_{\mathrm{max}} &= 1. \label{eq:rhobarmax} \end{align}
so we look for a root only between \(\bar{\rho}_{\mathrm{min}}\) and \(\bar{\rho}_{\mathrm{max}}\).
Sometimes a root is within roundoff of \(\bar{\rho}_{\mathrm{min}}\) This tends to happen at points on the boundary of the mapped region. In this case, the root might not be bracketed by \([\bar{\rho}_{\mathrm{min}},\bar{\rho}_{\mathrm{max}}]\) if the root is slightly outside that interval. If we find that \(Q(\bar{\rho}_{\mathrm{min}})\) is near zero but has the wrong sign, then we slightly expand the interval as follows:
\begin{align} \bar{\rho}_{\mathrm{min}} \to \bar{\rho}_{\mathrm{min}} - 2 \frac{Q(\bar{\rho}_{\mathrm{min}})}{Q'(\bar{\rho}_{\mathrm{min}})}, \end{align}
where \(Q'(\bar{\rho}_{\mathrm{min}})\) is the derivative of the function in Eq. ( \(\ref{eq:defQ}\)). Note that without the factor of 2, this is a Newton-Raphson step; the factor of 2 is there to overcompensate so that the new \(\bar{\rho}_{\mathrm{min}}\) brackets the root. Sometimes, if the derivative is large enough so that the correction above amounts to a value less than roundoff; in that case, we increase the correction to a value larger than roundoff.
Note that by differentiating Eqs. ( \(\ref{eq:defQ}\)) and ( \(\ref{eq:lambda_from_rho}\)), one obtains
\begin{align} Q'(\bar{\rho}) =& -2 \frac{d\lambda}{d\bar{\rho}}\left[ \left(x-C_1^x-\lambda(C_2^x-C_1^x)\right)(C_2^x-C_1^x)+ \left(y-C_1^y-\lambda(C_2^y-C_1^y)\right)(C_2^y-C_1^y) \right]\nonumber \\ & -2 \left((1-\lambda)R_1\sin\theta_1 + \lambda \bar\rho R_2\right) \left[ \frac{d\lambda}{d\bar{\rho}} (R_2\bar\rho - R_1\sin\theta_1) +(1-\lambda)R_1\theta_{1 \mathrm{max}}\cos\theta_1 +\lambda R_2 \right], \label{eq:defQderiv} \end{align}
where
\begin{align} \frac{d\lambda}{d\bar{\rho}} &= \frac{(1-\lambda)R_1\theta_{1 \mathrm{max}}\sin\theta_1} {C_2^z-C_1^z -R_1\cos\theta_1} \label{eq:dlambda_drhobar}. \end{align}
For some points on the boundary of the mapped domain, the root will be within roundoff of \(\bar{\rho}=0\) or \(\bar{\rho}=1\). Here it does not always make sense to expand the range of the map if the root fails (by roundoff) to be bracketed, as is done above. Furthermore, when \(\bar{\rho}=0\) is a root it turns out that both \(Q(\bar{\rho})=0\) and \(Q'(\bar{\rho})=0\) for \(\bar{\rho}=0\), so some root-finders (e.g. Newton-Raphson) have difficulty converging. Therefore the cases where the root is within roundoff of \(\bar{\rho}=0\) or \(\bar{\rho}=1\) are treated separately.
These cases are detected by comparing terms in the first-order power series of \(Q(\bar{\rho})=0\) when expanded about \(\bar{\rho}=0\) or \(\bar{\rho}=1\). When one of these cases is recognized, the root is returned as either \(\bar{\rho}=0\) or \(\bar{\rho}=1\).
It is expected that inverse() will often be passed points \((x,y,z)\) that are out of the range of the map; in this case inverse() returns a std::nullopt. To avoid the difficulty and expense of attempting to solve \(Q(\bar{\rho})=0\) numerically for such points (and then having this solution fail), it is useful to quickly reject points \((x,y,z)\) that are outside the range of the map.
Any point in the range of the map must be below the disk and it must be outside or on sphere 1, so the inverse map can immediately return a std::nullopt for a point that does not satisfy these conditions.
Likewise, the inverse map can immediately reject any point with \(z < z_{\mathrm{P}1}\).
Finally, consider the circle \(S_1\) defining the intersection of sphere 1 and the plane \(z = z_{\mathrm{P}1}\); this circle has radius \(r_1 = R_1 \sin\theta_{1\mathrm{max}}\). Now consider the cone that passes through both \(S_1\) and the circle \(S_2\) bounding the upper disk. A point in the range of the map must be inside or on this cone. The cone can be defined parametrically as
\begin{align} x_c &= C_1^x + \tilde{\lambda}(C_2^x-C_1^x) + \cos\varphi (r_1 + \tilde{\lambda} (R_2 -r_1)),\\ y_c &= C_1^y + \tilde{\lambda}(C_2^y-C_1^y),+ \sin\varphi (r_1 + \tilde{\lambda} (R_2 -r_1)),\\ z_c &= C_1^z + R_1 \cos\theta_{1\mathrm{max}} + \tilde{\lambda}(C_2^z - C_1^z - R_1 \cos\theta_{1\mathrm{max}}), \end{align}
where \((x_c,y_c,z_c)\) is a point on the cone, and the two parameters defining a point on the cone are the angle \(\varphi\) around the cone and the parameter \(\tilde{\lambda}\), which is defined to be zero on \(S_1\) and unity on \(S_2\).
Given an arbitrary point \((x, y, z)\), we can determine whether or not that point is inside the cone as follows. First determine
\begin{align} \tilde{\lambda} &= \frac{z - C_1^z - R_1 \cos\theta_{1\mathrm{max}}} {C_2^z - C_1^z- R_1 \cos\theta_{1\mathrm{max}}}, \\ \tilde{x} &= x - C_1^x - \tilde{\lambda} (C_2^x-C_1^x),\\ \tilde{y} &= y - C_1^y - \tilde{\lambda} (C_2^y-C_1^y).\\ \end{align}
Then the condition for the point to be inside or on the cone is
\begin{align} \sqrt{\tilde{x}^2+\tilde{y}^2} \le r_1 + (R_2-r_1)\tilde{\lambda}. \end{align}
The inverse map can therefore reject any points that do not satisfy this criterion.
One can rewrite Eqs.( \(\ref{eq:x0}\)) through ( \(\ref{eq:x2}\)) as
\begin{align} x &= \frac{1}{2}\left((1-\bar{z})C_1^x+ (1+\bar{z})C_2^x\right) + \frac{\bar{x}}{2}\left( (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 \right), \label{eq:x0alt} \\ y &= \frac{1}{2}\left((1-\bar{z})C_1^y + (1+\bar{z})C_2^y\right) + \frac{\bar{y}}{2}\left( (1-\bar{z})R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z})R_2 \right), \label{eq:x1alt} \\ z &= \frac{1}{2}\left((1-\bar{z})C_1^z + (1+\bar{z})C_2^z\right) + \frac{1}{2} (1-\bar{z})R_1 \cos\theta_1, \label{eq:x2alt} \\ \end{align}
where we have used Eq. ( \(\ref{eq:lambdafromzbar}\)) to eliminate \(\lambda\) in favor of \(\bar{z}\), and where we have defined the function
\begin{align} S(\bar{\rho},a) = \frac{\sin(\bar{\rho} a)}{\bar{\rho}}. \label{eq:Sdef} \end{align}
Note that \(S(\bar{\rho},a)\) is finite as \(\bar{\rho}\) approaches zero, and in the code we must take care that everything remains well-behaved in that limit.
Then differentiating Eqs. ( \(\ref{eq:x0alt}\)) and ( \(\ref{eq:x1alt}\)) with respect to \(\bar{x}\) and \(\bar{y}\), taking into account the dependence of \(\bar{\rho}\) on \(\bar{x}\) and \(\bar{y}\) from Eq. ( \(\ref{eq:defrhobar}\)), we find:
\begin{align} \frac{\partial x^0}{\partial \bar{x}} &= \frac{1}{2}\left( (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 \right) + \frac{\bar{x}^2}{2\bar{\rho}} (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}), \\ \frac{\partial x^1}{\partial \bar{y}} &= \frac{1}{2}\left( (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 \right) + \frac{\bar{y}^2}{2\bar{\rho}} (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}), \\ \frac{\partial x^0}{\partial \bar{y}} &= \frac{\bar{x}\bar{y}}{2\bar{\rho}} (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}), \\ \frac{\partial x^1}{\partial \bar{x}} &= \frac{\partial x^0}{\partial \bar{y}}, \end{align}
where \(S'(\bar{\rho},a)\) means the derivative of \(S(\bar{\rho},a)\) with respect to \(\bar\rho\). Note that \(S'(\bar{\rho},a)/\bar{\rho}\) approaches a constant value as \(\bar{\rho}\) approaches zero.
Differentiating Eq. ( \(\ref{eq:x2alt}\)) with respect to \(\bar{x}\) and \(\bar{y}\) we find
\begin{align} \frac{\partial z}{\partial \bar{x}} &= - \frac{\bar{x}}{2} (1-\bar{z}) R_1 \theta_{1 \mathrm{max}} S(\bar{\rho},\theta_{1 \mathrm{max}}),\\ \frac{\partial z}{\partial \bar{y}} &= - \frac{\bar{y}}{2} (1-\bar{z}) R_1 \theta_{1 \mathrm{max}} S(\bar{\rho},\theta_{1 \mathrm{max}}). \end{align}
Differentiating Eqs. ( \(\ref{eq:x0alt}\)) through ( \(\ref{eq:x2alt}\)) with respect to \(\bar{z}\) yields
\begin{align} \frac{\partial x}{\partial \bar{z}} &= \frac{1}{2}\left[ C_2^x-C_1^x + \bar{x}\left(R_2 - R_1 S(\bar{\rho},\theta_{1 \mathrm{max}})\right) \right],\\ \frac{\partial y}{\partial \bar{z}} &= \frac{1}{2}\left[ C_2^y-C_1^y + \bar{y}\left(R_2 - R_1 S(\bar{\rho},\theta_{1 \mathrm{max}})\right) \right],\\ \frac{\partial z}{\partial \bar{z}} &= \frac{1}{2}\left( C_2^z-C_1^z - R_1\cos\theta_1 \right). \end{align}
The inverse Jacobian is computed by numerically inverting the Jacobian.
We demand that \(C^2_1 + 1.05 R_1 \leq C^2_2 \leq C^2_1 + 5 R_1\). It is possible to construct a valid map without this assumption, but the assumption simplifies the code, and the expected use cases obey this restriction.
We also demand that the z plane in the above figure lies above the center of the sphere and is not too close to the center or edge of the sphere; specifically, we demand that
Here 0.075 and 0.35 are safety factors. These restrictions are not strictly necessary but are made for simplicity and to ensure the accuracy of the inverse map (the inverse map becomes less accurate if the map parameters are extreme).
Consider the line segment \(L\) that connects a point on the circle \(S_1\) (the circle formed by the intersection of sphere 1 and the plane \(z=z_{\mathrm{P}1}\)) with the center of the circle \(S_1\). Consider another line segment \(L'\) that connects the same point on the circle \(S_1\) with the corresponding point on the circle \(S_2\) (the circle bounding the disk with center \(C_2\) and radius \(R_2\)). Now consider the angle between \(L\) and \(L'\), as measured from the interior of sphere 1, and Let \(\alpha\) be the minimum value of this angle over the circle. \(\alpha\) is shown in the figure above. If \(\alpha < \theta_{1 \mathrm{max}}\), then the line segment \(L'\) intersects the mapped portion of sphere 1 twice, so the map is multi-valued. Therefore we demand that the map parameters are such that
where 1.1 is a safety factor.
The condition on \(\alpha\) is guaranteed to provide an invertible map if \(C_1^x=C_2^x\) and \(C_1^y=C_2^y\). However, for \(C_1^x \neq C_2^x\) or \(C_1^y\neq C_2^y\), even if the \(\alpha\) condition is satisfied, it is possible for two lines of constant \((\bar{x},\bar{y})\) (each line has different values of \((\bar{x},\bar{y})\)) to pass through the same point \((x,y,z)\) if those lines are not coplanar. This condition is difficult to check analytically, so we check it numerically. We have found empirically that if \(Q(\bar{\rho})\) from Eq. ( \(\ref{eq:defQ}\)) has only a single root between \(\bar{\rho}_{\mathrm{min}}\) and \(\bar{\rho}_{\mathrm{max}}\) for all points \((x,y,z)\) on the surface of sphere 1 with \(z\geq z_{\mathrm{P}1}\) and with \((x-C_1^x)/(y-C_1^y) = (C_1^x-C_2^x)/(C_1^y-C_2^y)\), then the map is single-valued everywhere. We cannot numerically check every point in this one-parameter family of points, but we demand that this condition is satisfied for a reasonably large number of points (currently 1000) in this family. This check is not very expensive since it is done only once, in the constructor.