Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines the class CylindricalSide.
       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/FocallyLiftedMap.hpp"
      16             : #include "Domain/CoordinateMaps/FocallyLiftedSide.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 a 3D unit right cylindrical shell to a volume that connects
      31             :  *  portions of two spherical surfaces.
      32             :  *
      33             :  * \image html CylindricalSide.svg "2D slice showing mapped (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 two
      37             :  * planes normal to the \f$z\f$ axis and located at \f$z = z_\mathrm{L}\f$
      38             :  * and \f$z = z_\mathrm{U}\f$, with \f$z_\mathrm{L} < z_\mathrm{U}\f$.
      39             :  * Also let there be a projection point \f$P\f$.
      40             :  *
      41             :  * Note that Sphere 1 and Sphere 2 are not equivalent, because the
      42             :  * mapped portion of Sphere 1 is bounded by planes of constant
      43             :  * \f$z\f$ but the mapped portion of Sphere 2 is not (except for
      44             :  * special choices of \f$C_1\f$, \f$C_2\f$, and \f$P\f$).
      45             :  *
      46             :  * CylindricalSide maps a 3D unit right cylindrical shell (with
      47             :  * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
      48             :  * \f$-1\leq\bar{z}\leq 1\f$ and \f$1 \leq \bar{x}^2+\bar{y}^2 \leq
      49             :  * 4\f$) to the shaded area in each panel of the figure above (with
      50             :  * coordinates \f$(x,y,z)\f$).  The figure shows two different allowed
      51             :  * possibilities: \f$R_1 > R_2\f$ and \f$R_2 > R_1\f$. Note that the
      52             :  * two portions of the shaded region in each panel of the figure
      53             :  * represent different portions of the same block; each panel of the
      54             :  * figure is to be understood as rotated around the \f$z\f$ axis. The
      55             :  * inner boundary of the cylindrical shell \f$\bar{x}^2+\bar{y}^2=1\f$
      56             :  * is mapped to the portion of sphere 1 that has \f$z_\mathrm{L} \leq
      57             :  * z \leq z_\mathrm{U}\f$.  Curves of constant \f$(\bar{z})\f$ along
      58             :  * the vector \f$(\bar{x},\bar{y})\f$ are mapped to portions of lines
      59             :  * that pass through \f$P\f$. Along each of these curves,
      60             :  * \f$\bar{x}^2+\bar{y}^2=1\f$ is mapped to a point on sphere 1 and
      61             :  * \f$\bar{x}^2+\bar{y}^2=4\f$ is mapped to a point on sphere 2.
      62             :  *
      63             :  * CylindricalSide is described briefly in the Appendix of
      64             :  * \cite Buchman:2012dw.  CylindricalSide is used to construct the blocks
      65             :  * labeled 'CA cylinder', 'EA cylinder', 'CB cylinder', 'EE cylinder',
      66             :  * and 'EB cylinder' in Figure 20 of that paper.  Note that 'CA
      67             :  * cylinder', 'CB cylinder', and 'EE cylinder' have Sphere 1 contained
      68             :  * in Sphere2, and 'EA cylinder' and 'EB cylinder' have Sphere 2
      69             :  * contained in Sphere 1.
      70             :  *
      71             :  * CylindricalSide is implemented using `FocallyLiftedMap`
      72             :  * and `FocallyLiftedInnerMaps::Side`; see those classes for
      73             :  * details.
      74             :  *
      75             :  * ### Restrictions on map parameters.
      76             :  *
      77             :  * We demand that:
      78             :  * - Either Sphere 1 is fully contained inside Sphere 2, or
      79             :  *   Sphere 2 is fully contained inside Sphere 1.
      80             :  * - \f$P\f$ is contained inside the smaller sphere, and
      81             :  *   between (or on) the planes defined by \f$z_\mathrm{L}\f$ and
      82             :  *   \f$z_\mathrm{U}\f$.
      83             :  * - If sphere 1 is contained in sphere 2:
      84             :  *   - \f$C_1^z - 0.95 R_1 \leq z_\mathrm{L}\f$
      85             :  *   - \f$z_\mathrm{U} \leq C_1^z + 0.95 R_1\f$
      86             :  * - If sphere 2 is contained in sphere 1:
      87             :  *   - \f$C_1^z - 0.95 R_1 \leq z_\mathrm{L} \leq C_1^z - 0.2 R_1\f$
      88             :  *   - \f$C_1^z + 0.2 R_1 \leq z_\mathrm{U} \leq C_1^z + 0.95 R_1\f$
      89             :  *
      90             :  */
      91           1 : class CylindricalSide {
      92             :  public:
      93           0 :   static constexpr size_t dim = 3;
      94           0 :   CylindricalSide(const std::array<double, 3>& center_one,
      95             :                   const std::array<double, 3>& center_two,
      96             :                   const std::array<double, 3>& proj_center, double radius_one,
      97             :                   double radius_two, double z_lower, double z_upper);
      98             : 
      99           0 :   CylindricalSide() = default;
     100           0 :   ~CylindricalSide() = default;
     101           0 :   CylindricalSide(CylindricalSide&&) = default;
     102           0 :   CylindricalSide(const CylindricalSide&) = default;
     103           0 :   CylindricalSide& operator=(const CylindricalSide&) = default;
     104           0 :   CylindricalSide& operator=(CylindricalSide&&) = default;
     105             : 
     106             :   template <typename T>
     107           0 :   std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
     108             :       const std::array<T, 3>& source_coords) const;
     109             : 
     110           0 :   std::optional<std::array<double, 3>> inverse(
     111             :       const std::array<double, 3>& target_coords) const;
     112             : 
     113             :   template <typename T>
     114           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
     115             :       const std::array<T, 3>& source_coords) const;
     116             : 
     117             :   template <typename T>
     118           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
     119             :       const std::array<T, 3>& source_coords) const;
     120             : 
     121             :   // NOLINTNEXTLINE(google-runtime-references)
     122           0 :   void pup(PUP::er& p);
     123             : 
     124           0 :   static bool is_identity() { return false; }
     125             : 
     126             :  private:
     127           0 :   friend bool operator==(const CylindricalSide& lhs,
     128             :                          const CylindricalSide& rhs);
     129           0 :   FocallyLiftedMap<FocallyLiftedInnerMaps::Side> impl_;
     130             : };
     131           0 : bool operator!=(const CylindricalSide& lhs, const CylindricalSide& rhs);
     132             : 
     133             : }  // namespace domain::CoordinateMaps