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