Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <optional>
       9             : 
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : 
      12             : /// \cond
      13             : namespace PUP {
      14             : class er;
      15             : }  // namespace PUP
      16             : /// \endcond
      17             : 
      18             : namespace domain::CoordinateMaps {
      19             : 
      20             : /*!
      21             :  * \ingroup CoordinateMapsGroup
      22             :  *
      23             :  * \brief Maps a cylindrical shell block to a region bounded by an inner
      24             :  * right cylinder and an outer spherical surface.
      25             :  *
      26             :  * \image html CylindricalSphericalShell.png "The shaded region is the image."
      27             :  *
      28             :  * \details The logical coordinates are:
      29             :  * - radial \f$\xi \in [-1,1]\f$: interpolates between the inner right
      30             :  *   cylinder and the outer sphere.
      31             :  * - angular \f$\eta \in (-\pi, \pi]\f$: azimuthal angle \f$\phi\f$ (S1
      32             :  *   periodic direction).  The inverse map returns
      33             :  *   \f$\eta = \mathrm{atan2}(z, y)\f$, which lies in \f$(-\pi, \pi]\f$.
      34             :  * - axial \f$\zeta \in [-1,1]\f$: blends between the lower and upper ends
      35             :  *   of the shell.
      36             :  *
      37             :  * The physical coordinates \f$(x, y, z)\f$ are computed as follows.
      38             :  * Let \f$\alpha = (\xi+1)/2\f$ and \f$\beta = (\zeta+1)/2\f$.  Define
      39             :  *
      40             :  * \f{align}{
      41             :  *   x_\mathrm{inner}(\beta) &= x^\mathrm{inner}_\mathrm{lower}
      42             :  *       + \beta\,(x^\mathrm{inner}_\mathrm{upper}
      43             :  *                - x^\mathrm{inner}_\mathrm{lower}), \\
      44             :  *   x_\mathrm{outer}(\beta) &= x^\mathrm{outer}_\mathrm{lower}
      45             :  *       + \beta\,(x^\mathrm{outer}_\mathrm{upper}
      46             :  *                - x^\mathrm{outer}_\mathrm{lower}), \\
      47             :  *   r_\mathrm{outer}(\beta) &= \sqrt{r_\mathrm{sphere}^2
      48             :  *       - x_\mathrm{outer}(\beta)^2}.
      49             :  * \f}
      50             :  *
      51             :  * Then
      52             :  *
      53             :  * \f{align}{
      54             :  *   x &= (1-\alpha)\,x_\mathrm{inner}(\beta)
      55             :  *         + \alpha\,x_\mathrm{outer}(\beta), \\
      56             :  *   r &= (1-\alpha)\,r_\mathrm{inner}
      57             :  *         + \alpha\,r_\mathrm{outer}(\beta), \\
      58             :  *   y &= r\cos\eta, \quad z = r\sin\eta.
      59             :  * \f}
      60             :  *
      61             :  * The six block faces have the following geometry:
      62             :  * - \f$\xi = -1\f$ (\f$\alpha=0\f$): the inner right cylinder at
      63             :  *   \f$r = r_\mathrm{inner}\f$, extending in \f$x\f$ from
      64             :  *   \f$x^\mathrm{inner}_\mathrm{lower}\f$ to
      65             :  *   \f$x^\mathrm{inner}_\mathrm{upper}\f$.
      66             :  * - \f$\xi = +1\f$ (\f$\alpha=1\f$): a portion of the sphere
      67             :  *   \f$r^2 + x^2 = r_\mathrm{sphere}^2\f$.
      68             :  * - \f$\zeta = \pm 1\f$: generically curved surfaces in Cartesian
      69             :  *   coordinates (ruled surfaces that blend linearly in \f$(x,r)\f$ between
      70             :  *   the inner cylinder edge and the outer sphere edge at the corresponding
      71             :  *   axial end).
      72             :  * - \f$\eta\f$ direction: periodic (connected by the cylindrical_shell
      73             :  *   topology).
      74             :  *
      75             :  * ### Requirements
      76             :  * - \f$x^\mathrm{inner}_\mathrm{lower} < x^\mathrm{inner}_\mathrm{upper}\f$
      77             :  * - \f$x^\mathrm{outer}_\mathrm{lower} < x^\mathrm{outer}_\mathrm{upper}\f$
      78             :  * - \f$r_\mathrm{inner} > 0\f$
      79             :  * - \f$|x^\mathrm{outer}_\mathrm{lower}|,
      80             :  *   |x^\mathrm{outer}_\mathrm{upper}| < r_\mathrm{sphere}\f$
      81             :  * - \f$r_\mathrm{inner} < r_\mathrm{outer}(\beta)\f$ for all
      82             :  *   \f$\beta \in [0,1]\f$ (equivalently,
      83             :  *   \f$r_\mathrm{inner} < \min(r_\mathrm{outer}(0), r_\mathrm{outer}(1))\f$).
      84             :  */
      85           1 : class CylindricalSphericalShell {
      86             :  public:
      87           0 :   static constexpr size_t dim = 3;
      88             : 
      89             :   /*!
      90             :    * \brief Construct the map.
      91             :    *
      92             :    * \param x_inner_lower Axial coordinate \f$x\f$ at the lower end of the
      93             :    *   inner cylinder (\f$\xi=-1, \zeta=-1\f$).
      94             :    * \param x_inner_upper Axial coordinate \f$x\f$ at the upper end of the
      95             :    *   inner cylinder (\f$\xi=-1, \zeta=+1\f$).
      96             :    * \param x_outer_lower Axial coordinate \f$x\f$ at the lower end of the
      97             :    *   outer spherical face (\f$\xi=+1, \zeta=-1\f$).
      98             :    * \param x_outer_upper Axial coordinate \f$x\f$ at the upper end of the
      99             :    *   outer spherical face (\f$\xi=+1, \zeta=+1\f$).
     100             :    * \param r_inner Radius of the inner right cylinder.
     101             :    * \param r_sphere Radius of the outer bounding sphere (centered at the
     102             :    *   origin).
     103             :    */
     104           1 :   CylindricalSphericalShell(double x_inner_lower, double x_inner_upper,
     105             :                             double x_outer_lower, double x_outer_upper,
     106             :                             double r_inner, double r_sphere);
     107             : 
     108           0 :   CylindricalSphericalShell() = default;
     109           0 :   ~CylindricalSphericalShell() = default;
     110           0 :   CylindricalSphericalShell(CylindricalSphericalShell&&) = default;
     111           0 :   CylindricalSphericalShell(const CylindricalSphericalShell&) = default;
     112           0 :   CylindricalSphericalShell& operator=(const CylindricalSphericalShell&) =
     113             :       default;
     114           0 :   CylindricalSphericalShell& operator=(CylindricalSphericalShell&&) = default;
     115             : 
     116             :   template <typename T>
     117           0 :   std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
     118             : 
     119           0 :   std::optional<std::array<double, 3>> inverse(
     120             :       const std::array<double, 3>& target_coords) const;
     121             : 
     122             :   template <typename T>
     123           0 :   tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
     124             :       const std::array<T, 3>& source_coords) const;
     125             : 
     126             :   template <typename T>
     127           0 :   tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
     128             :       const std::array<T, 3>& source_coords) const;
     129             : 
     130             :   // NOLINTNEXTLINE(google-runtime-references)
     131           0 :   void pup(PUP::er& p);
     132             : 
     133           0 :   static bool is_identity() { return false; }
     134             : 
     135           0 :   static constexpr bool supports_hessian{false};
     136             : 
     137             :  private:
     138           0 :   friend bool operator==(const CylindricalSphericalShell& lhs,
     139             :                          const CylindricalSphericalShell& rhs);
     140           0 :   double x_inner_lower_{};
     141           0 :   double x_inner_upper_{};
     142           0 :   double x_outer_lower_{};
     143           0 :   double x_outer_upper_{};
     144           0 :   double r_inner_{};
     145           0 :   double r_sphere_{};
     146             : };
     147             : 
     148           0 : bool operator!=(const CylindricalSphericalShell& lhs,
     149             :                 const CylindricalSphericalShell& rhs);
     150             : 
     151             : }  // namespace domain::CoordinateMaps