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