sxs-collaboration · wthrowe · May 12, 2025 · Apr 14, 2025 · Apr 8, 2025 · May 7, 2025
diff --git a/src/Domain/ElementLogicalCoordinates.cpp b/src/Domain/ElementLogicalCoordinates.cpp
@@ -4,7 +4,9 @@
 #include "Domain/ElementLogicalCoordinates.hpp"
 
 #include <array>
+#include <map>
 #include <memory>
+#include <optional>
 #include <unordered_map>
 #include <utility>
 #include <vector>
@@ -45,6 +47,30 @@ element_logical_coordinates(
   return x_element_logical;
 }
 
+template <size_t Dim>
+std::optional<
+    std::pair<ElementId<Dim>, tnsr::I<double, Dim, Frame::ElementLogical>>>
+element_logical_coordinates(
+    const IdPair<domain::BlockId, tnsr::I<double, Dim, Frame::BlockLogical>>&
+        block_logical_coords,
+    const std::map<size_t, domain::ElementSearchTree<Dim>>& search_tree) {
+  const auto block_search_tree_it =
+      search_tree.find(block_logical_coords.id.get_index());
+  if (block_search_tree_it == search_tree.end()) {
+    return std::nullopt;
+  }
+  const auto& block_search_tree = block_search_tree_it->second;
+  const auto found_element_id = block_search_tree.qbegin(
+      boost::geometry::index::covers(block_logical_coords.data));
+  if (found_element_id == block_search_tree.qend()) {
+    return std::nullopt;
+  }
+  const auto& element_id = *found_element_id;
+  return std::make_pair(element_id, element_logical_coordinates(
+                                        block_logical_coords.data, element_id)
+                                        .value());
+}
+
 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
@@ -152,6 +178,15 @@ element_logical_coordinates(
   element_logical_coordinates(                                                \
       const tnsr::I<double, DIM(data), Frame::BlockLogical>& x_block_logical, \
       const ElementId<DIM(data)>& element_id);                                \
+  template std::optional<                                                     \
+      std::pair<ElementId<DIM(data)>,                                         \
+                tnsr::I<double, DIM(data), Frame::ElementLogical>>>           \
+  element_logical_coordinates(                                                \
+      const IdPair<domain::BlockId,                                           \
+                   tnsr::I<double, DIM(data), Frame::BlockLogical>>&          \
+          block_logical_coords,                                               \
+      const std::map<size_t, domain::ElementSearchTree<DIM(data)>>&           \
+          search_tree);                                                       \
   template std::unordered_map<ElementId<DIM(data)>,                           \
                               ElementLogicalCoordHolder<DIM(data)>>           \
   element_logical_coordinates(                                                \

diff --git a/src/Domain/ElementLogicalCoordinates.hpp b/src/Domain/ElementLogicalCoordinates.hpp
@@ -4,12 +4,14 @@
 #pragma once
 
 #include <cstddef>
+#include <map>
 #include <optional>
 #include <unordered_map>
 #include <vector>
 
 #include "DataStructures/Tensor/TypeAliases.hpp"
 #include "Domain/BlockLogicalCoordinates.hpp"
+#include "Domain/Structure/ElementSearchTree.hpp"
 
 /// \cond
 namespace domain {
@@ -39,6 +41,29 @@ element_logical_coordinates(
     const tnsr::I<double, Dim, Frame::BlockLogical>& x_block_logical,
     const ElementId<Dim>& element_id);
 
+/*!
+ * \brief Map block logical coordinates to element logical coordinates using an
+ * indexed set of elements
+ *
+ * The indexing is done by the `domain::index_element_ids` function, which
+ * places element IDs into a search tree for fast lookup. If the point is on a
+ * shared element boundary then no guarantee is made which of the abutting
+ * element IDs is returned (the returned element ID is deterministic but depends
+ * on the order in which the elements were inserted into the search tree). See
+ * the other function overload for handling multiple points and disambiguating
+ * points on shared element boundaries.
+ *
+ * \param block_logical_coords the block logical coordinates of the point
+ * \param search_tree the result of `domain::index_element_ids`
+ */
+template <size_t Dim>
+std::optional<
+    std::pair<ElementId<Dim>, tnsr::I<double, Dim, Frame::ElementLogical>>>
+element_logical_coordinates(
+    const IdPair<domain::BlockId, tnsr::I<double, Dim, Frame::BlockLogical>>&
+        block_logical_coords,
+    const std::map<size_t, domain::ElementSearchTree<Dim>>& search_tree);
+
 /// \ingroup ComputationalDomainGroup
 ///
 /// Holds element logical coordinates of an arbitrary set of points on

diff --git a/src/Domain/Structure/CMakeLists.txt b/src/Domain/Structure/CMakeLists.txt
@@ -45,6 +45,7 @@ spectre_target_headers(
   DirectionMap.hpp
   Element.hpp
   ElementId.hpp
+  ElementSearchTree.hpp
   HasBoundary.hpp
   Hypercube.hpp
   IndexToSliceAt.hpp

diff --git a/src/Domain/Structure/ElementSearchTree.hpp b/src/Domain/Structure/ElementSearchTree.hpp
@@ -0,0 +1,122 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <boost/geometry.hpp>
+#include <boost/geometry/index/rtree.hpp>
+#include <cstddef>
+#include <map>
+#include <utility>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/Structure/ElementId.hpp"
+
+namespace boost::geometry::traits {
+// Make Tensor compatible with Boost.Geometry
+// This is needed to search for a block-logical coordinate given as a Tensor
+// in the `ElementSearchTree`.
+template <typename DataType, size_t Dim, typename Frame>
+struct tag<tnsr::I<DataType, Dim, Frame>> {
+  using type = boost::geometry::point_tag;
+};
+template <typename DataType, size_t Dim, typename Frame>
+struct coordinate_type<tnsr::I<DataType, Dim, Frame>> {
+  using type = DataType;
+};
+template <typename DataType, size_t Dim, typename Frame>
+struct coordinate_system<tnsr::I<DataType, Dim, Frame>> {
+  using type = boost::geometry::cs::cartesian;
+};
+template <typename DataType, size_t Dim, typename Frame>
+struct dimension<tnsr::I<DataType, Dim, Frame>>
+    : std::integral_constant<size_t, Dim> {};
+template <typename DataType, size_t Dim, typename Frame, size_t GetDimension>
+struct access<tnsr::I<DataType, Dim, Frame>, GetDimension> {
+  static constexpr DataType get(const tnsr::I<DataType, Dim, Frame>& point) {
+    return ::get<GetDimension>(point);
+  }
+  static void set(tnsr::I<DataType, Dim, Frame>& point, const DataType& value) {
+    ::get<GetDimension>(point) = value;
+  }
+};
+// Make ElementId compatible with Boost.Geometry
+// Each ElementId defines a bounding box in block-logical coordinates, which is
+// used to search for elements in the `ElementSearchTree`.
+template <size_t Dim>
+struct tag<ElementId<Dim>> {
+  using type = boost::geometry::box_tag;
+};
+template <size_t Dim>
+struct coordinate_type<ElementId<Dim>> {
+  using type = double;
+};
+template <size_t Dim>
+struct coordinate_system<ElementId<Dim>> {
+  using type = boost::geometry::cs::cartesian;
+};
+template <size_t Dim>
+struct dimension<ElementId<Dim>> : std::integral_constant<size_t, Dim> {};
+template <size_t Dim>
+struct point_type<ElementId<Dim>> {
+  using type = tnsr::I<double, Dim, ::Frame::BlockLogical>;
+};
+template <size_t Dim, size_t Index, size_t GetDimension>
+struct indexed_access<ElementId<Dim>, Index, GetDimension> {
+  static constexpr double get(const ElementId<Dim>& element_id) {
+    if constexpr (Index == 0) {
+      return element_id.segment_id(GetDimension).endpoint(Side::Lower);
+    } else {
+      return element_id.segment_id(GetDimension).endpoint(Side::Upper);
+    }
+  }
+};
+}  // namespace boost::geometry::traits
+
+namespace domain {
+
+/*!
+ * \brief Search tree for efficiently looking up elements by their bounding
+ * boxes in block-logical coordinates.
+ *
+ * The search tree is constructed from a list of `ElementId`s, which define
+ * bounding boxes in block-logical coordinates. All `ElementId`s must be in the
+ * same block. Then, the search tree can be used to efficiently search for the
+ * element that contains a given block-logical coordinate (see e.g.
+ * `element_logical_coordinates`).
+ *
+ * Example usage:
+ * \snippet Test_ElementSearchTree.cpp element_search_tree_example
+ *
+ * Use `domain::index_element_ids()` to create one search tree for each block in
+ * the domain given the full list of `ElementId`s.
+ *
+ * \details The search tree is a `boost::geometry::rtree` with a choice of
+ * quadratic splitting algorithm and a maximum of 16 elements per node. These
+ * choices work well, but haven't been extensively tuned.
+ */
+template <size_t Dim>
+using ElementSearchTree =
+    boost::geometry::index::rtree<ElementId<Dim>,
+                                  boost::geometry::index::quadratic<16>>;
+
+/*!
+ * \brief Sorts element IDs into one `ElementSearchTree` per block for efficient
+ * searching.
+ *
+ * Returns a map of search trees indexed by block ID. Each search tree contains
+ * all the `ElementId`s for that block.
+ *
+ * \see ElementSearchTree
+ */
+template <size_t Dim>
+std::map<size_t, ElementSearchTree<Dim>> index_element_ids(
+    const std::vector<ElementId<Dim>>& element_ids) {
+  std::map<size_t, ElementSearchTree<Dim>> trees{};
+  for (const auto& element_id : element_ids) {
+    trees[element_id.block_id()].insert(element_id);
+  }
+  return trees;
+}
+
+}  // namespace domain
diff --git a/src/IO/Exporter/PointwiseInterpolator.cpp b/src/IO/Exporter/PointwiseInterpolator.cpp
@@ -152,7 +152,7 @@ void interpolate_to_points_impl(
 template <size_t Dim>
 void interpolate_to_point_impl(
     const gsl::not_null<std::vector<double>*> result,
-    const tnsr::I<DataVector, Dim, Frame::ElementLogical>& x_element_logical,
+    const tnsr::I<double, Dim, Frame::ElementLogical>& x_element_logical,
     const Mesh<Dim>& mesh, const size_t offset, const size_t length,
     const std::vector<DataVector>& tensor_data) {
   const intrp::Irregular<Dim> interpolant(mesh, x_element_logical);
@@ -584,14 +584,21 @@ void interpolate_to_points(
     }
   }
 
-  if (extrapolate_into_excisions) {
+  if (extrapolate_into_excisions and not extrapolation_info.empty()) {
     extrapolate_into_excisions_impl(result, extrapolation_info,
                                     resolved_num_threads);
     // Clear the anchor points from the result
     for (size_t i = 0; i < result->size(); ++i) {
       ResultDataType& component = (*result)[i];
       if constexpr (std::is_same_v<ResultDataType, DataVector>) {
-        component.destructive_resize(num_target_points);
+        // Can't use `destructive_resize` because we want to preserve the
+        // leading elements of the DataVector
+        DataVector new_component(num_target_points);
+        std::copy(
+            component.begin(),
+            component.begin() + static_cast<std::ptrdiff_t>(num_target_points),
+            new_component.begin());
+        component = std::move(new_component);
       } else {
         component.resize(num_target_points);
       }
@@ -615,13 +622,16 @@ PointwiseInterpolator<Dim, Frame>::PointwiseInterpolator(
 
   // Load tensor data into memory
   element_ids_.reserve(filenames.size());
+  element_search_trees_.reserve(filenames.size());
   meshes_.reserve(filenames.size());
   tensor_data_.reserve(filenames.size());
   for (const auto& filename : filenames) {
     const h5::H5File<h5::AccessType::ReadOnly> h5file(filename);
     const auto& volfile = h5file.get<h5::VolumeData>(subfile_name);
     const auto [element_ids, meshes] = load_grids<Dim>(volfile, obs_id_);
     element_ids_.push_back(std::move(element_ids));
+    element_search_trees_.push_back(
+        domain::index_element_ids(element_ids_.back()));
     meshes_.push_back(std::move(meshes));
     tensor_data_.push_back(
         load_tensor_data(volfile, obs_id_, tensor_components));
@@ -708,19 +718,27 @@ void PointwiseInterpolator<Dim, Frame>::interpolate_to_point(
   if (not block_logical_coords.has_value()) {
     ERROR("Point is not in any block:\n" << target_point);
   }
-  // Process all volume files in serial
+  interpolate_to_point(result, block_logical_coords.value());
+}
+
+template <size_t Dim, typename Frame>
+void PointwiseInterpolator<Dim, Frame>::interpolate_to_point(
+    const gsl::not_null<std::vector<double>*> result,
+    const IdPair<domain::BlockId, tnsr::I<double, Dim, ::Frame::BlockLogical>>&
+        target_point) const {
+  ASSERT(not tensor_data_.empty(),
+         "PointwiseInterpolator has not been initialized with tensor data.");
   for (size_t file_id = 0; file_id < element_ids_.size(); ++file_id) {
-    const auto& element_ids = element_ids_[file_id];
+    const auto& element_search_tree = element_search_trees_[file_id];
     const auto& meshes = meshes_[file_id];
     const auto& tensor_data = tensor_data_[file_id];
     const auto element_logical_coords =
-        element_logical_coordinates(element_ids, {block_logical_coords});
-    if (element_logical_coords.empty()) {
+        element_logical_coordinates(target_point, element_search_tree);
+    if (not element_logical_coords.has_value()) {
       continue;
     }
-    const auto& element_id = element_logical_coords.begin()->first;
-    const auto& x_element_logical =
-        element_logical_coords.begin()->second.element_logical_coords;
+    const auto& element_id = element_logical_coords.value().first;
+    const auto& x_element_logical = element_logical_coords.value().second;
     const auto& mesh_offset_length = meshes.at(element_id);
     interpolate_to_point_impl(
         result, x_element_logical, get<0>(mesh_offset_length),

diff --git a/src/IO/Exporter/PointwiseInterpolator.hpp b/src/IO/Exporter/PointwiseInterpolator.hpp
@@ -13,7 +13,9 @@
 
 #include "DataStructures/DataBox/TagName.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/BlockLogicalCoordinates.hpp"
 #include "Domain/Domain.hpp"
+#include "Domain/Structure/ElementSearchTree.hpp"
 #include "IO/Exporter/Exporter.hpp"
 #include "Utilities/TaggedTuple.hpp"
 
@@ -182,6 +184,20 @@ struct PointwiseInterpolator {
                             std::optional<gsl::not_null<std::vector<size_t>*>>
                                 block_order = std::nullopt) const;
 
+  /*!
+   * \brief Interpolate to a single point in block-logical coordinates
+   *
+   * \param result the interpolated data at the target point. The vector is over
+   * the number of components. Will be resized automatically.
+   * \param target_point the point to interpolate to in block-logical
+   * coordinates.
+   */
+  void interpolate_to_point(
+      gsl::not_null<std::vector<double>*> result,
+      const IdPair<domain::BlockId,
+                   tnsr::I<double, Dim, ::Frame::BlockLogical>>& target_point)
+      const;
+
   size_t obs_id() const { return obs_id_; }
   double time() const { return time_; }
   const Domain<Dim>& domain() const { return domain_; }
@@ -197,6 +213,8 @@ struct PointwiseInterpolator {
   domain::FunctionsOfTimeMap functions_of_time_;
   // Outer vector is the source data file
   std::vector<std::vector<ElementId<Dim>>> element_ids_;
+  std::vector<std::map<size_t, domain::ElementSearchTree<Dim>>>
+      element_search_trees_;
   std::vector<
       std::unordered_map<ElementId<Dim>, std::tuple<Mesh<Dim>, size_t, size_t>>>
       meshes_;