fixup

src/Domain/BlockLogicalCoordinates.cpp

@@ -153,22 +153,26 @@ BlockLogicalCoords<Dim> block_logical_coordinates_in_excision(
       continue;
     }
     // Discard block if the point has angular logical coordinates outside the
-    // range [-1, 1]
+    // range [-1, 1]. Also discard if the point is on the upper angular face of
+    // the block, to disambiguate which of the two blocks sharing the face the
+    // point is in.
     for (size_t d = 0; d < Dim; ++d) {
-      if (d != direction.dimension() and std::abs(x_logical->get(d)) > 1.) {
+      if (d != direction.dimension() and
+          (x_logical->get(d) >= 1. or 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
+    // Discard block if the point is radially inside the block, or outside the
+    // block but on the side that doesn't abut the excision sphere. Return early
+    // here because this means the point isn't in the excision sphere at all.
     const double radial_distance_this_block =
         x_logical->get(direction.dimension()) * direction.sign();
     if (radial_distance_this_block < 1.) {
-      continue;
+      return std::nullopt;
     }
     // The checks above should leave only 1 valid block, so return that
     return make_id_pair(domain::BlockId(block_id),

src/Domain/BlockLogicalCoordinates.hpp

@@ -91,6 +91,8 @@ block_logical_coordinates_single_point(
  * excision sphere and will have a logical coordinate outside the range [-1, 1]
  * in the radial direction. This is useful for extrapolating data into excision
  * regions.
+ *
+ * Returns `std::nullopt` if the point is not in the given `excision_sphere`.
  */
 template <size_t Dim, typename Frame>
 BlockLogicalCoords<Dim> block_logical_coordinates_in_excision(

src/IO/Exporter/Exporter.cpp

@@ -188,8 +188,8 @@ void add_extrapolation_anchors(
     for (size_t s = 0; s < block_logical_coords->size(); ++s) {
       auto& block_logical_coord = (*block_logical_coords)[s];
       // Process only points inside the excision sphere. These points were found
-      // by `block_logical_coordinates` above, have a radial logical coordinate
-      // outside [-1, 1], and abut the excision sphere.
+      // by `block_logical_coordinates_in_excision`, have a radial logical
+      // coordinate outside [-1, 1], and abut the excision sphere.
       if (not block_logical_coord.has_value()) {
         continue;
       }

src/IO/Exporter/Exporter.hpp

@@ -49,11 +49,11 @@ struct ObservationStep {
  * of the domain (default is `false`). This can be useful to fill the excision
  * region with (constraint-violating but smooth) data so it can be imported into
  * moving puncture codes. Specifically, we implement the strategy used in
- * \cite Etienne2008re adjusted for distorted excisions: we choose uniformly
+ * \cite Etienne2008re adjusted for distorted excisions: we choose 8 uniformly
  * spaced radial anchor points spaced as $\Delta r = 0.3 r_\mathrm{AH}$ in the
  * grid frame (where the excision is spherical), then map the anchor points to
- * the distorted frame (where we have the target point) and do a 7th order
- * polynomial extrapolation into the excision region.
+ * the distorted frame (where we have the target point) and do a polynomial
+ * extrapolation into the excision region.
  * \param num_threads The number of threads to use if OpenMP is linked in. If
  * not specified, OpenMP will determine the number of threads automatically.
  * It's also possible to set the number of threads using the environment