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 <optional>
       8              : #include <pup.h>
       9              : 
      10              : #include "DataStructures/DataVector.hpp"
      11              : #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
      12              : 
      13              : namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
      14              : 
      15              : /*!
      16              :  * \ingroup CoordMapsTimeDependentGroup
      17              :  * \brief A transition function that falls off as $G(r,\theta,\phi) =
      18              :  * \frac{f(r)}{r} = \frac{ar + b}{r}$.
      19              :  *
      20              :  * \details The coefficients $a$ and $b$ are chosen so that the function $f(r) =
      21              :  * ar + b$ falls off linearly from 1 at \p r_min to 0 at \p r_max. The
      22              :  * coefficients are
      23              :  *
      24              :  * \f{align}{
      25              :  * \label{eq:transition_func}
      26              :  * a &= \frac{-1}{r_{\text{max}} - r_{\text{min}}} \\
      27              :  * b &= \frac{r_{\text{max}}}{r_{\text{max}} - r_{\text{min}}} = -a
      28              :  * r_{\text{max}}
      29              :  * \f}
      30              :  *
      31              :  * If \p reverse is set to `true`, then the function falls off from 0 at
      32              :  * \p r_min to 1 at \p r_max. To do this, the coefficients are modified as
      33              :  * $a \rightarrow -a$ and $b \rightarrow 1-b$.
      34              :  *
      35              :  * The function can be called within \p r_min, but only if \p interior is `true`
      36              :  * and \p reverse is `false`. Within \p r_min,
      37              :  *
      38              :  * \begin{equation}
      39              :  * G(r,\theta,\phi) = \frac{r^2}{r_{\text{min}}^3}.
      40              :  * \end{equation}
      41              :  *
      42              :  * which is chosen to match Eq. $\ref{eq:transition_func}$ at $r_{\text{min}}$
      43              :  * and go to 0 at $r=0$.
      44              :  *
      45              :  * This function cannot be called beyond \p r_max.
      46              :  *
      47              :  * If \p interior is `false` and if the `operator()` or `gradient()` is called
      48              :  * with a point within \p r_min, an error will occur and
      49              :  * `original_radius_over_radius` will return `std::nullopt`.
      50              :  *
      51              :  * This is a special, simplified, case of the
      52              :  * `domain::CoordinateMaps::ShapeMapTransitionFunctions::Wedge` class where both
      53              :  * the inner and outer surface are spheres centered on the same center.
      54              :  */
      55            1 : class SphereTransition final : public ShapeMapTransitionFunction {
      56              :  public:
      57            0 :   explicit SphereTransition() = default;
      58            0 :   SphereTransition(double r_min, double r_max, bool reverse = false,
      59              :                    bool interior = false);
      60              : 
      61            1 :   double operator()(
      62              :       const std::array<double, 3>& source_coords,
      63              :       const std::optional<size_t>& one_over_radius_power) const override;
      64            1 :   DataVector operator()(
      65              :       const std::array<DataVector, 3>& source_coords,
      66              :       const std::optional<size_t>& one_over_radius_power) const override;
      67              : 
      68            1 :   std::optional<double> original_radius_over_radius(
      69              :       const std::array<double, 3>& target_coords,
      70              :       double radial_distortion) const override;
      71              : 
      72            1 :   std::array<double, 3> gradient(
      73              :       const std::array<double, 3>& source_coords) const override;
      74            1 :   std::array<DataVector, 3> gradient(
      75              :       const std::array<DataVector, 3>& source_coords) const override;
      76              : 
      77            0 :   WRAPPED_PUPable_decl_template(SphereTransition);
      78            0 :   explicit SphereTransition(CkMigrateMessage* msg);
      79            0 :   void pup(PUP::er& p) override;
      80              : 
      81            0 :   std::unique_ptr<ShapeMapTransitionFunction> get_clone() const override {
      82              :     return std::make_unique<SphereTransition>(*this);
      83              :   }
      84              : 
      85            0 :   bool operator==(const ShapeMapTransitionFunction& other) const override;
      86            0 :   bool operator!=(const ShapeMapTransitionFunction& other) const override;
      87              : 
      88              :  private:
      89              :   template <typename T>
      90            0 :   T call_impl(const std::array<T, 3>& source_coords,
      91              :               const std::optional<size_t>& one_over_radius_power) const;
      92              : 
      93              :   template <typename T>
      94            0 :   std::array<T, 3> gradient_impl(const std::array<T, 3>& source_coords) const;
      95              : 
      96            0 :   double r_min_{};
      97            0 :   double inverse_cube_r_min_{};
      98            0 :   double r_max_{};
      99            0 :   double a_{};
     100            0 :   double b_{};
     101            0 :   bool interior_{};
     102            0 :   static constexpr double eps_ = std::numeric_limits<double>::epsilon() * 100;
     103              : };
     104              : }  // namespace domain::CoordinateMaps::ShapeMapTransitionFunctions