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