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 "Utilities/Gsl.hpp"
      13             : 
      14             : /// \cond
      15             : namespace PUP {
      16             : class er;
      17             : }  // namespace PUP
      18             : /// \endcond
      19             : 
      20             : /// Contains FocallyLiftedInnerMaps
      21             : namespace domain::CoordinateMaps::FocallyLiftedInnerMaps {
      22             : /*!
      23             :  * \brief A FocallyLiftedInnerMap that maps a 3D unit right cylinder
      24             :  *  to a volume that connects a 2D annulus to a spherical
      25             :  *  surface.
      26             :  *
      27             :  * \details The domain of the map is a 3D unit right cylinder with
      28             :  * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
      29             :  * \f$-1\leq\bar{z}\leq 1\f$ and \f$1\leq \bar{x}^2+\bar{y}^2 \leq
      30             :  * 4\f$.  The range of the map has coordinates \f$(x,y,z)\f$.
      31             :  *
      32             :  * Consider a 2D annulus in 3D space
      33             :  * oriented normal to the \f$z\f$ axis.  The inner and outer radii
      34             :  * of the annulus are \f$R_\mathrm{in}\f$ and \f$R_\mathrm{out}\f$, and the
      35             :  * (3D) center of the annulus is \f$C^i\f$.
      36             :  * `FlatSide` provides the following functions:
      37             :  *
      38             :  * ### forward_map()
      39             :  * `forward_map()` maps \f$(\bar{x},\bar{y},\bar{z}=-1)\f$ to the interior
      40             :  * of the annulus.  The arguments to `forward_map()`
      41             :  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
      42             :  * `forward_map()` returns \f$x_0^i\f$,
      43             :  * the 3D coordinates on the annulus, which are given by
      44             :  *
      45             :  * \f{align}
      46             :  * x_0^0 &= \left(R_\mathrm{in}+(R_\mathrm{out}-R_\mathrm{in})
      47             :  *    (\bar{\rho}-1)\right)
      48             :  *    \frac{\bar{x}}{\bar{\rho}} + C^0, \label{eq:forward_map_x}\\
      49             :  * x_0^1 &= \left(R_\mathrm{in}+(R_\mathrm{out}-R_\mathrm{in})
      50             :  *    (\bar{\rho}-1)\right)
      51             :  *    \frac{\bar{y}}{\bar{\rho}} + C^1,\\
      52             :  * x_0^2 &= C^2 \label{eq:forward_map_z},
      53             :  * \f}
      54             :  *
      55             :  * where
      56             :  *
      57             :  * \f{align} \bar{\rho} = \sqrt{\bar{x}^2+\bar{y}^2}.\label{eq:rhobar}\f}
      58             :  *
      59             :  * ### sigma
      60             :  *
      61             :  * \f$\sigma\f$ is a function that is zero on the sphere
      62             :  * \f$x^i=x_0^i\f$ and unity at \f$\bar{z}=+1\f$ (corresponding to the
      63             :  * upper surface of the FocallyLiftedMap). We define
      64             :  *
      65             :  * \f{align}
      66             :  *  \sigma &= \frac{\bar{z}+1}{2}.
      67             :  * \f}
      68             :  *
      69             :  * ### deriv_sigma
      70             :  *
      71             :  * `deriv_sigma` returns
      72             :  *
      73             :  * \f{align}
      74             :  *  \frac{\partial \sigma}{\partial \bar{x}^j} &= (0,0,1/2).
      75             :  *  \label{eq:deriv_sigma}
      76             :  * \f}
      77             :  *
      78             :  * ### jacobian
      79             :  *
      80             :  * `jacobian` returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
      81             :  * The arguments to `jacobian`
      82             :  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
      83             :  *
      84             :  * Differentiating
      85             :  * Eqs. (\f$\ref{eq:forward_map_x}\f$--\f$\ref{eq:forward_map_z}\f$)
      86             :  * above yields
      87             :  *
      88             :  * \f{align*}
      89             :  * \frac{\partial x_0^0}{\partial \bar{x}} &=
      90             :  * R_\mathrm{out}-R_\mathrm{in} + (2 R_\mathrm{in}-R_\mathrm{out})
      91             :  * \frac{\bar{y}^2}{\bar{\rho}^3},\\
      92             :  * \frac{\partial x_0^0}{\partial \bar{y}} &=
      93             :  * -(2 R_\mathrm{in}-R_\mathrm{out})
      94             :  * \frac{\bar{x}\bar{y}}{\bar{\rho}^3},\\
      95             :  * \frac{\partial x_0^1}{\partial \bar{x}} &=
      96             :  * -(2 R_\mathrm{in}-R_\mathrm{out})
      97             :  * \frac{\bar{x}\bar{y}}{\bar{\rho}^3},\\
      98             :  * \frac{\partial x_0^1}{\partial \bar{y}} &=
      99             :  * R_\mathrm{out}-R_\mathrm{in} + (2 R_\mathrm{in}-R_\mathrm{out})
     100             :  * \frac{\bar{x}^2}{\bar{\rho}^3},\\
     101             :  * \f}
     102             :  * and all other components are zero.
     103             :  *
     104             :  * ### inverse
     105             :  *
     106             :  * `inverse` takes \f$x_0^i\f$ and \f$\sigma\f$ as arguments, and
     107             :  * returns \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed
     108             :  * `std::optional<std::array<double, 3>>` if
     109             :  * \f$x_0^i\f$ or \f$\sigma\f$ are outside the range of the map.
     110             :  *
     111             :  * Let
     112             :  * \f{align}
     113             :  *  \rho = \sqrt{(x_0^0-C^0)^2+(x_0^1-C^1)^2}.
     114             :  *  \label{eq:rho}
     115             :  * \f}
     116             :  *
     117             :  * Then
     118             :  * \f{align}
     119             :  *  \bar{x} &= \frac{x_0^0-C^0}{\rho}
     120             :  *  \frac{\rho+R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
     121             :  *  \bar{y} &= \frac{x_0^1-C^1}{\rho}
     122             :  *  \frac{\rho+R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
     123             :  *  \bar{z} &= 2\sigma - 1.
     124             :  * \f}
     125             :  *
     126             :  * Note that \f$\rho\f$ in Eq. (\f$\ref{eq:rho}\f$) can be written
     127             :  * \f{align}
     128             :  * \rho = R_\mathrm{in}+(R_\mathrm{out}-R_\mathrm{in})(\bar{\rho}-1),
     129             :  * \label{eq:rho_from_rhobar}
     130             :  * \f}
     131             :  * where \f$\bar{\rho}\f$ is given by Eq. (\f$\ref{eq:rhobar}\f$).
     132             :  *
     133             :  * If \f$\bar{z}\f$ is outside the range \f$[-1,1]\f$ or
     134             :  * if \f$\bar{x}^2+\bar{y}^2\f$ is less than 1 or greater than 4
     135             :  * then we return a default-constructed
     136             :  * `std::optional<std::array<double, 3>>`.
     137             :  *
     138             :  * ### lambda_tilde
     139             :  *
     140             :  * `lambda_tilde` takes as arguments a point \f$x^i\f$ and a projection point
     141             :  *  \f$P^i\f$, and computes \f$\tilde{\lambda}\f$, the solution to
     142             :  *
     143             :  * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
     144             :  *
     145             :  * Since \f$x_0^i\f$ must lie on the plane \f$x_0^3=C^3\f$,
     146             :  *
     147             :  * \f{align} \tilde{\lambda} &= \frac{C^3-P^3}{x^3-P^3}.\f}
     148             :  *
     149             :  * If \f$\tilde{\lambda}\f$ is less than unity (indicating that the
     150             :  * supplied point is outside the range of the map), then a
     151             :  * default-constructed `std::optional<double>` is returned.
     152             :  *
     153             :  * ### deriv_lambda_tilde
     154             :  *
     155             :  * `deriv_lambda_tilde` takes as arguments \f$x_0^i\f$, a projection point
     156             :  *  \f$P^i\f$, and \f$\tilde{\lambda}\f$, and
     157             :  *  returns \f$\partial \tilde{\lambda}/\partial x^i\f$. We have
     158             :  *
     159             :  * \f{align}
     160             :  * \frac{\partial\tilde{\lambda}}{\partial x^3} =
     161             :  * -\frac{C^3-P^3}{(x^3-P^3)^2} = -\frac{\tilde{\lambda}^2}{C^3-P^3},
     162             :  * \f}
     163             :  * and other components are zero.
     164             :  *
     165             :  * ### inv_jacobian
     166             :  *
     167             :  * `inv_jacobian` returns \f$\partial \bar{x}^i/\partial x_0^k\f$,
     168             :  *  where \f$\sigma\f$ is held fixed.
     169             :  *  The arguments to `inv_jacobian`
     170             :  *  are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
     171             :  *
     172             :  * The nonzero components are
     173             :  * \f{align}
     174             :  * \frac{\partial \bar{x}}{\partial x_0^0} &=
     175             :  *  \frac{1}{R_\mathrm{out}-R_\mathrm{in}} + \frac{\bar{y}^2}{\bar{\rho}^2\rho}
     176             :  *  \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
     177             :  * \frac{\partial \bar{x}}{\partial x_0^1} &=
     178             :  *  - \frac{\bar{x}\bar{y}}{\bar{\rho}^2\rho}
     179             :  *  \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
     180             :  * \frac{\partial \bar{y}}{\partial x_0^1} &=
     181             :  *  \frac{1}{R_\mathrm{out}-R_\mathrm{in}} + \frac{\bar{x}^2}{\bar{\rho}^2\rho}
     182             :  *  \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
     183             :  * \frac{\partial \bar{y}}{\partial x_0^0} &=
     184             :  *  - \frac{\bar{x}\bar{y}}{\bar{\rho}^2\rho}
     185             :  *  \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},
     186             :  * \f}
     187             :  * where \f$\rho\f$ is computed from Eq. (\f$\ref{eq:rho_from_rhobar}\f$).
     188             :  *
     189             :  * ### dxbar_dsigma
     190             :  *
     191             :  * `dxbar_dsigma` returns \f$\partial \bar{x}^i/\partial \sigma\f$,
     192             :  *  where \f$x_0^i\f$ is held fixed.
     193             :  *
     194             :  * From Eq. (\f$\ref{eq:deriv_sigma}\f$) we have
     195             :  *
     196             :  * \f{align}
     197             :  * \frac{\partial \bar{x}^i}{\partial \sigma} &= (0,0,2).
     198             :  * \f}
     199             :  *
     200             :  */
     201           1 : class FlatSide {
     202             :  public:
     203           0 :   FlatSide(const std::array<double, 3>& center, const double inner_radius,
     204             :            const double outer_radius);
     205             : 
     206           0 :   FlatSide() = default;
     207           0 :   ~FlatSide() = default;
     208           0 :   FlatSide(FlatSide&&) = default;
     209           0 :   FlatSide(const FlatSide&) = default;
     210           0 :   FlatSide& operator=(const FlatSide&) = default;
     211           0 :   FlatSide& operator=(FlatSide&&) = default;
     212             : 
     213             :   template <typename T>
     214           0 :   void forward_map(gsl::not_null<std::array<T, 3>*> target_coords,
     215             :                    const std::array<T, 3>& source_coords) const;
     216             : 
     217           0 :   std::optional<std::array<double, 3>> inverse(
     218             :       const std::array<double, 3>& target_coords, double sigma_in) const;
     219             : 
     220             :   template <typename T>
     221           0 :   void jacobian(gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> jacobian_out,
     222             :                 const std::array<T, 3>& source_coords) const;
     223             : 
     224             :   template <typename T>
     225           0 :   void inv_jacobian(
     226             :       gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> inv_jacobian_out,
     227             :       const std::array<T, 3>& source_coords) const;
     228             : 
     229             :   template <typename T>
     230           0 :   void sigma(gsl::not_null<T*> sigma_out,
     231             :              const std::array<T, 3>& source_coords) const;
     232             : 
     233             :   template <typename T>
     234           0 :   void deriv_sigma(gsl::not_null<std::array<T, 3>*> deriv_sigma_out,
     235             :                    const std::array<T, 3>& source_coords) const;
     236             : 
     237             :   template <typename T>
     238           0 :   void dxbar_dsigma(gsl::not_null<std::array<T, 3>*> dxbar_dsigma_out,
     239             :                     const std::array<T, 3>& source_coords) const;
     240             : 
     241           0 :   std::optional<double> lambda_tilde(
     242             :       const std::array<double, 3>& parent_mapped_target_coords,
     243             :       const std::array<double, 3>& projection_point,
     244             :       bool source_is_between_focus_and_target) const;
     245             : 
     246             :   template <typename T>
     247           0 :   void deriv_lambda_tilde(
     248             :       gsl::not_null<std::array<T, 3>*> deriv_lambda_tilde_out,
     249             :       const std::array<T, 3>& target_coords, const T& lambda_tilde,
     250             :       const std::array<double, 3>& projection_point) const;
     251             : 
     252             :   // NOLINTNEXTLINE(google-runtime-references)
     253           0 :   void pup(PUP::er& p);
     254             : 
     255           0 :   static bool is_identity() { return false; }
     256             : 
     257           0 :   static constexpr bool supports_hessian{false};
     258             : 
     259             :  private:
     260           0 :   friend bool operator==(const FlatSide& lhs, const FlatSide& rhs);
     261           0 :   std::array<double, 3> center_{};
     262           0 :   double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
     263           0 :   double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
     264             : };
     265           0 : bool operator!=(const FlatSide& lhs, const FlatSide& rhs);
     266             : }  // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps