Support blocks representing spherical shells

src/Domain/Block.cpp

@@ -27,6 +27,31 @@ Block<VolumeDim>::Block(
     std::string name)
     : stationary_map_(std::move(stationary_map)),
       id_(id),
+      geometry_(domain::BlockGeometry::Cube),
+      name_(std::move(name)) {
+  for (auto& [direction, neighbor] : neighbors) {
+    neighbors_[direction].insert(std::move(neighbor));
+  }
+  // Loop over Directions to search which Directions were not set to neighbors_,
+  // set these Directions to external_boundaries_.
+  for (const auto& direction : Direction<VolumeDim>::all_directions()) {
+    if (neighbors_.find(direction) == neighbors_.end()) {
+      external_boundaries_.emplace(direction);
+    }
+  }
+}
+
+template <size_t VolumeDim>
+Block<VolumeDim>::Block(
+    std::unique_ptr<domain::CoordinateMapBase<
+        Frame::BlockLogical, Frame::Inertial, VolumeDim>>&& stationary_map,
+    const size_t id,
+    DirectionMap<VolumeDim, std::unordered_set<BlockNeighbor<VolumeDim>>>
+        neighbors,
+    std::string name)
+    : stationary_map_(std::move(stationary_map)),
+      id_(id),
+      geometry_(domain::BlockGeometry::SphericalShell),
       neighbors_(std::move(neighbors)),
       name_(std::move(name)) {
   // Loop over Directions to search which Directions were not set to neighbors_,
@@ -118,7 +143,7 @@ void Block<VolumeDim>::inject_time_dependent_map(
 
 template <size_t VolumeDim>
 void Block<VolumeDim>::pup(PUP::er& p) {
-  size_t version = 1;
+  size_t version = 2;
   p | version;
   // Remember to increment the version number when making changes to this
   // function. Retain support for unpacking data written by previous versions
@@ -130,7 +155,18 @@ void Block<VolumeDim>::pup(PUP::er& p) {
     p | moving_mesh_grid_to_distorted_map_;
     p | moving_mesh_distorted_to_inertial_map_;
     p | id_;
-    p | neighbors_;
+    if (version < 2) {
+      geometry_ = domain::BlockGeometry::Cube;
+      DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>> neighbors;
+      p | neighbors;
+      neighbors_.clear();
+      for (auto& [direction, neighbor] : neighbors) {
+        neighbors_[direction].insert(std::move(neighbor));
+      }
+    } else {
+      p | geometry_;
+      p | neighbors_;
+    }
     p | external_boundaries_;
   }
   if (version >= 1) {
@@ -141,6 +177,7 @@ void Block<VolumeDim>::pup(PUP::er& p) {
 template <size_t VolumeDim>
 std::ostream& operator<<(std::ostream& os, const Block<VolumeDim>& block) {
   os << "Block " << block.id() << " (" << block.name() << "):\n";
+  os << "Geometry: " << block.geometry() << '\n';
   os << "Neighbors: " << block.neighbors() << '\n';
   os << "External boundaries: " << block.external_boundaries() << '\n';
   os << "Is time dependent: " << std::boolalpha << block.is_time_dependent();
@@ -150,7 +187,8 @@ std::ostream& operator<<(std::ostream& os, const Block<VolumeDim>& block) {
 template <size_t VolumeDim>
 bool operator==(const Block<VolumeDim>& lhs, const Block<VolumeDim>& rhs) {
   bool blocks_are_equal =
-      (lhs.id() == rhs.id() and lhs.neighbors() == rhs.neighbors() and
+      (lhs.id() == rhs.id() and lhs.geometry() == rhs.geometry() and
+       lhs.neighbors() == rhs.neighbors() and
        lhs.external_boundaries() == rhs.external_boundaries() and
        lhs.name() == rhs.name() and
        lhs.is_time_dependent() == rhs.is_time_dependent() and

src/Domain/Block.hpp

@@ -13,6 +13,7 @@
 #include <unordered_set>
 
 #include "Domain/CoordinateMaps/CoordinateMap.hpp"
+#include "Domain/Structure/BlockGeometry.hpp"
 #include "Domain/Structure/BlockNeighbor.hpp"
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/DirectionMap.hpp"
@@ -33,7 +34,10 @@ class er;
 /// Elements that cover a region of the computational domain.
 ///
 /// Each codimension 1 boundary of a Block<VolumeDim> is either an external
-/// boundary or identical to a boundary of one other Block.
+/// boundary or an internal boundary to one or more neighboring blocks. The only
+/// currently supported case where blocks can have more than one neighbor are
+/// interfaces between a spherical shell and wedges. In all other cases the
+/// internal boundaries between neighboring blocks must be identical.
 ///
 /// A Block has logical coordinates that go from -1 to +1 in each
 /// dimension.  The global coordinates are obtained from the logical
@@ -44,6 +48,8 @@ class er;
 template <size_t VolumeDim>
 class Block {
  public:
+  /// Block with one neighbor per direction. Currently always a deformed cube.
+  ///
   /// \param stationary_map the CoordinateMap.
   /// \param id a unique ID.
   /// \param neighbors info about the Blocks that share a codimension 1
@@ -54,13 +60,24 @@ class Block {
         size_t id, DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>> neighbors,
         std::string name = "");
 
+  /// Block with multiple neighbors per direction. Currently only supports
+  /// spherical shells.
+  Block(std::unique_ptr<domain::CoordinateMapBase<
+            Frame::BlockLogical, Frame::Inertial, VolumeDim>>&& stationary_map,
+        size_t id,
+        DirectionMap<VolumeDim, std::unordered_set<BlockNeighbor<VolumeDim>>>
+            neighbors,
+        std::string name = "");
+
   Block() = default;
   ~Block() = default;
   Block(const Block&) = delete;
   Block(Block&&) = default;
   Block& operator=(const Block&) = delete;
   Block& operator=(Block&&) = default;
 
+  domain::BlockGeometry geometry() const { return geometry_; }
+
   /// \brief The map used when the coordinate map is time-independent.
   ///
   /// \see is_time_dependent()
@@ -148,7 +165,8 @@ class Block {
   size_t id() const { return id_; }
 
   /// Information about the neighboring Blocks.
-  const DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>>& neighbors() const {
+  const DirectionMap<VolumeDim, std::unordered_set<BlockNeighbor<VolumeDim>>>&
+  neighbors() const {
     return neighbors_;
   }
 
@@ -186,7 +204,9 @@ class Block {
       moving_mesh_distorted_to_inertial_map_{nullptr};
 
   size_t id_{0};
-  DirectionMap<VolumeDim, BlockNeighbor<VolumeDim>> neighbors_;
+  domain::BlockGeometry geometry_{domain::BlockGeometry::Cube};
+  DirectionMap<VolumeDim, std::unordered_set<BlockNeighbor<VolumeDim>>>
+      neighbors_;
   std::unordered_set<Direction<VolumeDim>> external_boundaries_;
   std::string name_;
 };

src/Domain/CreateInitialElement.cpp

@@ -12,6 +12,7 @@
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/Element.hpp"
 #include "Domain/Structure/ElementId.hpp"
+#include "Domain/Structure/InitialElementIds.hpp"
 #include "Domain/Structure/Neighbors.hpp"
 #include "Domain/Structure/OrientationMap.hpp"
 #include "Domain/Structure/SegmentId.hpp"
@@ -25,18 +26,82 @@
 namespace domain::Initialization {
 template <size_t VolumeDim>
 Element<VolumeDim> create_initial_element(
-    const ElementId<VolumeDim>& element_id, const Block<VolumeDim>& block,
+    const ElementId<VolumeDim>& element_id,
+    const std::vector<Block<VolumeDim>>& blocks,
     const std::vector<std::array<size_t, VolumeDim>>&
         initial_refinement_levels) {
+  const auto& block = blocks[element_id.block_id()];
   const auto& neighbors_of_block = block.neighbors();
   const auto segment_ids = element_id.segment_ids();
 
   // Declare two helper lambdas for setting the neighbors of an element
   const auto compute_element_neighbor_in_other_block =
-      [&block, &initial_refinement_levels, &neighbors_of_block, &segment_ids,
-       grid_index =
-           element_id.grid_index()](const Direction<VolumeDim>& direction) {
-        const auto& block_neighbor = neighbors_of_block.at(direction);
+      [&block, &blocks, &initial_refinement_levels, &neighbors_of_block,
+       &segment_ids, grid_index = element_id.grid_index()](
+          const Direction<VolumeDim>& direction) {
+        const auto& block_neighbors = neighbors_of_block.at(direction);
+
+        // Handle spherical shells. We assume that spherical shells have no
+        // angular h-refinement. That makes the logic here quite simple:
+        // The single-element shell is the neighbor of all radial neighbors.
+        if (block.geometry() == domain::BlockGeometry::SphericalShell or
+            (block_neighbors.size() == 1 and
+             blocks[block_neighbors.begin()->id()].geometry() ==
+                 domain::BlockGeometry::SphericalShell)) {
+#ifdef SPECTRE_DEBUG
+          const auto check_angular_refinement =
+              [&initial_refinement_levels,
+               &direction](const Block<VolumeDim>& test_block) {
+                if (test_block.geometry() !=
+                    domain::BlockGeometry::SphericalShell) {
+                  return;
+                }
+                const auto& refinement_levels =
+                    initial_refinement_levels[test_block.id()];
+                for (size_t d = 0; d < VolumeDim; ++d) {
+                  if (d == direction.dimension()) {
+                    continue;
+                  }
+                  ASSERT(refinement_levels[d] == 0,
+                         "Spherical shells are assumed here to have no angular "
+                         "refinement in angular directions.");
+                }
+              };
+          check_angular_refinement(block);
+#endif  // SPECTRE_DEBUG
+          std::unordered_set<ElementId<VolumeDim>> neighbor_ids;
+          for (const auto& block_neighbor : block_neighbors) {
+#ifdef SPECTRE_DEBUG
+            check_angular_refinement(blocks[block_neighbor.id()]);
+#endif  // SPECTRE_DEBUG
+            const auto& orientation = block_neighbor.orientation();
+            const auto direction_from_neighbor =
+                orientation(direction.opposite());
+            const auto& refinement_of_neighbor =
+                initial_refinement_levels[block_neighbor.id()];
+            const size_t index_of_neighbor =
+                direction_from_neighbor.side() == Side::Lower
+                    ? 0
+                    : (two_to_the(refinement_of_neighbor[direction_from_neighbor
+                                                             .dimension()]) -
+                       1);
+            for (const auto& neighbor_id : initial_element_ids(
+                     block_neighbor.id(), refinement_of_neighbor, grid_index)) {
+              if (neighbor_id.segment_id(direction_from_neighbor.dimension())
+                      .index() == index_of_neighbor) {
+                neighbor_ids.insert(neighbor_id);
+              }
+            }
+          }
+          return std::make_pair(
+              direction, Neighbors<VolumeDim>(
+                             std::move(neighbor_ids),
+                             // TODO: may have to set this orientation
+                             OrientationMap<VolumeDim>::create_aligned(),
+                             blocks[block_neighbors.begin()->id()].geometry()));
+        }
+
+        const auto& block_neighbor = *block_neighbors.begin();
         const auto& orientation = block_neighbor.orientation();
         const auto direction_in_neighbor = orientation(direction);
 
@@ -131,29 +196,28 @@ Element<VolumeDim> create_initial_element(
         }
         return std::make_pair(
             direction,
-            Neighbors<VolumeDim>(std::move(neighbor_ids), orientation));
+            Neighbors<VolumeDim>(std::move(neighbor_ids), orientation,
+                                 blocks[block_neighbor.id()].geometry()));
       };
 
-  const auto compute_element_neighbor_in_same_block = [&element_id,
-                                                       &segment_ids](
-                                                          const Direction<
-                                                              VolumeDim>&
-                                                              direction) {
-    auto segment_ids_of_neighbor = segment_ids;
-    auto& perpendicular_segment_id =
-        gsl::at(segment_ids_of_neighbor, direction.dimension());
-    const auto index = perpendicular_segment_id.index();
-    perpendicular_segment_id =
-        SegmentId(perpendicular_segment_id.refinement_level(),
-                  direction.side() == Side::Upper ? index + 1 : index - 1);
-    return std::make_pair(
-        direction,
-        Neighbors<VolumeDim>(
-            {{ElementId<VolumeDim>{element_id.block_id(),
-                                   std::move(segment_ids_of_neighbor),
-                                   element_id.grid_index()}}},
-            OrientationMap<VolumeDim>::create_aligned()));
-  };
+  const auto compute_element_neighbor_in_same_block =
+      [&element_id, &segment_ids,
+       geometry = block.geometry()](const Direction<VolumeDim>& direction) {
+        auto segment_ids_of_neighbor = segment_ids;
+        auto& perpendicular_segment_id =
+            gsl::at(segment_ids_of_neighbor, direction.dimension());
+        const auto index = perpendicular_segment_id.index();
+        perpendicular_segment_id =
+            SegmentId(perpendicular_segment_id.refinement_level(),
+                      direction.side() == Side::Upper ? index + 1 : index - 1);
+        return std::make_pair(
+            direction,
+            Neighbors<VolumeDim>(
+                {{ElementId<VolumeDim>{element_id.block_id(),
+                                       std::move(segment_ids_of_neighbor),
+                                       element_id.grid_index()}}},
+                OrientationMap<VolumeDim>::create_aligned(), geometry));
+      };
 
   typename Element<VolumeDim>::Neighbors_t neighbors_of_element;
   for (size_t d = 0; d < VolumeDim; ++d) {
@@ -186,10 +250,10 @@ Element<VolumeDim> create_initial_element(
 
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
-#define INSTANTIATE(_, data)                                 \
-  template Element<DIM(data)>                                \
-  domain::Initialization::create_initial_element<DIM(data)>( \
-      const ElementId<DIM(data)>&, const Block<DIM(data)>&,  \
+#define INSTANTIATE(_, data)                                             \
+  template Element<DIM(data)>                                            \
+  domain::Initialization::create_initial_element<DIM(data)>(             \
+      const ElementId<DIM(data)>&, const std::vector<Block<DIM(data)>>&, \
       const std::vector<std::array<size_t, DIM(data)>>&);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))

src/Domain/CreateInitialElement.hpp

@@ -29,7 +29,8 @@ namespace Initialization {
  */
 template <size_t VolumeDim>
 Element<VolumeDim> create_initial_element(
-    const ElementId<VolumeDim>& element_id, const Block<VolumeDim>& block,
+    const ElementId<VolumeDim>& element_id,
+    const std::vector<Block<VolumeDim>>& blocks,
     const std::vector<std::array<size_t, VolumeDim>>&
         initial_refinement_levels);
 }  // namespace Initialization

src/Domain/ElementDistribution.cpp

@@ -57,14 +57,15 @@ namespace {
 // local time stepping.
 template <size_t Dim>
 double get_num_points_and_grid_spacing_cost(
-    const ElementId<Dim>& element_id, const Block<Dim>& block,
+    const ElementId<Dim>& element_id, const std::vector<Block<Dim>>& blocks,
     const std::vector<std::array<size_t, Dim>>& initial_refinement_levels,
     const std::vector<std::array<size_t, Dim>>& initial_extents,
     const Spectral::Quadrature quadrature) {
+  const auto& block = blocks[element_id.block_id()];
   Mesh<Dim> mesh = ::domain::Initialization::create_initial_mesh(
       initial_extents, element_id, quadrature);
   Element<Dim> element = ::domain::Initialization::create_initial_element(
-      element_id, block, initial_refinement_levels);
+      element_id, blocks, initial_refinement_levels);
   ElementMap<Dim, Frame::Grid> element_map{element_id, block};
   const tnsr::I<DataVector, Dim, Frame::ElementLogical> logical_coords =
       logical_coordinates(mesh);
@@ -122,7 +123,7 @@ std::unordered_map<ElementId<Dim>, double> get_element_costs(
 
         element_costs.insert(
             {element_id, get_num_points_and_grid_spacing_cost(
-                             element_id, block, initial_refinement_levels,
+                             element_id, blocks, initial_refinement_levels,
                              initial_extents, quadrature.value())});
       }
     }
@@ -323,7 +324,7 @@ size_t BlockZCurveProcDistribution<Dim>::get_proc_for_element(
   template class BlockZCurveProcDistribution<GET_DIM(data)>;                 \
   double get_num_points_and_grid_spacing_cost(                               \
       const ElementId<GET_DIM(data)>& element_id,                            \
-      const Block<GET_DIM(data)>& block,                                     \
+      const std::vector<Block<GET_DIM(data)>>& blocks,                       \
       const std::vector<std::array<size_t, GET_DIM(data)>>&                  \
           initial_refinement_levels,                                         \
       const std::vector<std::array<size_t, GET_DIM(data)>>& initial_extents, \