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