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