CylindricalEndcap.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines the class CylindricalEndcap.
6 
7 #pragma once
8 
9 #include <array>
10 #include <cstddef>
11 #include <limits>
12 #include <optional>
13 
15 #include "Domain/CoordinateMaps/FocallyLiftedEndcap.hpp"
16 #include "Domain/CoordinateMaps/FocallyLiftedMap.hpp"
17 #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
18 
19 /// \cond
20 namespace PUP {
21 class er;
22 } // namespace PUP
23 /// \endcond
24 
25 namespace domain::CoordinateMaps {
26 
27 /*!
28  * \ingroup CoordinateMapsGroup
29  *
30  * \brief Map from 3D unit right cylinder to a volume that connects
31  * portions of two spherical surfaces.
32  *
33  * \image html CylindricalEndcap.svg "A cylinder maps to the shaded region."
34  *
35  * \details Consider two spheres with centers \f$C_1\f$ and \f$C_2\f$,
36  * and radii \f$R_1\f$ and \f$R_2\f$. Let sphere 1 be intersected by a
37  * plane normal to the \f$z\f$ axis and located at \f$z = z_\mathrm{P}\f$.
38  * Also let there be a projection point \f$P\f$.
39  *
40  * CylindricalEndcap maps a 3D unit right cylinder (with coordinates
41  * \f$(\bar{x},\bar{y},\bar{z})\f$ such that \f$-1\leq\bar{z}\leq 1\f$
42  * and \f$\bar{x}^2+\bar{y}^2 \leq 1\f$) to the shaded area
43  * in the figure above (with coordinates \f$(x,y,z)\f$). The "bottom"
44  * of the cylinder \f$\bar{z}=-1\f$ is mapped to the portion of sphere
45  * 1 that has \f$z \geq z_\mathrm{P}\f$. Curves of constant
46  * \f$(\bar{x},\bar{y})\f$ are mapped to portions of lines that pass
47  * through \f$P\f$. Along each of these curves, \f$\bar{z}=-1\f$ is
48  * mapped to a point on sphere 1 and \f$\bar{z}=+1\f$ is mapped to a
49  * point on sphere 2.
50  *
51  * Note that Sphere 1 and Sphere 2 are not equivalent, because the
52  * mapped portion of Sphere 1 is bounded by a plane of constant
53  * \f$z\f$ but the mapped portion of Sphere 2 is not (except for
54  * special choices of \f$C_1\f$, \f$C_2\f$, and \f$P\f$).
55  *
56  * CylindricalEndcap is intended to be composed with Wedge2D maps to
57  * construct a portion of a cylindrical domain for a binary system.
58  *
59  * CylindricalEndcap is described briefly in the Appendix of
60  * \cite Buchman:2012dw.
61  * CylindricalEndcap is used to construct the blocks labeled 'CA
62  * wedge', 'EA wedge', 'CB wedge', 'EE wedge', and 'EB wedge' in
63  * Figure 20 of that paper. Note that 'CA wedge', 'CB wedge', and
64  * 'EE wedge' have Sphere 1 contained in Sphere 2, and 'EA wedge'
65  * and 'EB wedge' have Sphere 2 contained in Sphere 1.
66  *
67  * CylindricalEndcap is implemented using `FocallyLiftedMap`
68  * and `FocallyLiftedInnerMaps::Endcap`; see those classes for details.
69  *
70  * ### Restrictions on map parameters.
71  *
72  * We demand that either Sphere 1 is fully contained inside Sphere 2, or
73  * that Sphere 2 is fully contained inside Sphere 1. It is
74  * possible to construct a valid map without this assumption, but the
75  * assumption simplifies the code, and the expected use cases obey
76  * this restriction.
77  *
78  * We also demand that \f$z_\mathrm{P} > C_1^2\f$, that is, the plane
79  * in the above and below figures lies to the right of \f$C_1^2\f$.
80  * This restriction not strictly necessary but is made for simplicity.
81  *
82  * The map is invertible only for some choices of the projection point
83  * \f$P\f$. Given the above restrictions, the allowed values of
84  * \f$P\f$ are illustrated by the following diagram:
85  *
86  * \image html CylindricalEndcap_Allowed.svg "Allowed region for P." width=75%
87  *
88  * The plane \f$z=z_\mathrm{P}\f$ intersects sphere 1 on a circle. The
89  * cone with apex \f$C_1\f$ that intersects that circle has opening
90  * angle \f$2\theta\f$ as shown in the above figure. Construct another
91  * cone, the "invertibility cone", with apex \f$S\f$ chosen such that
92  * the two cones intersect at right angles on the circle; thus the
93  * opening angle of the invertibility cone is \f$\pi-2\theta\f$. A
94  * necessary condition for invertibility is that the projection point
95  * \f$P\f$ lies inside the invertibility cone, but not between \f$S\f$
96  * and sphere 1. (If \f$P\f$ does not obey this condition, then from the
97  * diagram one can find at least one line through \f$P\f$
98  * that twice intersects the surface of sphere 1 with \f$z>z_\mathrm{P}\f$;
99  * the inverse map is thus double-valued at those intersection points.)
100  * Placing the projection point \f$P\f$ to the
101  * right of \f$S\f$ (but inside the invertibility cone) is ok for
102  * invertibility.
103  *
104  * In addition to invertibility and the two additional restrictions
105  * already mentioned above, we demand a few more restrictions on the
106  * map parameters to simplify the logic for the expected use cases and
107  * to ensure that jacobians do not get too large. The numbers in the
108  * restrictions below were chosen empirically so that the unit tests pass
109  * with errors less than 100 times machine roundoff; we do not expect to
110  * run into these restrictions in normal usage. We demand:
111  *
112  * - \f$P\f$ is not too close to the edge of the invertibility cone.
113  * Here we demand that the angle between \f$P\f$ and \f$S\f$
114  * is less than \f$0.85 (\pi/2-\theta)\f$ (note that if this
115  * angle is exactly \f$\pi/2-\theta\f$ it is exactly on the
116  * invertibility cone); The 0.85 was chosen empirically based on
117  * unit tests.
118  * - \f$z_\mathrm{P}\f$ is not too close to the center or the edge of sphere 1.
119  * Here we demand that \f$0.2 \leq \cos(\theta) \leq 0.95\f$, where
120  * the values 0.2 and 0.95 were chosen empirically based on unit tests.
121  * - \f$P\f$ is contained in sphere 2.
122  * - If sphere 2 is contained in sphere 1, then
123  * - \f$0.25 R_1 \leq R_2 \leq 0.85 (R_1 - |C_1-C_2|)\f$,
124  * This prevents the two spheres from having a very narrow space between
125  * them, and it prevents sphere 2 from being very small.
126  * - \f$|P - C_2| < 0.1 R_2\f$, i.e. \f$P\f$ is near the center of sphere 2.
127  * - If sphere 1 is contained in sphere 2, then
128  * - \f$ R_2 \geq 1.01 (R_1 + |C_1-C_2|)\f$, where the 1.01
129  * prevents the spheres from (barely) touching.
130  * - If a line segment is drawn between \f$P\f$ and any point on the
131  * intersection circle (the circle where sphere 1 intersects the
132  * plane \f$z=z_\mathrm{P}\f$), the angle between the line segment
133  * and the z-axis is smaller than \f$\pi/3\f$.
134  */
136  public:
137  static constexpr size_t dim = 3;
138  CylindricalEndcap(const std::array<double, 3>& center_one,
139  const std::array<double, 3>& center_two,
140  const std::array<double, 3>& proj_center,
141  double radius_one, double radius_two,
142  double z_plane) noexcept;
143 
144  CylindricalEndcap() = default;
145  ~CylindricalEndcap() = default;
147  CylindricalEndcap(const CylindricalEndcap&) = default;
148  CylindricalEndcap& operator=(const CylindricalEndcap&) = default;
149  CylindricalEndcap& operator=(CylindricalEndcap&&) = default;
150 
151  template <typename T>
153  const std::array<T, 3>& source_coords) const noexcept;
154 
156  const std::array<double, 3>& target_coords) const noexcept;
157 
158  template <typename T>
159  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
160  const std::array<T, 3>& source_coords) const noexcept;
161 
162  template <typename T>
163  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
164  const std::array<T, 3>& source_coords) const noexcept;
165 
166  // clang-tidy: google runtime references
167  void pup(PUP::er& p) noexcept; // NOLINT
168 
169  static bool is_identity() noexcept { return false; }
170 
171  private:
172  friend bool operator==(const CylindricalEndcap& lhs,
173  const CylindricalEndcap& rhs) noexcept;
175 };
176 bool operator!=(const CylindricalEndcap& lhs,
177  const CylindricalEndcap& rhs) noexcept;
178 
179 } // namespace domain::CoordinateMaps
domain::CoordinateMaps::FocallyLiftedMap
Map from to the volume contained between a 2D surface and the surface of a 2-sphere.
Definition: FocallyLiftedMap.hpp:23
Frame::NoFrame
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
cstddef
array
domain::CoordinateMaps
Contains all coordinate maps.
Definition: Affine.cpp:14
limits
TypeAliases.hpp
optional
domain::CoordinateMaps::CylindricalEndcap
Map from 3D unit right cylinder to a volume that connects portions of two spherical surfaces.
Definition: CylindricalEndcap.hpp:135