Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines the class CylindricalFlatSide.
       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/FocallyLiftedFlatSide.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 cylindrical shell to a volume that connects
      30             :  *  a portion of an annulus to a portion of a spherical surface.
      31             :  *
      32             :  * \image html CylindricalFlatSide.svg "A cylinder maps to the shaded region."
      33             :  *
      34             :  * \details Consider a 2D annulus in 3D space that is normal to the
      35             :  * \f$z\f$ axis and has (3D) center \f$C_1\f$, inner radius
      36             :  * \f$R_\mathrm{in}\f$ and outer radius \f$R_\mathrm{out}\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             :  * CylindricalFlatSide maps a 3D unit right cylindrical shell (with
      41             :  * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
      42             :  * \f$-1\leq\bar{z}\leq 1\f$ and \f$1 \leq \bar{x}^2+\bar{y}^2 \leq
      43             :  * 4\f$) to the shaded area in the figure above (with coordinates
      44             :  * \f$(x,y,z)\f$).  The "bottom" of the cylinder \f$\bar{z}=-1\f$ is
      45             :  * mapped to the interior of the annulus with radii
      46             :  * \f$R_\mathrm{in}\f$ and \f$R_\mathrm{out}\f$.  Curves of constant
      47             :  * \f$(\bar{x},\bar{y})\f$ are mapped to portions of lines that pass
      48             :  * through \f$P\f$. Along each of these curves, \f$\bar{z}=-1\f$ is
      49             :  * mapped to a point inside the annulus and \f$\bar{z}=+1\f$ is mapped to a
      50             :  * point on the sphere.
      51             :  *
      52             :  * CylindricalFlatSide is intended to be composed with Wedge2D maps to
      53             :  * construct a portion of a cylindrical domain for a binary system.
      54             :  *
      55             :  * CylindricalFlatSide is described briefly in the Appendix of
      56             :  * \cite Buchman:2012dw.
      57             :  * CylindricalFlatSide is used to construct the blocks labeled 'ME
      58             :  * cylinder' in Figure 20 of that paper.
      59             :  *
      60             :  * CylindricalFlatSide is implemented using `FocallyLiftedMap`
      61             :  * and `FocallyLiftedInnerMaps::FlatSide`; see those classes for
      62             :  * details.
      63             :  *
      64             :  * ### Restrictions on map parameters.
      65             :  *
      66             :  * We demand that:
      67             :  * - The sphere is at a larger value of \f$z\f$ (plus 5
      68             :  *   percent of the sphere radius) than the plane containing the
      69             :  *   annulus.
      70             :  * - The projection point \f$z_\mathrm{P}\f$ is
      71             :  *   inside the sphere and more than 15 percent away from the boundary
      72             :  *   of the sphere.
      73             :  * - The center of the annulus is contained in the circle that results from
      74             :  *   projecting the sphere into the \f$xy\f$ plane.
      75             :  * - The outer radius of the annulus is larger than 5 percent of the distance
      76             :  *   between the center of the annulus and the projection point.
      77             :  * - The inner radius of the annulus is less than 95 percent of the outer
      78             :  *   radius, larger than 5 percent of the outer radius, and larger than one
      79             :  *   percent of the distance between the center of the annulus and the
      80             :  *   projection point.  The last condition means that the angle subtended by
      81             :  *   the inner radius with respect to the projection point is not too small.
      82             :  *
      83             :  * It is possible to construct a valid map without these assumptions,
      84             :  * but some of these assumptions simplify the code and others eliminate
      85             :  * edge cases where Jacobians become large or small.
      86             :  *
      87             :  */
      88           1 : class CylindricalFlatSide {
      89             :  public:
      90           0 :   static constexpr size_t dim = 3;
      91           0 :   CylindricalFlatSide(const std::array<double, 3>& center_one,
      92             :                       const std::array<double, 3>& center_two,
      93             :                       const std::array<double, 3>& proj_center,
      94             :                       const double inner_radius, const double outer_radius,
      95             :                       const double radius_two);
      96             : 
      97           0 :   CylindricalFlatSide() = default;
      98           0 :   ~CylindricalFlatSide() = default;
      99           0 :   CylindricalFlatSide(CylindricalFlatSide&&) = default;
     100           0 :   CylindricalFlatSide(const CylindricalFlatSide&) = default;
     101           0 :   CylindricalFlatSide& operator=(const CylindricalFlatSide&) = default;
     102           0 :   CylindricalFlatSide& operator=(CylindricalFlatSide&&) = default;
     103             : 
     104             :   template <typename T>
     105           0 :   std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
     106             : 
     107           0 :   std::optional<std::array<double, 3>> inverse(
     108             :       const std::array<double, 3>& target_coords) const;
     109             : 
     110             :   template <typename T>
     111           0 :   tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
     112             :       const std::array<T, 3>& source_coords) const;
     113             : 
     114             :   template <typename T>
     115           0 :   tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
     116             :       const std::array<T, 3>& source_coords) const;
     117             : 
     118             :   // NOLINTNEXTLINE(google-runtime-references)
     119           0 :   void pup(PUP::er& p);
     120             : 
     121           0 :   static bool is_identity() { return false; }
     122             : 
     123           0 :   static constexpr bool supports_hessian{false};
     124             : 
     125             :  private:
     126           0 :   friend bool operator==(const CylindricalFlatSide& lhs,
     127             :                          const CylindricalFlatSide& rhs);
     128           0 :   FocallyLiftedMap<FocallyLiftedInnerMaps::FlatSide> impl_;
     129             : };
     130           0 : bool operator!=(const CylindricalFlatSide& lhs, const CylindricalFlatSide& rhs);
     131             : 
     132             : }  // namespace domain::CoordinateMaps