Exporter: allow extrapolation into excisions

This can be useful to fill the excision region with
(non-constraint-satisfying) data so it can be imported into
moving puncture codes.

docs/References.bib

@@ -843,6 +843,21 @@ @article{Dumbser2017okk
     year = "2018"
 }
 
+@article{Etienne2008re,
+  author        = {Etienne, Zachariah B. and Liu, Yuk Tung and Shapiro, Stuart
+                   L. and Baumgarte, Thomas W.},
+  title         = {General relativistic simulations of black-hole-neutron-star
+                   mergers: Effects of black-hole spin},
+  eprint        = {0812.2245},
+  archiveprefix = {arXiv},
+  primaryclass  = {astro-ph},
+  doi           = {10.1103/PhysRevD.79.044024},
+  journal       = {Phys. Rev. D},
+  volume        = {79},
+  pages         = {044024},
+  year          = {2009}
+}
+
 @article{Etienne2010ui,
   author        = "Etienne, Zachariah B. and Liu, Yuk Tung and Shapiro, Stuart
                   L.",

src/Domain/BlockLogicalCoordinates.cpp

@@ -21,7 +21,8 @@ template <size_t Dim, typename Frame>
 std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>>
 block_logical_coordinates_single_point(
     const tnsr::I<double, Dim, Frame>& 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<Frame, ::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 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, 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 Frame>
 std::vector<BlockLogicalCoords<Dim>> block_logical_coordinates(
     const Domain<Dim>& domain, const tnsr::I<DataVector, Dim, Frame>& 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,30 @@ 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 data into excision regions.
 template <size_t Dim, typename Frame>
 std::optional<tnsr::I<double, Dim, ::Frame::BlockLogical>>
 block_logical_coordinates_single_point(
     const tnsr::I<double, Dim, Frame>& 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 one of the blocks that abut the
+ * 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.
+ */
+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,
+    const domain::FunctionsOfTimeMap& functions_of_time);