Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines the class CylindricalFlatEndcap.
       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/FocallyLiftedFlatEndcap.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             :  *  a portion of a circle to a portion of a spherical surface.
      31             :  *
      32             :  * \image html CylindricalFlatEndcap.svg "A cylinder maps to the shaded region."
      33             :  *
      34             :  * \details Consider a 2D circle in 3D space that is normal to the
      35             :  * \f$z\f$ axis and has (3D) center \f$C_1\f$ and radius \f$R_1\f$.
      36             :  * Also consider a sphere with center \f$C_2\f$, and radius \f$R_2\f$.
      37             :  * Also let there be a projection point \f$P\f$.
      38             :  *
      39             :  * CylindricalFlatEndcap 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 interior of the
      44             :  * circle of radius \f$R_1\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 the circle and \f$\bar{z}=+1\f$ is mapped to a
      48             :  * point on the sphere.
      49             :  *
      50             :  * CylindricalFlatEndcap is intended to be composed with Wedge2D maps to
      51             :  * construct a portion of a cylindrical domain for a binary system.
      52             :  *
      53             :  * CylindricalFlatEndcap is described briefly in the Appendix of
      54             :  * \cite Buchman:2012dw.
      55             :  * CylindricalFlatEndcap is used to construct the blocks labeled 'MA
      56             :  * wedge' and 'MB wedge' in Figure 20 of that paper.
      57             :  *
      58             :  * CylindricalFlatEndcap is implemented using `FocallyLiftedMap`
      59             :  * and `FocallyLiftedInnerMaps::FlatEndcap`; see those classes for
      60             :  * details.
      61             :  *
      62             :  * ### Restrictions on map parameters.
      63             :  *
      64             :  * The following restrictions are made so that the map is not singular
      65             :  * or close to singular. It is possible to construct a valid map
      66             :  * without these assumptions, but the assumptions simplify the code and
      67             :  * avoid problematic edge cases, and the expected use cases obey
      68             :  * these restrictions.
      69             :  *
      70             :  * We demand that
      71             :  * - The plane containing the circle is below (i.e. at a smaller value
      72             :  *   of \f$z\f$ than) the sphere, by an amount at least 5% of the sphere
      73             :  *   radius \f$R_2\f$ but not more than 5 times the sphere radius \f$R_2\f$.
      74             :  * - \f$P\f$ is inside the sphere but not too close to its surface;
      75             :  *   specifically, we demand that \f$|P-C_2|\leq 0.95 R_2\f$.
      76             :  * - The ratio \f$R_1/R_2\f$ is between 10 and 1/10, inclusive.
      77             :  * - The x and y components of \f$C_2-C_1\f$ both have magnitudes
      78             :  *   smaller than or equal to \f$R_1+R_2\f$.
      79             :  *
      80             :  */
      81           1 : class CylindricalFlatEndcap {
      82             :  public:
      83           0 :   static constexpr size_t dim = 3;
      84           0 :   CylindricalFlatEndcap(const std::array<double, 3>& center_one,
      85             :                         const std::array<double, 3>& center_two,
      86             :                         const std::array<double, 3>& proj_center,
      87             :                         double radius_one, double radius_two);
      88             : 
      89           0 :   CylindricalFlatEndcap() = default;
      90           0 :   ~CylindricalFlatEndcap() = default;
      91           0 :   CylindricalFlatEndcap(CylindricalFlatEndcap&&) = default;
      92           0 :   CylindricalFlatEndcap(const CylindricalFlatEndcap&) = default;
      93           0 :   CylindricalFlatEndcap& operator=(const CylindricalFlatEndcap&) = default;
      94           0 :   CylindricalFlatEndcap& operator=(CylindricalFlatEndcap&&) = default;
      95             : 
      96             :   template <typename T>
      97           0 :   std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
      98             : 
      99           0 :   std::optional<std::array<double, 3>> inverse(
     100             :       const std::array<double, 3>& target_coords) const;
     101             : 
     102             :   template <typename T>
     103           0 :   tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
     104             :       const std::array<T, 3>& source_coords) const;
     105             : 
     106             :   template <typename T>
     107           0 :   tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
     108             :       const std::array<T, 3>& source_coords) const;
     109             : 
     110             :   // NOLINTNEXTLINE(google-runtime-references)
     111           0 :   void pup(PUP::er& p);
     112             : 
     113           0 :   static bool is_identity() { return false; }
     114             : 
     115           0 :   static constexpr bool supports_hessian{false};
     116             : 
     117             :  private:
     118           0 :   friend bool operator==(const CylindricalFlatEndcap& lhs,
     119             :                          const CylindricalFlatEndcap& rhs);
     120           0 :   FocallyLiftedMap<FocallyLiftedInnerMaps::FlatEndcap> impl_;
     121             : };
     122           0 : bool operator!=(const CylindricalFlatEndcap& lhs,
     123             :                 const CylindricalFlatEndcap& rhs);
     124             : 
     125             : }  // namespace domain::CoordinateMaps