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             : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
      14             : 
      15             : /// \cond
      16             : namespace PUP {
      17             : class er;
      18             : }  // namespace PUP
      19             : /// \endcond
      20             : 
      21             : /// Contains FocallyLiftedInnerMaps
      22             : namespace domain::CoordinateMaps::FocallyLiftedInnerMaps {
      23             : /*!
      24             :  * \brief A FocallyLiftedInnerMap that maps a 3D unit right cylinder
      25             :  *  to a volume that connects a portion of a plane and a spherical
      26             :  *  surface.
      27             :  *
      28             :  * \details The domain of the map is a 3D unit right cylinder with
      29             :  * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
      30             :  * \f$-1\leq\bar{z}\leq 1\f$ and \f$\bar{x}^2+\bar{y}^2 \leq
      31             :  * 1\f$.  The range of the map has coordinates \f$(x,y,z)\f$.
      32             :  *
      33             :  * Consider a 2D circle in 3D space that is normal to the \f$z\f$ axis
      34             :  * and has (3D) center \f$C^i\f$ and radius \f$R\f$.  `FlatEndcap`
      35             :  * provides the following functions:
      36             :  *
      37             :  * ### forward_map()
      38             :  * `forward_map()` maps \f$(\bar{x},\bar{y},\bar{z}=-1)\f$ to the interior
      39             :  * of the circle.  The arguments to `forward_map()`
      40             :  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
      41             :  * `forward_map()` returns \f$x_0^i\f$,
      42             :  * the 3D coordinates on the circle, which are given by
      43             :  *
      44             :  * \f{align}
      45             :  * x_0^0 &= R \bar{x} + C^0,\\
      46             :  * x_0^1 &= R \bar{y} + C^1,\\
      47             :  * x_0^2 &= C^2.
      48             :  * \f}
      49             :  *
      50             :  * ### sigma
      51             :  *
      52             :  * \f$\sigma\f$ is a function that is zero on the plane
      53             :  * \f$x^i=x_0^i\f$ and unity at \f$\bar{z}=+1\f$ (corresponding to the
      54             :  * upper surface of the FocallyLiftedMap). We define
      55             :  *
      56             :  * \f{align}
      57             :  *  \sigma &= \frac{\bar{z}+1}{2}.
      58             :  * \f}
      59             :  *
      60             :  * ### deriv_sigma
      61             :  *
      62             :  * `deriv_sigma` returns
      63             :  *
      64             :  * \f{align}
      65             :  *  \frac{\partial \sigma}{\partial \bar{x}^j} &= (0,0,1/2).
      66             :  * \f}
      67             :  *
      68             :  * ### jacobian
      69             :  *
      70             :  * `jacobian` returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
      71             :  * The arguments to `jacobian`
      72             :  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
      73             :  *
      74             :  * Differentiating Eqs.(1--3) above yields
      75             :  *
      76             :  * \f{align*}
      77             :  * \frac{\partial x_0^0}{\partial \bar{x}} &= R,\\
      78             :  * \frac{\partial x_0^1}{\partial \bar{y}} &= R,
      79             :  * \f}
      80             :  * and all other components are zero.
      81             :  *
      82             :  * ### inverse
      83             :  *
      84             :  * `inverse` takes \f$x_0^i\f$ and \f$\sigma\f$ as arguments, and
      85             :  * returns \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed
      86             :  * `std::optional<std::array<double, 3>>` if
      87             :  * \f$x_0^i\f$ or \f$\sigma\f$ are outside the range of the map.
      88             :  * The formula for the inverse is straightforward:
      89             :  *
      90             :  * \f{align}
      91             :  *  \bar{x} &= \frac{x_0^0-C^0}{R},\\
      92             :  *  \bar{y} &= \frac{x_0^1-C^1}{R},\\
      93             :  *  \bar{z} &= 2\sigma - 1.
      94             :  * \f}
      95             :  *
      96             :  * If \f$\bar{z}\f$ is outside the range \f$[-1,1]\f$ or
      97             :  * if \f$\bar{x}^2+\bar{y}^2 > 1\f$ then we return
      98             :  * a default-constructed `std::optional<std::array<double, 3>>`
      99             :  *
     100             :  * ### lambda_tilde
     101             :  *
     102             :  * `lambda_tilde` takes as arguments a point \f$x^i\f$ and a projection point
     103             :  *  \f$P^i\f$, and computes \f$\tilde{\lambda}\f$, the solution to
     104             :  *
     105             :  * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
     106             :  *
     107             :  * Since \f$x_0^i\f$ must lie on the plane \f$x_0^3=C^3\f$,
     108             :  *
     109             :  * \f{align} \tilde{\lambda} &= \frac{C^3-P^3}{x^3-P^3}.\f}
     110             :  *
     111             :  * The valid range of \f$\tilde{\lambda}\f$ depends on whether the source
     112             :  * lies between the focus and the target:
     113             :  * - Non-interior case (`source_is_between_focus_and_target`=`false`):
     114             :  * \f$x_0^i\f$ lies beyond \f$x^i\f$ from \f$P^i\f$, so \f$\tilde{\lambda}\ge
     115             :  * 1\f$.  A default-constructed `std::optional<double>` is returned if
     116             :  * \f$\tilde{\lambda}<1\f$.
     117             :  * - Interior case (`source_is_between_focus_and_target`=`true`): \f$x_0^i\f$
     118             :  * lies between \f$P^i\f$ and \f$x^i\f$, so \f$\tilde{\lambda}\in(0,1]\f$.  A
     119             :  * default-constructed `std::optional<double>` is returned if
     120             :  * \f$\tilde{\lambda}\le 0\f$ or \f$\tilde{\lambda}>1\f$. In both cases a
     121             :  * default-constructed `std::optional<double>` is also returned if \f$x^3 =
     122             :  * P^3\f$ (the ray from \f$P\f$ is parallel to the disk and never intersects
     123             :  * it).
     124             :  *
     125             :  * ### deriv_lambda_tilde
     126             :  *
     127             :  * `deriv_lambda_tilde` takes as arguments \f$x_0^i\f$, a projection point
     128             :  *  \f$P^i\f$, and \f$\tilde{\lambda}\f$, and
     129             :  *  returns \f$\partial \tilde{\lambda}/\partial x^i\f$. We have
     130             :  *
     131             :  * \f{align}
     132             :  * \frac{\partial\tilde{\lambda}}{\partial x^3} =
     133             :  * -\frac{C^3-P^3}{(x^3-P^3)^2} = -\frac{\tilde{\lambda}^2}{C^3-P^3},
     134             :  * \f}
     135             :  * and other components are zero.
     136             :  *
     137             :  * ### inv_jacobian
     138             :  *
     139             :  * `inv_jacobian` returns \f$\partial \bar{x}^i/\partial x_0^k\f$,
     140             :  *  where \f$\sigma\f$ is held fixed.
     141             :  *  The arguments to `inv_jacobian`
     142             :  *  are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
     143             :  *
     144             :  * The nonzero components are
     145             :  * \f{align}
     146             :  * \frac{\partial \bar{x}}{\partial x_0^0} &= \frac{1}{R},\\
     147             :  * \frac{\partial \bar{y}}{\partial x_0^1} &= \frac{1}{R}.
     148             :  * \f}
     149             :  *
     150             :  * ### dxbar_dsigma
     151             :  *
     152             :  * `dxbar_dsigma` returns \f$\partial \bar{x}^i/\partial \sigma\f$,
     153             :  *  where \f$x_0^i\f$ is held fixed.
     154             :  *
     155             :  * From Eq. (6) we have
     156             :  *
     157             :  * \f{align}
     158             :  * \frac{\partial \bar{x}^i}{\partial \sigma} &= (0,0,2).
     159             :  * \f}
     160             :  *
     161             :  */
     162           1 : class FlatEndcap {
     163             :  public:
     164           0 :   FlatEndcap(const std::array<double, 3>& center, double radius);
     165             : 
     166           0 :   FlatEndcap() = default;
     167           0 :   ~FlatEndcap() = default;
     168           0 :   FlatEndcap(FlatEndcap&&) = default;
     169           0 :   FlatEndcap(const FlatEndcap&) = default;
     170           0 :   FlatEndcap& operator=(const FlatEndcap&) = default;
     171           0 :   FlatEndcap& operator=(FlatEndcap&&) = default;
     172             : 
     173             :   template <typename T>
     174           0 :   void forward_map(
     175             :       const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
     176             :           target_coords,
     177             :       const std::array<T, 3>& source_coords) const;
     178             : 
     179           0 :   std::optional<std::array<double, 3>> inverse(
     180             :       const std::array<double, 3>& target_coords, double sigma_in) const;
     181             : 
     182             :   template <typename T>
     183           0 :   void jacobian(const gsl::not_null<
     184             :                     tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
     185             :                     jacobian_out,
     186             :                 const std::array<T, 3>& source_coords) const;
     187             : 
     188             :   template <typename T>
     189           0 :   void inv_jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3,
     190             :                                                  Frame::NoFrame>*>
     191             :                         inv_jacobian_out,
     192             :                     const std::array<T, 3>& source_coords) const;
     193             : 
     194             :   template <typename T>
     195           0 :   void sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,
     196             :              const std::array<T, 3>& source_coords) const;
     197             : 
     198             :   template <typename T>
     199           0 :   void deriv_sigma(
     200             :       const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
     201             :           deriv_sigma_out,
     202             :       const std::array<T, 3>& source_coords) const;
     203             : 
     204             :   template <typename T>
     205           0 :   void dxbar_dsigma(
     206             :       const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
     207             :           dxbar_dsigma_out,
     208             :       const std::array<T, 3>& source_coords) const;
     209             : 
     210           0 :   std::optional<double> lambda_tilde(
     211             :       const std::array<double, 3>& parent_mapped_target_coords,
     212             :       const std::array<double, 3>& projection_point,
     213             :       bool source_is_between_focus_and_target) const;
     214             : 
     215             :   template <typename T>
     216           0 :   void deriv_lambda_tilde(
     217             :       const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
     218             :           deriv_lambda_tilde_out,
     219             :       const std::array<T, 3>& target_coords, const T& lambda_tilde,
     220             :       const std::array<double, 3>& projection_point) const;
     221             : 
     222             :   // NOLINTNEXTLINE(google-runtime-references)
     223           0 :   void pup(PUP::er& p);
     224             : 
     225           0 :   static bool is_identity() { return false; }
     226             : 
     227           0 :   static constexpr bool supports_hessian{false};
     228             : 
     229             :  private:
     230           0 :   friend bool operator==(const FlatEndcap& lhs, const FlatEndcap& rhs);
     231           0 :   std::array<double, 3> center_{};
     232           0 :   double radius_{std::numeric_limits<double>::signaling_NaN()};
     233             : };
     234           0 : bool operator!=(const FlatEndcap& lhs, const FlatEndcap& rhs);
     235             : }  // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps