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 <limits>
       9             : #include <optional>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Domain/CoordinateMaps/FocallyLiftedFlatEndcap.hpp"
      13             : #include "Domain/CoordinateMaps/FocallyLiftedMap.hpp"
      14             : 
      15             : /// \cond
      16             : namespace PUP {
      17             : class er;
      18             : }  // namespace PUP
      19             : /// \endcond
      20             : 
      21             : namespace domain::CoordinateMaps {
      22             : 
      23             : /*!
      24             :  * \ingroup CoordinateMapsGroup
      25             :  *
      26             :  * \brief Map from a 3D unit right cylinder to a volume that connects
      27             :  *  a flat circular disk (lying inside a sphere) to the far wall of
      28             :  *  that sphere.
      29             :  *
      30             :  * \details This is the "interior" counterpart to `CylindricalFlatEndcap`.
      31             :  * The two maps are identical in structure — both use `FocallyLiftedMap`
      32             :  * with `FocallyLiftedInnerMaps::FlatEndcap` — but differ in which
      33             :  * intersection of the projecting ray with the sphere is chosen:
      34             :  *
      35             :  * - `CylindricalFlatEndcap`: the flat disk lies *outside* the sphere
      36             :  *   (the sphere is between \f$P\f$ and the flat disk), so
      37             :  *   `source_is_between_focus_and_target = false`.
      38             :  * - `CylindricalFlatEndcapInterior`: the flat disk lies *inside* the
      39             :  *   sphere (the flat disk is between \f$P\f$ and the sphere's far wall),
      40             :  *   so `source_is_between_focus_and_target = true`.
      41             :  *
      42             :  * Consider a 2D circle in 3D space normal to the \f$z\f$ axis with
      43             :  * (3D) center \f$C_1\f$, a sphere with center \f$C_2\f$ and radius
      44             :  * \f$R_2\f$, and a projection point \f$P\f$.
      45             :  *
      46             :  * The parameter \f$z_\mathrm{extent}\f$ specifies the \f$z\f$-coordinate
      47             :  * (in the map's frame) of the rim circle where the spherical face of the
      48             :  * block meets the adjacent hollow-cylinder block.  This single number
      49             :  * determines the radius \f$R_1\f$ of the flat disk:
      50             :  *
      51             :  * \f{align}{
      52             :  *   t &= \frac{z_\mathrm{extent} - P_z}{C_1^z - P_z}, \\
      53             :  *   r_\mathrm{rim} &= \sqrt{R_2^2 - (z_\mathrm{extent} - C_2^z)^2}, \\
      54             :  *   R_1 &= \frac{r_\mathrm{rim}}{t}.
      55             :  * \f}
      56             :  *
      57             :  * CylindricalFlatEndcapInterior maps a 3D unit right cylinder
      58             :  * \f$(\bar{x},\bar{y},\bar{z})\f$ with \f$-1\leq\bar{z}\leq 1\f$ and
      59             :  * \f$\bar{x}^2+\bar{y}^2\leq 1\f$ so that:
      60             :  * - \f$\bar{z}=+1\f$ maps to the interior of the disk of radius
      61             :  *   \f$R_1\f$ centred at \f$C_1\f$.
      62             :  * - \f$\bar{z}=-1\f$ maps to the portion of the sphere on the *far* side
      63             :  *   of the flat disk from \f$P\f$.
      64             :  * - Curves of constant \f$(\bar{x},\bar{y})\f$ are portions of lines
      65             :  *   passing through \f$P\f$.
      66             :  * - The rim of the disk (\f$\bar{x}^2+\bar{y}^2=1\f$ on \f$\bar{z}=-1\f$)
      67             :  *   maps to the circle at \f$z = z_\mathrm{extent}\f$ on the sphere.
      68             :  *
      69             :  * Note that the \f$\bar{z}\f$ orientation is the *opposite* of
      70             :  * `CylindricalFlatEndcap`: here \f$\bar{z}=+1\f$ is the flat disk and
      71             :  * \f$\bar{z}=-1\f$ is the sphere.  This reversal is necessary to keep the
      72             :  * Jacobian determinant positive, because the flat disk sits at a larger
      73             :  * physical \f$z\f$ than the far sphere wall.
      74             :  *
      75             :  * CylindricalFlatEndcapInterior is intended for the Pill domain, where
      76             :  * the filled-cylinder endcap blocks have their flat face (\f$\bar{z}=+1\f$)
      77             :  * at the end of the inner cubed-cylinder region (inside the outer domain
      78             :  * sphere) and their spherical face (\f$\bar{z}=-1\f$) on the outer domain
      79             :  * sphere.
      80             :  *
      81             :  * ### Requirements on map parameters
      82             :  *
      83             :  * - \f$P\f$ is sufficiently inside the sphere:
      84             :  *   \f$|P - C_2| \leq 0.95\,R_2\f$.
      85             :  * - The flat disk lies inside the sphere, at least 5 % of \f$R_2\f$ below
      86             :  *   \f$C_2^z\f$ and at most 95 % of \f$R_2\f$ below \f$C_2^z\f$:
      87             :  *   \f[
      88             :  *     C_2^z - 0.95\,R_2 \;\leq\; C_1^z \;\leq\; C_2^z - 0.05\,R_2.
      89             :  *   \f]
      90             :  * - The flat disk is below the projection point: \f$C_1^z < P^z\f$.
      91             :  * - \f$z_\mathrm{extent}\f$ lies strictly on the sphere:
      92             :  *   \f$|z_\mathrm{extent} - C_2^z| < R_2\f$.
      93             :  * - The focal parameter satisfies \f$t > 1\f$ (the sphere is beyond the
      94             :  *   disk from \f$P\f$).
      95             :  * - The ratio \f$R_1/R_2\f$ is between 1/10 and 10.
      96             :  * - The entire flat disk rim lies strictly inside the sphere:
      97             :  *   \f[
      98             :  *     \sqrt{\bigl(\sqrt{(C_1^x-C_2^x)^2+(C_1^y-C_2^y)^2}+R_1\bigr)^2
      99             :  *           +(C_1^z-C_2^z)^2} < R_2.
     100             :  *   \f]
     101             :  *   This is required for the inverse to be defined everywhere: if any rim
     102             :  *   point were outside the sphere, the corresponding ray from \f$P\f$ would
     103             :  *   hit the sphere before reaching the disk (\f$t_\mathrm{sphere}<1\f$),
     104             :  *   violating the `source_is_between_focus_and_target` assumption.
     105             :  */
     106           1 : class CylindricalFlatEndcapInterior {
     107             :  public:
     108           0 :   static constexpr size_t dim = 3;
     109             : 
     110             :   /*!
     111             :    * \brief Construct the map.
     112             :    *
     113             :    * \param center_one Center of the flat disk (\f$C_1\f$).
     114             :    * \param center_two Center of the outer sphere (\f$C_2\f$).
     115             :    * \param proj_center Projection point \f$P\f$.
     116             :    * \param z_sphere_extent z-coordinate of the rim circle where the spherical
     117             :    *   face meets the adjacent hollow-cylinder block.  This determines the flat
     118             :    *   disk radius \f$R_1\f$ via the focal projection.
     119             :    * \param radius_two Radius of the outer sphere \f$R_2\f$.
     120             :    */
     121           1 :   CylindricalFlatEndcapInterior(const std::array<double, 3>& center_one,
     122             :                                 const std::array<double, 3>& center_two,
     123             :                                 const std::array<double, 3>& proj_center,
     124             :                                 double z_sphere_extent, double radius_two);
     125             : 
     126           0 :   CylindricalFlatEndcapInterior() = default;
     127           0 :   ~CylindricalFlatEndcapInterior() = default;
     128           0 :   CylindricalFlatEndcapInterior(CylindricalFlatEndcapInterior&&) = default;
     129           0 :   CylindricalFlatEndcapInterior(const CylindricalFlatEndcapInterior&) = default;
     130           0 :   CylindricalFlatEndcapInterior& operator=(
     131             :       const CylindricalFlatEndcapInterior&) = default;
     132           0 :   CylindricalFlatEndcapInterior& operator=(CylindricalFlatEndcapInterior&&) =
     133             :       default;
     134             : 
     135             :   template <typename T>
     136           0 :   std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
     137             : 
     138           0 :   std::optional<std::array<double, 3>> inverse(
     139             :       const std::array<double, 3>& target_coords) const;
     140             : 
     141             :   template <typename T>
     142           0 :   tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
     143             :       const std::array<T, 3>& source_coords) const;
     144             : 
     145             :   template <typename T>
     146           0 :   tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
     147             :       const std::array<T, 3>& source_coords) const;
     148             : 
     149             :   // NOLINTNEXTLINE(google-runtime-references)
     150           0 :   void pup(PUP::er& p);
     151             : 
     152           0 :   static bool is_identity() { return false; }
     153             : 
     154           0 :   static constexpr bool supports_hessian{false};
     155             : 
     156             :  private:
     157           0 :   friend bool operator==(const CylindricalFlatEndcapInterior& lhs,
     158             :                          const CylindricalFlatEndcapInterior& rhs);
     159           0 :   FocallyLiftedMap<FocallyLiftedInnerMaps::FlatEndcap> impl_;
     160             : };
     161             : 
     162           0 : bool operator!=(const CylindricalFlatEndcapInterior& lhs,
     163             :                 const CylindricalFlatEndcapInterior& rhs);
     164             : 
     165             : }  // namespace domain::CoordinateMaps