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 : #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 `Wedge<2>` 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.15 \leq \cos(\theta) \leq 0.95\f$, where 120 : * the values 0.15 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.1 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 : */ 135 1 : class CylindricalEndcap { 136 : public: 137 0 : static constexpr size_t dim = 3; 138 0 : CylindricalEndcap(const std::array<double, 3>& center_one, 139 : const std::array<double, 3>& center_two, 140 : const std::array<double, 3>& proj_center, double radius_one, 141 : double radius_two, double z_plane); 142 : 143 0 : CylindricalEndcap() = default; 144 0 : ~CylindricalEndcap() = default; 145 0 : CylindricalEndcap(CylindricalEndcap&&) = default; 146 0 : CylindricalEndcap(const CylindricalEndcap&) = default; 147 0 : CylindricalEndcap& operator=(const CylindricalEndcap&) = default; 148 0 : CylindricalEndcap& operator=(CylindricalEndcap&&) = default; 149 : 150 : template <typename T> 151 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()( 152 : const std::array<T, 3>& source_coords) const; 153 : 154 : /// The inverse function is only callable with doubles because the inverse 155 : /// might fail if called for a point out of range, and it is unclear 156 : /// what should happen if the inverse were to succeed for some points in a 157 : /// DataVector but fail for other points. 158 1 : std::optional<std::array<double, 3>> inverse( 159 : const std::array<double, 3>& target_coords) const; 160 : 161 : template <typename T> 162 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian( 163 : const std::array<T, 3>& source_coords) const; 164 : 165 : template <typename T> 166 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian( 167 : const std::array<T, 3>& source_coords) const; 168 : 169 : // NOLINTNEXTLINE(google-runtime-references) 170 0 : void pup(PUP::er& p); 171 : 172 0 : static bool is_identity() { return 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