fixup. Wedge

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp

@@ -99,19 +99,14 @@ T SphereTransition::call_impl(const std::array<T, 3>& source_coords) const {
 template <typename T>
 std::array<T, 3> SphereTransition::gradient_impl(
     const std::array<T, 3>& source_coords) const {
+  if (UNLIKELY(extend_)) {
+    ERROR(
+        "The gradient of the Sphere transition isn't supported when we are "
+        "extending the transition beyond the inner and outer surfaces.");
+  }
   const T mag = magnitude(source_coords);
   check_magnitudes(mag);
-  std::array<T, 3> result = a_ * source_coords / mag;
-  if (extend_) {
-    for (size_t i = 0; i < get_size(mag); i++) {
-      if (get_element(mag, i) < r_min_ or get_element(mag, i) > r_max_) {
-        for (size_t j = 0; j < 3; j++) {
-          get_element(gsl::at(result, j), i) = 0.0;
-        }
-      }
-    }
-  }
-  return result;
+  return a_ * source_coords / mag;
 }
 
 bool SphereTransition::operator==(

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.hpp

@@ -32,7 +32,8 @@ namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
  *
  * If \p extend is set to `true`, then the function can be called beyond `r_min`
  * and `r_max`. Within `r_min` the value is 1, and outside `r_max` the value is
- * 0. This is reversed if \p reverse is true.
+ * 0. This is reversed if \p reverse is true. However, the gradient function
+ * cannot be called while \p extend is true.
  */
 class SphereTransition final : public ShapeMapTransitionFunction {
  public:

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.cpp

@@ -23,6 +23,7 @@
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/ErrorHandling/CaptureForError.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/Gsl.hpp"
 #include "Utilities/MakeArray.hpp"
 #include "Utilities/MakeString.hpp"
 #include "Utilities/MakeWithValue.hpp"
@@ -62,12 +63,13 @@ Wedge::Wedge(const std::array<double, 3>& inner_center,
              const double inner_radius, const double inner_sphericity,
              const std::array<double, 3>& outer_center,
              const double outer_radius, const double outer_sphericity,
-             const Axis axis, const bool reverse)
+             const Axis axis, const bool reverse, const bool extend)
     : inner_surface_(inner_center, inner_radius, inner_sphericity),
       outer_surface_(outer_center, outer_radius, outer_sphericity),
       projection_center_(inner_surface_.center - outer_surface_.center),
       axis_(axis),
-      reverse_(reverse) {
+      reverse_(reverse),
+      extend_(extend) {
   if (projection_center_ != std::array{0.0, 0.0, 0.0} and
       inner_surface_.sphericity != 1.0) {
     ERROR(
@@ -76,6 +78,14 @@ Wedge::Wedge(const std::array<double, 3>& inner_center,
         << inner_surface_.center
         << "\nouter center: " << outer_surface_.center);
   }
+  if (projection_center_ != std::array{0.0, 0.0, 0.0} and extend_) {
+    ERROR(
+        "You cannot extend the transition outside the inner and outer surfaces "
+        "if centers of the inner and outer object are different.\ninner "
+        "center: "
+        << inner_surface_.center
+        << "\nouter center: " << outer_surface_.center);
+  }
   for (size_t i = 0; i < inner_surface_.center.size(); i++) {
     if (abs(gsl::at(projection_center_, i)) >=
         outer_surface_.half_cube_length) {
@@ -342,8 +352,13 @@ T Wedge::call_impl(const std::array<T, 3>& source_coords) const {
   check_distances(inner_distance, outer_distance, centered_coords_magnitude,
                   source_coords);
 
-  const auto result = (outer_distance - centered_coords_magnitude) /
-                      (outer_distance - inner_distance);
+  T result = (outer_distance - centered_coords_magnitude) /
+             (outer_distance - inner_distance);
+
+  if (extend_) {
+    result = blaze::clamp(result, 0.0, 1.0);
+  }
+
   if (reverse_) {
     return 1.0 - result;
   } else {
@@ -384,7 +399,7 @@ std::optional<double> Wedge::original_radius_over_radius(
   // already checked above. This logic is reversed if we are in reverse mode.
   if ((not reverse_ and (centered_coords_magnitude > outer_distance + eps_)) or
       (reverse_ and (centered_coords_magnitude < inner_distance - eps_))) {
-    return std::nullopt;
+    return extend_ ? std::optional{1.0} : std::nullopt;
   }
 
   // If distorted radius is 0, this means the map is the identity so the radius
@@ -423,10 +438,27 @@ std::optional<double> Wedge::original_radius_over_radius(
   const double original_radius =
       original_radius_over_radius * centered_coords_magnitude;
 
-  return (original_radius + eps_) >= inner_distance and
-                 (original_radius - eps_) <= outer_distance
-             ? std::optional<double>{original_radius_over_radius}
-             : std::nullopt;
+  // If we are within the inner distance and the outer distance, doesn't matter
+  // if we are reverse or extending, we just return the value we calculated. If
+  // we are extending, then if we are reverse and the point is outside the outer
+  // distance or we aren't reversed and the point is inside the inner distance,
+  // then we return a simplified formula since the transition function is 1 in
+  // this region. If we are extending, but the above conditions aren't true,
+  // then we are in the region where the transition function is 0, so we
+  // return 1.0. Otherwise we return nullopt.
+  if ((original_radius + eps_) >= inner_distance and
+      (original_radius - eps_) <= outer_distance) {
+    return std::optional{original_radius_over_radius};
+  } else if (extend_) {
+    if ((not reverse_ and original_radius < inner_distance) or
+        (reverse_ and original_radius > outer_distance)) {
+      return std::optional{1.0 + radial_distortion / centered_coords_magnitude};
+    } else {
+      return std::optional{1.0};
+    }
+  } else {
+    return std::nullopt;
+  }
 }
 
 std::array<double, 3> Wedge::gradient(
@@ -441,6 +473,11 @@ std::array<DataVector, 3> Wedge::gradient(
 template <typename T>
 std::array<T, 3> Wedge::gradient_impl(
     const std::array<T, 3>& source_coords) const {
+  if (UNLIKELY(extend_)) {
+    ERROR(
+        "The gradient of the Wedge transition isn't supported when we are "
+        "extending the transition beyond the inner and outer surfaces.");
+  }
   // The source coords are centered
   const T centered_coords_magnitude = magnitude(source_coords);
   const T one_over_centered_coords_magnitude = 1.0 / centered_coords_magnitude;
@@ -538,7 +575,8 @@ bool Wedge::operator==(const ShapeMapTransitionFunction& other) const {
   const Wedge& other_ref = *dynamic_cast<const Wedge*>(&other);
   return surface_equal(inner_surface_, other_ref.inner_surface_) and
          surface_equal(outer_surface_, other_ref.outer_surface_) and
-         axis_ == other_ref.axis_ and reverse_ == other_ref.reverse_;
+         axis_ == other_ref.axis_ and reverse_ == other_ref.reverse_ and
+         extend_ == other_ref.extend_;
 }
 
 bool Wedge::operator!=(const ShapeMapTransitionFunction& other) const {
@@ -555,8 +593,11 @@ void Wedge::check_distances(
   const T result = (outer_distance - centered_coords_magnitude) /
                    (outer_distance - inner_distance);
   for (size_t i = 0; i < get_size(centered_coords_magnitude); ++i) {
-    if (get_element(result, i) + eps_ < 0.0 or
-        get_element(result, i) - eps_ > 1.0) {
+    const bool point_is_bad =
+        extend_ ? get_element(centered_coords_magnitude, i) <= 0.0
+                : (get_element(result, i) + eps_ < 0.0 or
+                   get_element(result, i) - eps_ > 1.0);
+    if (point_is_bad) {
       ERROR(
           "The Wedge transition map was called with coordinates outside "
           "the set inner and outer surfaces.\nThe requested (centered) point "
@@ -598,7 +639,7 @@ struct LegacySurface {
 
 void Wedge::pup(PUP::er& p) {
   ShapeMapTransitionFunction::pup(p);
-  size_t version = 2;
+  size_t version = 3;
   p | version;
   // Remember to increment the version number when making changes to this
   // function. Retain support for unpacking data written by previous versions
@@ -613,6 +654,11 @@ void Wedge::pup(PUP::er& p) {
     } else {
       reverse_ = false;
     }
+    if (version >= 3) {
+      p | extend_;
+    } else {
+      extend_ = false;
+    }
   } else if (p.isUnpacking()) {
     LegacySurface inner_surface{};
     LegacySurface outer_surface{};

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.hpp

@@ -47,6 +47,14 @@ namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
  * \label{eq:transition_func_reverse}
  * \end{equation}
  *
+ * \parblock
+ *
+ * \note If \p extend is true, then within $D_{\text{in}}(r, \theta, \phi)$,
+ * $f(r, \theta, \phi) = 1$ and beyond $D_{\text{out}}(r, \theta, \phi)$, $f(r,
+ * \theta, \phi) = 0$. If \p reverse is true, this logic is flipped.
+ *
+ * \endparblock
+ *
  * Here, $s$ is the sphericity of the surface which goes from 0 (flat) to 1
  * (spherical), $R$ is the radius of the spherical surface, $\text{out}$ is the
  * outer surface, and $\text{in}$ is the inner surface. If the sphericity is 1,
@@ -242,6 +250,19 @@ namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
  *
  * \endparblock
  *
+ * \parblock
+ *
+ * \note If \p extend is true and we are inside the inner surface, then this
+ * simplifies to
+ * \begin{equation}
+ * \frac{r}{\tilde{r}} = 1 + \frac{\Sigma(\theta, \phi)}{\tilde{r}}
+ * \end{equation}
+ * because $f(r, \theta, \phi) = 1$. If we are outside the outer surface, then
+ * $r/\tilde{r} = 1$ because $f(r, \theta, \phi) = 0$. If \p reverse is true,
+ * this logic is reversed.
+ *
+ * \endparblock
+ *
  * ## Gradient
  *
  * The cartesian gradient of the transition function is
@@ -257,6 +278,12 @@ namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
  *
  * \note If \p reverse is true, the gradient picks up an overall factor of -1.0.
  *
+ * \parblock
+ *
+ * \note The gradient is not supported if \p extend is true.
+ *
+ * \endparblock
+ *
  * Therefore, we need to compute the gradients of $\vec x_0$ and $\vec x_1$.
  *
  * ### Gradient of vector to inner surface
@@ -413,11 +440,15 @@ class Wedge final : public ShapeMapTransitionFunction {
    * boundary and 1 at the outer boundary (useful for deforming star surfaces).
    * Otherwise, it will be 1 at the inner boundary and 0 at the outer boundary
    * (useful for deforming black hole excision surfaces).
+   * \param extend If true, allows the transition function to be called inside
+   * the inner surface (with a value of 1) and beyond the outer surface (with a
+   * value of 0). This logic is reversed if \p reverse is true. The gradient is
+   * not supported when \p extend is true.
    */
   Wedge(const std::array<double, 3>& inner_center, double inner_radius,
         double inner_sphericity, const std::array<double, 3>& outer_center,
         double outer_radius, double outer_sphericity, Axis axis,
-        bool reverse = false);
+        bool reverse = false, bool extend = false);
 
   double operator()(const std::array<double, 3>& source_coords) const override;
   DataVector operator()(
@@ -499,6 +530,7 @@ class Wedge final : public ShapeMapTransitionFunction {
   std::array<double, 3> projection_center_{};
   Axis axis_{};
   bool reverse_{false};
+  bool extend_{false};
 
   static constexpr double eps_ = std::numeric_limits<double>::epsilon() * 100;
 };