Support spherical shell blocks, create Sphere domain

src/Domain/Block.cpp

@@ -27,6 +27,31 @@
     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 @@
 
 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 @@
     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 @@
 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 @@
 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/CMakeLists.txt

@@ -40,6 +40,7 @@ spectre_target_headers(
   CreateInitialElement.hpp
   Domain.hpp
   DomainHelpers.hpp
+  DomainHelpers.tpp
   ElementDistribution.hpp
   ElementLogicalCoordinates.hpp
   ElementMap.hpp

src/Domain/CoordinateMaps/CMakeLists.txt

@@ -31,6 +31,7 @@ spectre_target_sources(
   Interval.cpp
   KerrHorizonConforming.cpp
   Rotation.cpp
+  ShellType.cpp
   SpecialMobius.cpp
   SphericalToCartesianPfaffian.cpp
   UniformCylindricalEndcap.cpp
@@ -72,6 +73,7 @@ spectre_target_headers(
   ProductMaps.hpp
   ProductMaps.tpp
   Rotation.hpp
+  ShellType.hpp
   SpecialMobius.hpp
   SphericalToCartesianPfaffian.hpp
   Tags.hpp

src/Domain/CoordinateMaps/ShellType.cpp

@@ -0,0 +1,39 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Domain/CoordinateMaps/ShellType.hpp"
+
+#include <ostream>
+#include <string>
+
+#include "Options/Options.hpp"
+#include "Options/ParseOptions.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+
+namespace domain::CoordinateMaps {
+
+std::ostream& operator<<(std::ostream& os, const ShellType shell_type) {
+  switch (shell_type) {
+    case ShellType::Cubed:
+      return os << "Cubed";
+    case ShellType::Spherical:
+      return os << "Spherical";
+    default:
+      ERROR("Unknown domain::CoordinateMaps::ShellType");
+  }
+}
+
+}  // namespace domain::CoordinateMaps
+
+template <>
+domain::CoordinateMaps::ShellType
+Options::create_from_yaml<domain::CoordinateMaps::ShellType>::create<void>(
+    const Options::Option& options) {
+  const auto shell_type = options.parse_as<std::string>();
+  if (shell_type == "Cubed") {
+    return domain::CoordinateMaps::ShellType::Cubed;
+  } else if (shell_type == "Spherical") {
+    return domain::CoordinateMaps::ShellType::Spherical;
+  }
+  PARSE_ERROR(options.context(), "ShellType must be 'Cubed' or 'Spherical'.");
+}

src/Domain/CoordinateMaps/ShellType.hpp

@@ -0,0 +1,44 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <ostream>
+
+/// \cond
+namespace Options {
+class Option;
+template <typename T>
+struct create_from_yaml;
+}  // namespace Options
+/// \endcond
+
+namespace domain::CoordinateMaps {
+
+/*!
+ * \brief Type of spherical shell: built from four (2D) or six (3D) wedges
+ * ("cubed") or from a single spherical shell ("spherical").
+ *
+ * Used to select the shell type in the input file.
+ */
+enum class ShellType {
+  Cubed,
+  Spherical,
+};
+
+std::ostream& operator<<(std::ostream& os, ShellType shell_type);
+
+}  // namespace domain::CoordinateMaps
+
+template <>
+struct Options::create_from_yaml<domain::CoordinateMaps::ShellType> {
+  template <typename Metavariables>
+  static domain::CoordinateMaps::ShellType create(
+      const Options::Option& options) {
+    return create<void>(options);
+  }
+};
+template <>
+domain::CoordinateMaps::ShellType
+Options::create_from_yaml<domain::CoordinateMaps::ShellType>::create<void>(
+    const Options::Option& options);