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