Extrapolate from nearest element

src/Domain/BlockLogicalCoordinates.cpp

@@ -21,7 +21,8 @@ template <size_t Dim, typename Fr>
 std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>>
 block_logical_coordinates_single_point(
     const tnsr::I<double, Dim, Fr>& input_point, const Block<Dim>& block,
-    const double time, const domain::FunctionsOfTimeMap& functions_of_time) {
+    const double time, const domain::FunctionsOfTimeMap& functions_of_time,
+    const bool allow_extrapolation) {
   std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>> logical_point{};
   if (block.is_time_dependent()) {
     if constexpr (std::is_same_v<Fr, ::Frame::Inertial>) {
@@ -125,14 +126,57 @@ block_logical_coordinates_single_point(
       logical_point->get(d) = -1.0;
       continue;
     }
-    if (abs(logical_point->get(d)) > 1.0) {
+    if (abs(logical_point->get(d)) > 1.0 and not allow_extrapolation) {
       return std::nullopt;
     }
   }
 
   return logical_point;
 }
 
+template <size_t Dim, typename Fr>
+BlockLogicalCoords<Dim> block_logical_coordinates_in_excision(
+    const tnsr::I<double, Dim, Fr>& input_point,
+    const ExcisionSphere<Dim>& excision_sphere,
+    const std::vector<Block<Dim>>& blocks, const double time,
+    const domain::FunctionsOfTimeMap& functions_of_time) {
+  // Comment by NV (Apr 2024): This function can be made more robust by first
+  // checking if the point is inside the excision sphere at all, but the
+  // excision sphere doesn't currently have this information. It needs the
+  // grid-to-inertial map including the deformation of the excision sphere to
+  // determine this.
+  for (const auto& [block_id, direction] :
+       excision_sphere.abutting_directions()) {
+    auto x_logical = block_logical_coordinates_single_point(
+        input_point, blocks[block_id], time, functions_of_time, true);
+    if (not x_logical.has_value()) {
+      continue;
+    }
+    // Discard block if the point has angular logical coordinates outside the
+    // range [-1, 1]
+    for (size_t d = 0; d < Dim; ++d) {
+      if (d != direction.dimension() and std::abs(x_logical->get(d)) > 1.) {
+        x_logical = std::nullopt;
+        break;
+      }
+    }
+    if (not x_logical.has_value()) {
+      continue;
+    }
+    // Discard block if the point is radially inside the block or on the other
+    // side of the excision sphere
+    const double radial_distance_this_block =
+        x_logical->get(direction.dimension()) * direction.sign();
+    if (radial_distance_this_block < 1.) {
+      continue;
+    }
+    // The checks above should leave only 1 valid block, so return that
+    return make_id_pair(domain::BlockId(block_id),
+                        std::move(x_logical.value()));
+  }
+  return std::nullopt;
+}
+
 template <size_t Dim, typename Fr>
 std::vector<BlockLogicalCoords<Dim>> block_logical_coordinates(
     const Domain<Dim>& domain, const tnsr::I<DataVector, Dim, Fr>& x,
@@ -174,11 +218,18 @@ std::vector<BlockLogicalCoords<Dim>> block_logical_coordinates(
   block_logical_coordinates_single_point(                                      \
       const tnsr::I<double, DIM(data), FRAME(data)>& input_point,              \
       const Block<DIM(data)>& block, const double time,                        \
-      const domain::FunctionsOfTimeMap& functions_of_time);                    \
+      const domain::FunctionsOfTimeMap& functions_of_time,                     \
+      bool allow_extrapolation);                                               \
   template std::vector<BlockLogicalCoords<DIM(data)>>                          \
   block_logical_coordinates(                                                   \
       const Domain<DIM(data)>& domain,                                         \
       const tnsr::I<DataVector, DIM(data), FRAME(data)>& x, const double time, \
+      const domain::FunctionsOfTimeMap& functions_of_time);                    \
+  template BlockLogicalCoords<DIM(data)>                                       \
+  block_logical_coordinates_in_excision(                                       \
+      const tnsr::I<double, DIM(data), FRAME(data)>& input_point,              \
+      const ExcisionSphere<DIM(data)>& excision_sphere,                        \
+      const std::vector<Block<DIM(data)>>& blocks, const double time,          \
       const domain::FunctionsOfTimeMap& functions_of_time);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3),

src/Domain/BlockLogicalCoordinates.hpp

@@ -21,6 +21,8 @@ template <size_t VolumeDim>
 class Domain;
 template <size_t VolumeDim>
 class Block;
+template <size_t VolumeDim>
+class ExcisionSphere;
 /// \endcond
 
 template <size_t Dim>
@@ -69,10 +71,32 @@ auto block_logical_coordinates(
     const domain::FunctionsOfTimeMap& functions_of_time = {})
     -> std::vector<BlockLogicalCoords<Dim>>;
 
+/// \par Extrapolation
+/// If `allow_extrapolation` is `true` (default is `false`), then the function
+/// can return block-logical coordinates outside of [-1, 1]. This is useful for
+/// extrapolating into excision regions.
 template <size_t Dim, typename Fr>
 std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>>
 block_logical_coordinates_single_point(
     const tnsr::I<double, Dim, Fr>& input_point, const Block<Dim>& block,
     double time = std::numeric_limits<double>::signaling_NaN(),
-    const domain::FunctionsOfTimeMap& functions_of_time = {});
+    const domain::FunctionsOfTimeMap& functions_of_time = {},
+    bool allow_extrapolation = false);
 /// @}
+
+/*!
+ * \brief Finds the block-logical coordinates of a point in an excision sphere.
+ *
+ * The point will be reported to belong to the nearest block that abuts the
+ * excision sphere and will have a logical coordinate outside the range [-1, 1]
+ * in the radial direction. This is useful for extrapolating into excision
+ * regions using the Lagrange polynomial basis of the nearest element. If the
+ * point is not within the excision sphere, the return value is std::nullopt.
+ */
+template <size_t Dim, typename Frame>
+BlockLogicalCoords<Dim> block_logical_coordinates_in_excision(
+    const tnsr::I<double, Dim, Frame>& input_point,
+    const ExcisionSphere<Dim>& excision_sphere,
+    const std::vector<Block<Dim>>& blocks,
+    double time = std::numeric_limits<double>::signaling_NaN(),
+    const domain::FunctionsOfTimeMap& functions_of_time = {});

src/Domain/ElementLogicalCoordinates.cpp

@@ -27,36 +27,29 @@ element_logical_coordinates(
     const ElementId<Dim>& element_id) {
   tnsr::I<double, Dim, Frame::ElementLogical> x_element_logical{};
   for (size_t d = 0; d < Dim; ++d) {
-    // Check if the point is outside the element
-    const double up = element_id.segment_id(d).endpoint(Side::Upper);
-    const double lo = element_id.segment_id(d).endpoint(Side::Lower);
-    if (x_block_logical.get(d) < lo or x_block_logical.get(d) > up) {
+    // Check if the point is inside the element in this dimension.
+    // The following conditions are valid:
+    // - The block-logical coordinates are within the element's bounds.
+    // - The block-logical coordinate is outside the block and we're in the
+    //   outermost element in that direction. This means we want to extrapolate
+    //   out of the block.
+    const auto& segment_id = element_id.segment_id(d);
+    const double up = segment_id.endpoint(Side::Upper);
+    const double lo = segment_id.endpoint(Side::Lower);
+    if ((x_block_logical.get(d) >= lo and x_block_logical.get(d) <= up) or
+        (x_block_logical.get(d) < -1. and segment_id.index() == 0) or
+        (x_block_logical.get(d) > 1. and
+         segment_id.index() == two_to_the(segment_id.refinement_level()) - 1)) {
+      // Map to element logical coords
+      x_element_logical.get(d) =
+          (2.0 * x_block_logical.get(d) - up - lo) / (up - lo);
+    } else {
       return std::nullopt;
     }
-    // Map to element logical coords
-    x_element_logical.get(d) =
-        (2.0 * x_block_logical.get(d) - up - lo) / (up - lo);
   }
   return x_element_logical;
 }
 
-namespace {
-// The segments bounds are binary fractions (i.e. the numerator is an
-// integer and the denominator is a power of 2) so these floating point
-// comparisons should be safe from roundoff problems
-// Need to return true if on upper face of block
-// Otherwise return true if point is within the segment or on the lower bound
-bool segment_contains(const double x_block_logical,
-                      const double lower_bound_block_logical,
-                      const double upper_bound_block_logical) {
-  if (UNLIKELY(x_block_logical == upper_bound_block_logical)) {
-    return (upper_bound_block_logical == 1.0);
-  }
-  return (x_block_logical >= lower_bound_block_logical and
-          x_block_logical < upper_bound_block_logical);
-}
-}  // namespace
-
 template <size_t Dim>
 std::unordered_map<ElementId<Dim>, ElementLogicalCoordHolder<Dim>>
 element_logical_coordinates(
@@ -90,13 +83,15 @@ element_logical_coordinates(
         if (not x_elem.has_value()) {
           continue;
         }
-        // Disambiguate points on shared element boundaries
+        // Disambiguate points on shared element boundaries. To do this,
+        // associate boundary points with the element with the lower index.
         bool is_contained = true;
         for (size_t d = 0; d < Dim; ++d) {
-          const double up = element_id.segment_id(d).endpoint(Side::Upper);
-          const double lo = element_id.segment_id(d).endpoint(Side::Lower);
-          const double x_block_log = x_block_logical.get(d);
-          if (not segment_contains(x_block_log, lo, up)) {
+          const auto& segment_id = element_id.segment_id(d);
+          const double up = segment_id.endpoint(Side::Upper);
+          if (x_block_logical.get(d) == up and
+              segment_id.index() !=
+                  two_to_the(segment_id.refinement_level()) - 1) {
             is_contained = false;
             break;
           }

src/Domain/ElementLogicalCoordinates.hpp

@@ -32,6 +32,11 @@ class IdPair;
  * have element logical coordinates of -1 or 1. See the other function overload
  * for handling multiple points and disambiguating points on shared element
  * boundaries.
+ *
+ * Points with block-logical coordinates outside [-1, 1] are considered to
+ * belong to the outermost element of the block in that direction, and will
+ * therefore return an element logical coordinate outside of [-1, 1] as well.
+ * This is useful for extrapolation.
  */
 template <size_t Dim>
 std::optional<tnsr::I<double, Dim, Frame::ElementLogical>>

tests/Unit/Domain/Test_BlockAndElementLogicalCoordinates.cpp

@@ -33,6 +33,7 @@
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Structure/InitialElementIds.hpp"
 #include "Framework/TestHelpers.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeArray.hpp"
 #include "Utilities/StdHelpers.hpp"
@@ -817,6 +818,86 @@ void test_block_logical_coordinates_with_roundoff_error() {
   CHECK(get<1>(block_logical_coords[5]->data) < 1.0);
   CHECK(get<0>(block_logical_coords[6]->data) < 1.0);
 }
+
+void test_excision_spheres(const bool time_dependent) {
+  CAPTURE(time_dependent);
+  const auto shell = domain::creators::Sphere(
+      1., 3., domain::creators::Sphere::Excision{}, 0_st, 3_st, true,
+      std::nullopt, {}, domain::CoordinateMaps::Distribution::Linear,
+      ShellWedges::All,
+      time_dependent
+          ? std::make_optional(domain::creators::Sphere::TimeDepOptionType{
+                domain::creators::sphere::TimeDependentMapOptions{
+                    0.0,
+                    {{12,
+                      domain::creators::time_dependent_options::
+                          KerrSchildFromBoyerLindquist{1.0, {{0.0, 0.0, 0.9}}},
+                      std::nullopt}},
+                    std::nullopt,
+                    std::nullopt,
+                    std::nullopt,
+                    false}})
+          : std::nullopt);
+  const auto domain = shell.create_domain();
+  const auto& excision_sphere = domain.excision_spheres().at("ExcisionSphere");
+  const double time = 0.0;
+  const auto functions_of_time = shell.functions_of_time();
+
+  std::vector<std::pair<tnsr::I<double, 3>, BlockLogicalCoords<3>>> test_data{};
+  // Outside the excision sphere
+  test_data.push_back({tnsr::I<double, 3>{{2., 0., 0.}}, std::nullopt});
+  if (time_dependent) {
+    // On the boundary of the excision sphere
+    test_data.push_back(
+        {tnsr::I<double, 3>{
+             {get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist<
+                  double>(1.0, {{M_PI_2, 0.0}}, 1.0, {{0.0, 0.0, 0.9}})),
+              0., 0.}},
+         {{domain::BlockId{4},
+           tnsr::I<double, 3, Frame::BlockLogical>{{0., 0., -1.}}}}});
+    // Inside the excision sphere
+    test_data.push_back(
+        {tnsr::I<double, 3>{
+             {get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist<
+                  double>(0.5, {{M_PI_2, 0.0}}, 1.0, {{0.0, 0.0, 0.9}})),
+              0., 0.}},
+         {{domain::BlockId{4},
+           // If the shape map preserved the BL-KS transformation within the
+           // excision, the result should be -1.5. However, apparently it
+           // doesn't.
+           tnsr::I<double, 3, Frame::BlockLogical>{
+               {0., 0., -1.31579939060867179}}}}});
+  } else {
+    // On the boundary of the excision sphere
+    test_data.push_back(
+        {tnsr::I<double, 3>{{1., 0., 0.}},
+         {{domain::BlockId{4},
+           tnsr::I<double, 3, Frame::BlockLogical>{{0., 0., -1.}}}}});
+    // Inside the excision sphere
+    test_data.push_back(
+        {tnsr::I<double, 3>{{0.5, 0., 0.}},
+         {{domain::BlockId{4},
+           tnsr::I<double, 3, Frame::BlockLogical>{{0., 0., -1.5}}}}});
+  }
+
+  for (const auto& [input_point, expected] : test_data) {
+    CAPTURE(input_point);
+    const auto result = block_logical_coordinates_in_excision(
+        input_point, excision_sphere, domain.blocks(), time, functions_of_time);
+    CHECK(result.has_value() == expected.has_value());
+    if (result.has_value()) {
+      CHECK(result->id.get_index() == expected->id.get_index());
+      CHECK_ITERABLE_APPROX(result->data, expected->data);
+      const auto element_logical_coords = element_logical_coordinates(
+          result->data, ElementId<3>{result->id.get_index()});
+      REQUIRE(element_logical_coords.has_value());
+      const tnsr::I<double, 3, Frame::ElementLogical> expected_element_logical{
+          {get<0>(expected->data), get<1>(expected->data),
+           get<2>(expected->data)}};
+      CHECK_ITERABLE_APPROX(*element_logical_coords, expected_element_logical);
+    }
+  }
+}
 }  // namespace
 
 SPECTRE_TEST_CASE("Unit.Domain.BlockAndElementLogicalCoords",
@@ -834,4 +915,6 @@ SPECTRE_TEST_CASE("Unit.Domain.BlockAndElementLogicalCoords",
   test_block_logical_coordinates1fail();
   test_element_ids_are_uniquely_determined();
   test_block_logical_coordinates_with_roundoff_error();
+  test_excision_spheres(false);
+  test_excision_spheres(true);
 }