SpECTRE  v2024.12.16
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
domain::CoordinateMaps::UniformCylindricalFlatEndcap Class Reference

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 > &center_one, const std::array< double, 3 > &center_two, double radius_one, double radius_two, double z_plane_one)
 
 UniformCylindricalFlatEndcap (UniformCylindricalFlatEndcap &&)=default
 
 UniformCylindricalFlatEndcap (const UniformCylindricalFlatEndcap &)=default
 
UniformCylindricalFlatEndcapoperator= (const UniformCylindricalFlatEndcap &)=default
 
UniformCylindricalFlatEndcapoperator= (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::NoFramejacobian (const std::array< T, 3 > &source_coords) const
 
template<typename T >
tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrameinv_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)
 

Detailed Description

Map from 3D unit right cylinder to a 3D volume that connects a portion of a spherical surface with a disk.

A cylinder maps to the shaded region.

Details

Consider a sphere with center C1 and radius R1, and a disk in the xy plane with center C2 and radius R2. Let sphere 1 be intersected by a plane normal to the z axis and located at z=zP1,

UniformCylindricalFlatEndcap maps a 3D unit right cylinder (with coordinates (x¯,y¯,z¯) such that 1z¯1 and x¯2+y¯21) to the shaded area in the figure above (with coordinates (x,y,z)). The "bottom" of the cylinder z¯=1 is mapped to the portion of sphere 1 that has zzP1, and on this portion of the sphere the angular coordinate θ1=acos((zC1z)/R1) is uniform in ρ¯=x¯2+y¯2 and the angular coordinate ϕ1=atan((yC1y)/(xC1x)) is the same as ϕ=atan(y¯/x¯). The "top" of the cylinder z¯=+1 is mapped to the disk, and the radial polar coordinate coordinate x2+y2 on this disk is equal to R2ρ¯ and the angular coordinate ϕ2=atan((yC2y)/(xC2x)) is the same as ϕ.

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 . 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()

operator() maps (x¯,y¯,z¯) to (x,y,z) according to

(1)x=C1x+λ(C2xC1x)+cosϕ(R1sinθ1+λ(R2ρ¯R1sinθ1)),(2)y=C1y+λ(C2yC1y)+sinϕ(R1sinθ1+λ(R2ρ¯R1sinθ1)),(3)z=C1z+λ(C2zC1z)+(1λ)R1cosθ1.

Here

(4)λ=z¯+12,(5)θ1=ρ¯θ1max,(6)ϕ=atan(y¯/x¯),

where θ1max is defined by

(7)cos(θ1max)=(zP1C1z)/R1,

and

(8)ρ¯=x¯2+y¯2/R¯,

where R¯ is the radius of the cylinder in barred coordinates, which is always unity.

inverse

Given (x,y,z) we want to find (x¯,y¯,z¯). From Eq. ( 3) we can write λ as a function of ρ¯:

(9)λ=zC1zR1cosθ1C2zC1zR1cosθ1.

Then by eliminating ϕ from Eqs. ( 1) and ( 2) we find that ρ¯ is the solution of Q(ρ¯)=0, where

(10)Q(ρ¯)=(xC1xλ(C2xC1x))2+(yC1yλ(C2yC1y))2((1λ)R1sinθ1+λρ¯R2)2.

Here λ and θ1, are functions of ρ¯ through Eqs. ( 9) and ( 5).

We solve Q(ρ¯)=0 numerically; it is a one-dimensional root-finding problem.

Once we have determined ρ¯, we then obtain λ from Eq. ( 9), and we obtain ϕ from

(11)tanϕ=yC1yλ(C2yC1y)xC1xλ(C2xC1x).

Then z¯ is obtained from Eq. ( 4) and x¯ and y¯ are obtained from

(12)x¯=ρ¯R¯cosϕ,(13)y¯=ρ¯R¯sinϕ.

Considerations when root-finding.

We solve Q(ρ¯)=0 numerically for ρ¯, where Q(ρ¯) is given by Eq. ( 10).

min/max values of \f$\bar{\rho}\f$:

Note that the root we care about must have 0λ1; therefore from Eq. ( 9) we have

(14)ρ¯min={0for zC1zR1,1θ1maxcos1(zC1zR1)otherwise(15)ρ¯max=1.

so we look for a root only between ρ¯min and ρ¯max.

Roots within roundoff of endpoints:

Sometimes a root is within roundoff of ρ¯min This tends to happen at points on the boundary of the mapped region. In this case, the root might not be bracketed by [ρ¯min,ρ¯max] if the root is slightly outside that interval. If we find that Q(ρ¯min) is near zero but has the wrong sign, then we slightly expand the interval as follows:

(16)ρ¯minρ¯min2Q(ρ¯min)Q(ρ¯min),

where Q(ρ¯min) is the derivative of the function in Eq. ( 10). 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 ρ¯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. ( 10) and ( 9), one obtains

Q(ρ¯)=2dλdρ¯[(xC1xλ(C2xC1x))(C2xC1x)+(yC1yλ(C2yC1y))(C2yC1y)](17)2((1λ)R1sinθ1+λρ¯R2)[dλdρ¯(R2ρ¯R1sinθ1)+(1λ)R1θ1maxcosθ1+λR2],

where

(18)dλdρ¯=(1λ)R1θ1maxsinθ1C2zC1zR1cosθ1.

Roots within roundoff of \f$\bar{\rho}=0\f$ or \f$\bar{\rho}=1\f$:

For some points on the boundary of the mapped domain, the root will be within roundoff of ρ¯=0 or ρ¯=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 ρ¯=0 is a root it turns out that both Q(ρ¯)=0 and Q(ρ¯)=0 for ρ¯=0, so some root-finders (e.g. Newton-Raphson) have difficulty converging. Therefore the cases where the root is within roundoff of ρ¯=0 or ρ¯=1 are treated separately.

These cases are detected by comparing terms in the first-order power series of Q(ρ¯)=0 when expanded about ρ¯=0 or ρ¯=1. When one of these cases is recognized, the root is returned as either ρ¯=0 or ρ¯=1.

Quick rejection of points out of range of the map.

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(ρ¯)=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<zP1.

Finally, consider the circle S1 defining the intersection of sphere 1 and the plane z=zP1; this circle has radius r1=R1sinθ1max. Now consider the cone that passes through both S1 and the circle S2 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

(19)xc=C1x+λ~(C2xC1x)+cosφ(r1+λ~(R2r1)),(20)yc=C1y+λ~(C2yC1y),+sinφ(r1+λ~(R2r1)),(21)zc=C1z+R1cosθ1max+λ~(C2zC1zR1cosθ1max),

where (xc,yc,zc) is a point on the cone, and the two parameters defining a point on the cone are the angle φ around the cone and the parameter λ~, which is defined to be zero on S1 and unity on S2.

Given an arbitrary point (x,y,z), we can determine whether or not that point is inside the cone as follows. First determine

(22)λ~=zC1zR1cosθ1maxC2zC1zR1cosθ1max,(23)x~=xC1xλ~(C2xC1x),(24)y~=yC1yλ~(C2yC1y).

Then the condition for the point to be inside or on the cone is

(25)x~2+y~2r1+(R2r1)λ~.

The inverse map can therefore reject any points that do not satisfy this criterion.

jacobian

One can rewrite Eqs.( 1) through ( 3) as

(26)x=12((1z¯)C1x+(1+z¯)C2x)+x¯2((1z¯)R1S(ρ¯,θ1max)+(1+z¯)R2),(27)y=12((1z¯)C1y+(1+z¯)C2y)+y¯2((1z¯)R1S(ρ¯,θ1max)+(1+z¯)R2),(28)z=12((1z¯)C1z+(1+z¯)C2z)+12(1z¯)R1cosθ1,

where we have used Eq. ( 4) to eliminate λ in favor of z¯, and where we have defined the function

(29)S(ρ¯,a)=sin(ρ¯a)ρ¯.

Note that S(ρ¯,a) is finite as ρ¯ approaches zero, and in the code we must take care that everything remains well-behaved in that limit.

Then differentiating Eqs. ( 26) and ( 27) with respect to x¯ and y¯, taking into account the dependence of ρ¯ on x¯ and y¯ from Eq. ( 8), we find:

(30)x0x¯=12((1z¯)R1S(ρ¯,θ1max)+(1+z¯)R2)+x¯22ρ¯(1z¯)R1S(ρ¯,θ1max),(31)x1y¯=12((1z¯)R1S(ρ¯,θ1max)+(1+z¯)R2)+y¯22ρ¯(1z¯)R1S(ρ¯,θ1max),(32)x0y¯=x¯y¯2ρ¯(1z¯)R1S(ρ¯,θ1max),(33)x1x¯=x0y¯,

where S(ρ¯,a) means the derivative of S(ρ¯,a) with respect to ρ¯. Note that S(ρ¯,a)/ρ¯ approaches a constant value as ρ¯ approaches zero.

Differentiating Eq. ( 28) with respect to x¯ and y¯ we find

(34)zx¯=x¯2(1z¯)R1θ1maxS(ρ¯,θ1max),(35)zy¯=y¯2(1z¯)R1θ1maxS(ρ¯,θ1max).

Differentiating Eqs. ( 26) through ( 28) with respect to z¯ yields

(36)xz¯=12[C2xC1x+x¯(R2R1S(ρ¯,θ1max))],(37)yz¯=12[C2yC1y+y¯(R2R1S(ρ¯,θ1max))],(38)zz¯=12(C2zC1zR1cosθ1).

inv_jacobian

The inverse Jacobian is computed by numerically inverting the Jacobian.

Restrictions on map parameters

We demand that C12+1.05R1C22C12+5R1. 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

  • 0.075π<θ1max<0.35π

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 S1 (the circle formed by the intersection of sphere 1 and the plane z=zP1) with the center of the circle S1. Consider another line segment L that connects the same point on the circle S1 with the corresponding point on the circle S2 (the circle bounding the disk with center C2 and radius R2). Now consider the angle between L and L, as measured from the interior of sphere 1, and Let α be the minimum value of this angle over the circle. α is shown in the figure above. If α<θ1max, 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

  • α>1.1θ1max

where 1.1 is a safety factor.

The condition on α is guaranteed to provide an invertible map if C1x=C2x and C1y=C2y. However, for C1xC2x or C1yC2y, even if the α condition is satisfied, it is possible for two lines of constant (x¯,y¯) (each line has different values of (x¯,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(ρ¯) from Eq. ( 10) has only a single root between ρ¯min and ρ¯max for all points (x,y,z) on the surface of sphere 1 with zzP1 and with (xC1x)/(yC1y)=(C1xC2x)/(C1yC2y), 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.


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