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