From a9e43057f5b398e29b299538ae72dc15aebca80d Mon Sep 17 00:00:00 2001
From: Nils Vu <owls@nilsvu.de>
Date: Mon, 9 Dec 2024 08:57:54 -0800
Subject: [PATCH] Create spherical shell in Sphere domain

---
 src/Domain/CMakeLists.txt                     |   1 +
 src/Domain/CoordinateMaps/CMakeLists.txt      |   2 +
 src/Domain/CoordinateMaps/ShellType.cpp       |  39 +++++
 src/Domain/CoordinateMaps/ShellType.hpp       |  44 ++++++
 src/Domain/Creators/BinaryCompactObject.cpp   |  60 ++++---
 src/Domain/Creators/Sphere.cpp                | 147 ++++++++++++++----
 src/Domain/Creators/Sphere.hpp                |  33 +++-
 src/Domain/Domain.cpp                         |  10 +-
 src/Domain/Domain.hpp                         |   7 +-
 src/Domain/DomainHelpers.cpp                  | 143 +----------------
 src/Domain/DomainHelpers.hpp                  |  11 +-
 src/Domain/DomainHelpers.tpp                  | 140 +++++++++++++++++
 .../ExportCoordinates/ExportCoordinates.hpp   |  62 ++++----
 .../InputFiles/ExportCoordinates/Input3D.yaml |  24 +--
 14 files changed, 469 insertions(+), 254 deletions(-)
 create mode 100644 src/Domain/CoordinateMaps/ShellType.cpp
 create mode 100644 src/Domain/CoordinateMaps/ShellType.hpp
 create mode 100644 src/Domain/DomainHelpers.tpp

diff --git a/src/Domain/CMakeLists.txt b/src/Domain/CMakeLists.txt
index 39833400bf6a..693bd08fc694 100644
--- a/src/Domain/CMakeLists.txt
+++ b/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
diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt
index 17fa590c14d8..8a495ecdc170 100644
--- a/src/Domain/CoordinateMaps/CMakeLists.txt
+++ b/src/Domain/CoordinateMaps/CMakeLists.txt
@@ -31,6 +31,7 @@ spectre_target_sources(
   Interval.cpp
   KerrHorizonConforming.cpp
   Rotation.cpp
+  ShellType.cpp
   SpecialMobius.cpp
   SphericalToCartesian.cpp
   UniformCylindricalEndcap.cpp
@@ -72,6 +73,7 @@ spectre_target_headers(
   ProductMaps.hpp
   ProductMaps.tpp
   Rotation.hpp
+  ShellType.hpp
   SpecialMobius.hpp
   SphericalToCartesian.hpp
   Tags.hpp
diff --git a/src/Domain/CoordinateMaps/ShellType.cpp b/src/Domain/CoordinateMaps/ShellType.cpp
new file mode 100644
index 000000000000..7da241165741
--- /dev/null
+++ b/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'.");
+}
diff --git a/src/Domain/CoordinateMaps/ShellType.hpp b/src/Domain/CoordinateMaps/ShellType.hpp
new file mode 100644
index 000000000000..a8f42112b75e
--- /dev/null
+++ b/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);
diff --git a/src/Domain/Creators/BinaryCompactObject.cpp b/src/Domain/Creators/BinaryCompactObject.cpp
index cda5a42579d2..50691a05ab1f 100644
--- a/src/Domain/Creators/BinaryCompactObject.cpp
+++ b/src/Domain/Creators/BinaryCompactObject.cpp
@@ -36,6 +36,7 @@
 #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp"
 #include "Domain/Domain.hpp"
 #include "Domain/DomainHelpers.hpp"
+#include "Domain/DomainHelpers.tpp"
 #include "Domain/ExcisionSphere.hpp"
 #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp"
 #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp"
@@ -459,20 +460,18 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
                   length_inner_cube_ * 0.5,
                   std::array<double, 3>{{offset_x_coord_a_, 0.0, 0.0}}));
     Maps maps_center_A =
-        domain::make_vector_coordinate_map_base<Frame::BlockLogical,
-                                                Frame::Inertial, 3>(
-            sph_wedge_coordinate_maps(
-                object_a.inner_radius, object_a.outer_radius,
-                inner_sphericity_A, 1.0, use_equiangular_map_,
-                offset_a_optional, false, {}, object_A_radial_distribution),
-            translation_A);
+        spherical_shells_coordinate_maps<Frame::BlockLogical, Frame::Inertial>(
+            object_a.inner_radius, object_a.outer_radius, inner_sphericity_A,
+            1.0, use_equiangular_map_, false, {}, object_A_radial_distribution,
+            {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All,
+            M_PI_2, translation_A);
     Maps maps_cube_A =
-        domain::make_vector_coordinate_map_base<Frame::BlockLogical,
-                                                Frame::Inertial, 3>(
-            sph_wedge_coordinate_maps(
-                object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_,
-                1.0, 0.0, use_equiangular_map_, offset_a_optional),
-            translation_A);
+        spherical_shells_coordinate_maps<Frame::BlockLogical, Frame::Inertial>(
+            object_a.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
+            0.0, use_equiangular_map_, false, {},
+            {domain::CoordinateMaps::Distribution::Linear},
+            {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All,
+            M_PI_2, translation_A);
     std::move(maps_center_A.begin(), maps_center_A.end(),
               std::back_inserter(maps));
     std::move(maps_cube_A.begin(), maps_cube_A.end(), std::back_inserter(maps));
@@ -502,20 +501,18 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
                   length_inner_cube_ * 0.5,
                   std::array<double, 3>{{offset_x_coord_b_, 0.0, 0.0}}));
     Maps maps_center_B =
-        domain::make_vector_coordinate_map_base<Frame::BlockLogical,
-                                                Frame::Inertial, 3>(
-            sph_wedge_coordinate_maps(
-                object_b.inner_radius, object_b.outer_radius,
-                inner_sphericity_B, 1.0, use_equiangular_map_,
-                offset_b_optional, false, {}, object_B_radial_distribution),
-            translation_B);
+        spherical_shells_coordinate_maps<Frame::BlockLogical, Frame::Inertial>(
+            object_b.inner_radius, object_b.outer_radius, inner_sphericity_B,
+            1.0, use_equiangular_map_, false, {}, object_B_radial_distribution,
+            {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All,
+            M_PI_2, translation_B);
     Maps maps_cube_B =
-        domain::make_vector_coordinate_map_base<Frame::BlockLogical,
-                                                Frame::Inertial, 3>(
-            sph_wedge_coordinate_maps(
-                object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_,
-                1.0, 0.0, use_equiangular_map_, offset_b_optional),
-            translation_B);
+        spherical_shells_coordinate_maps<Frame::BlockLogical, Frame::Inertial>(
+            object_b.outer_radius, sqrt(3.0) * 0.5 * length_inner_cube_, 1.0,
+            0.0, use_equiangular_map_, false, {},
+            {domain::CoordinateMaps::Distribution::Linear},
+            {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All,
+            M_PI_2, translation_B);
     std::move(maps_center_B.begin(), maps_center_B.end(),
               std::back_inserter(maps));
     std::move(maps_cube_B.begin(), maps_cube_B.end(), std::back_inserter(maps));
@@ -545,12 +542,11 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
 
   // --- Outer spherical shell (10 blocks) ---
   Maps maps_outer_shell =
-      domain::make_vector_coordinate_map_base<Frame::BlockLogical,
-                                              Frame::Inertial, 3>(
-          sph_wedge_coordinate_maps(envelope_radius_, outer_radius_, 1.0, 1.0,
-                                    use_equiangular_map_, std::nullopt, true,
-                                    {}, {radial_distribution_outer_shell_},
-                                    ShellWedges::All, opening_angle_));
+      spherical_shells_coordinate_maps<Frame::BlockLogical, Frame::Inertial>(
+          envelope_radius_, outer_radius_, 1.0, 1.0, use_equiangular_map_, true,
+          {}, {radial_distribution_outer_shell_},
+          {domain::CoordinateMaps::ShellType::Cubed}, ShellWedges::All,
+          opening_angle_);
   std::move(maps_outer_shell.begin(), maps_outer_shell.end(),
             std::back_inserter(maps));
 
diff --git a/src/Domain/Creators/Sphere.cpp b/src/Domain/Creators/Sphere.cpp
index 668afdc5fb22..a4d6335a2e87 100644
--- a/src/Domain/Creators/Sphere.cpp
+++ b/src/Domain/Creators/Sphere.cpp
@@ -30,6 +30,7 @@
 #include "Domain/Creators/TimeDependence/None.hpp"
 #include "Domain/Domain.hpp"
 #include "Domain/DomainHelpers.hpp"
+#include "Domain/DomainHelpers.tpp"
 #include "Domain/Structure/BlockNeighbor.hpp"
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/DirectionMap.hpp"
@@ -37,6 +38,7 @@
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/GetOutput.hpp"
 #include "Utilities/MakeArray.hpp"
+#include "Utilities/Overloader.hpp"
 
 namespace Frame {
 struct Inertial;
@@ -79,7 +81,7 @@ Sphere::Sphere(
     std::optional<EquatorialCompressionOptions> equatorial_compression,
     std::vector<double> radial_partitioning,
     const typename RadialDistribution::type& radial_distribution,
-    ShellWedges which_wedges,
+    const typename ShellType::type& shell_type, ShellWedges which_wedges,
     std::optional<TimeDepOptionType> time_dependent_options,
     std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
         outer_boundary_condition,
@@ -153,29 +155,84 @@ Sphere::Sphere(
                 "Add entries to 'RadialPartitioning' to add outer shells for "
                 "which you can select different radial distributions.");
   }
+  shell_types_ = std::visit(
+      Overloader{
+          [&num_shells = num_shells_](
+              const domain::CoordinateMaps::ShellType& uniform_shell_type) {
+            return std::vector<domain::CoordinateMaps::ShellType>(
+                num_shells, uniform_shell_type);
+          },
+          [](const std::vector<domain::CoordinateMaps::ShellType>&
+                 shell_types) { return shell_types; }},
+      shell_type);
+  if (shell_types_.size() != num_shells_) {
+    PARSE_ERROR(context,
+                "Specify a 'ShellType' for every spherical shell. You "
+                "specified "
+                    << shell_types_.size() << " items, but the domain has "
+                    << num_shells_ << " shells.");
+  }
+  if (fill_interior_ and
+      shell_types_.front() != domain::CoordinateMaps::ShellType::Cubed) {
+    PARSE_ERROR(context,
+                "The 'ShellType' must be 'Cubed' for the innermost "
+                "shell filled with a cube. "
+                "Add entries to 'RadialPartitioning' to add outer shells for "
+                "which you can select different shell types.");
+  }
+  if (which_wedges_ != ShellWedges::All and
+      shell_types_ !=
+          std::vector<domain::CoordinateMaps::ShellType>{
+              num_shells_, domain::CoordinateMaps::ShellType::Cubed}) {
+    PARSE_ERROR(context,
+                "To specify 'ShellWedges' other than 'All', you must specify "
+                "'ShellType' as 'Cubed' for all shells.");
+  }
 
   // Create grid anchors
   grid_anchors_["Center"] =
       tnsr::I<double, 3, Frame::Grid>{std::array{0.0, 0.0, 0.0}};
 
   // Determine number of blocks
-  num_blocks_per_shell_ = which_wedges_ == ShellWedges::All             ? 6
-                          : which_wedges_ == ShellWedges::FourOnEquator ? 4
-                                                                        : 1;
-  num_blocks_ = num_blocks_per_shell_ * num_shells_ + (fill_interior_ ? 1 : 0);
+  num_blocks_per_cubed_shell_ = which_wedges_ == ShellWedges::All ? 6
+                                : which_wedges_ == ShellWedges::FourOnEquator
+                                    ? 4
+                                    : 1;
+  num_blocks_ = alg::accumulate(
+      shell_types_, 0_st,
+      [num_blocks_per_cubed_shell = num_blocks_per_cubed_shell_](
+          const size_t num_blocks,
+          const domain::CoordinateMaps::ShellType shell_type) {
+        if (shell_type == domain::CoordinateMaps::ShellType::Cubed) {
+          return num_blocks + num_blocks_per_cubed_shell;
+        } else if (shell_type == domain::CoordinateMaps::ShellType::Spherical) {
+          return num_blocks + 1;
+        } else {
+          ERROR("Unknown ShellType");
+        }
+      });
+  if (fill_interior_) {
+    ++num_blocks_;
+  }
 
   // Create block names and groups
   static std::array<std::string, 6> wedge_directions{
       "UpperZ", "LowerZ", "UpperY", "LowerY", "UpperX", "LowerX"};
   for (size_t shell = 0; shell < num_shells_; ++shell) {
     std::string shell_prefix = "Shell" + std::to_string(shell);
-    for (size_t direction = which_wedge_index(which_wedges_); direction < 6;
-         ++direction) {
-      const std::string wedge_name =
-          shell_prefix + gsl::at(wedge_directions, direction);
-      block_names_.emplace_back(wedge_name);
-      block_groups_[shell_prefix].insert(wedge_name);
-      block_groups_["Wedges"].insert(wedge_name);
+    if (shell_types_[shell] == domain::CoordinateMaps::ShellType::Cubed) {
+      for (size_t direction = which_wedge_index(which_wedges_); direction < 6;
+           ++direction) {
+        const std::string wedge_name =
+            shell_prefix + gsl::at(wedge_directions, direction);
+        block_names_.emplace_back(wedge_name);
+        block_groups_[shell_prefix].insert(wedge_name);
+        block_groups_["Wedges"].insert(wedge_name);
+      }
+    } else if (shell_types_[shell] ==
+               domain::CoordinateMaps::ShellType::Spherical) {
+      block_names_.emplace_back(shell_prefix);
+      block_groups_["SphericalShells"].insert(shell_prefix);
     }
   }
   if (fill_interior_) {
@@ -264,10 +321,9 @@ Sphere::Sphere(
 }
 
 Domain<3> Sphere::create_domain() const {
-  std::vector<std::array<size_t, 8>> corners =
-      corners_for_radially_layered_domains(num_shells_, fill_interior_,
-                                           {{1, 2, 3, 4, 5, 6, 7, 8}},
-                                           which_wedges_);
+  std::vector<Block<3>> blocks{};
+  blocks.reserve(num_blocks_);
+
   double aspect_ratio = 1.0;
   size_t index_polar_axis = 2;
   if (equatorial_compression_.has_value()) {
@@ -277,14 +333,39 @@ Domain<3> Sphere::create_domain() const {
   const domain::CoordinateMaps::EquatorialCompression compression{
       aspect_ratio, index_polar_axis};
 
-  auto coord_maps = domain::make_vector_coordinate_map_base<Frame::BlockLogical,
-                                                            Frame::Inertial, 3>(
-      sph_wedge_coordinate_maps(
+  auto coord_maps =
+      spherical_shells_coordinate_maps<Frame::BlockLogical, Frame::Inertial>(
           inner_radius_, outer_radius_,
           fill_interior_ ? std::get<InnerCube>(interior_).sphericity : 1.0, 1.0,
-          use_equiangular_map_, std::nullopt, false, radial_partitioning_,
-          radial_distribution_, which_wedges_),
-      compression);
+          use_equiangular_map_, false, radial_partitioning_,
+          radial_distribution_, shell_types_, which_wedges_, M_PI, compression);
+
+  // Assuming only the outer shell is a spherical shell
+  const auto corners_of_wedges = corners_for_radially_layered_domains(
+      num_shells_ - 1, fill_interior_, {{1, 2, 3, 4, 5, 6, 7, 8}},
+      which_wedges_);
+  std::vector<DirectionMap<3, BlockNeighbor<3>>> neighbors_of_wedges{};
+  set_internal_boundaries<3>(make_not_null(&neighbors_of_wedges),
+                             corners_of_wedges);
+  std::unordered_set<BlockNeighbor<3>> neighbors_of_outer_shell{};
+  for (size_t i = 0; i < coord_maps.size() - 1; i++) {
+    if (fill_interior_ and i < 7) {
+      neighbors_of_wedges[i][Direction<3>::lower_zeta()] = BlockNeighbor<3>{
+          num_blocks_ - 1,
+          neighbors_of_wedges[i][Direction<3>::lower_zeta()].orientation()};
+    }
+    if (i >= coord_maps.size() - 7) {
+      neighbors_of_wedges[i][Direction<3>::upper_zeta()] =
+          BlockNeighbor<3>{coord_maps.size() - 1, OrientationMap<3>{}};
+      neighbors_of_outer_shell.emplace(i, OrientationMap<3>{});
+    }
+    blocks.emplace_back(std::move(coord_maps[i]), i,
+                        std::move(neighbors_of_wedges[i]), block_names_[i]);
+  }
+  blocks.emplace_back(Block<3>(
+      std::move(coord_maps.back()), coord_maps.size() - 1,
+      {{Direction<3>::lower_zeta(), std::move(neighbors_of_outer_shell)}},
+      block_names_[coord_maps.size() - 1]));
 
   std::unordered_map<std::string, ExcisionSphere<3>> excision_spheres{};
 
@@ -318,6 +399,9 @@ Domain<3> Sphere::create_domain() const {
               BulgedCube{inner_radius_, inner_cube_sphericity,
                          use_equiangular_map_}));
     }
+    blocks.emplace_back(std::move(coord_maps.back()), num_blocks_ - 1,
+                        std::move(neighbors_of_wedges.back()),
+                        block_names_.back());
   } else {
     // Set up excision sphere only for ShellWedges::All
     // - The first 6 blocks enclose the excised sphere, see
@@ -338,8 +422,8 @@ Domain<3> Sphere::create_domain() const {
     }
   }
 
-  Domain<3> domain{std::move(coord_maps),       corners,      {},
-                   std::move(excision_spheres), block_names_, block_groups_};
+  Domain<3> domain{std::move(blocks), std::move(excision_spheres),
+                   block_groups_};
   ASSERT(domain.blocks().size() == num_blocks_,
          "Unexpected number of blocks. Expected "
              << num_blocks_ << " but created " << domain.blocks().size()
@@ -364,16 +448,11 @@ Domain<3> Sphere::create_domain() const {
       const size_t first_block_outer_shell =
           (num_shells_ - 1) * num_blocks_per_shell_;
       for (size_t block_id = 0; block_id < num_blocks_; block_id++) {
-        const bool is_outer_shell =
-            block_id >= first_block_outer_shell and
-            block_id < (first_block_outer_shell + num_blocks_per_shell_);
-        const bool is_inner_cube =
-            fill_interior_ and (block_id == num_blocks_ - 1);
-        // Correct for 'which_wedges' option
-        const size_t shell = block_id / num_blocks_per_shell_;
-        const size_t block_number = shell * 6 +
-                                    which_wedge_index(which_wedges_) +
-                                    block_id % num_blocks_per_shell_;
+        const bool include_distorted_map_in_first_shell =
+            block_id <
+            (shell_types_[0] == domain::CoordinateMaps::ShellType::Cubed
+                 ? num_blocks_per_cubed_shell_
+                 : 1);
         block_maps_grid_to_distorted[block_id] =
             hard_coded_options.grid_to_distorted_map(block_number,
                                                      is_inner_cube);
diff --git a/src/Domain/Creators/Sphere.hpp b/src/Domain/Creators/Sphere.hpp
index f70a2e2dea5b..63f1d10bc9af 100644
--- a/src/Domain/Creators/Sphere.hpp
+++ b/src/Domain/Creators/Sphere.hpp
@@ -15,6 +15,9 @@
 #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
 #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
 #include "Domain/CoordinateMaps/Distribution.hpp"
+#include "Domain/CoordinateMaps/Identity.hpp"
+#include "Domain/CoordinateMaps/ShellType.hpp"
+#include "Domain/CoordinateMaps/SphericalToCartesian.hpp"
 #include "Domain/Creators/DomainCreator.hpp"
 #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
 #include "Domain/Creators/TimeDependentOptions/Sphere.hpp"
@@ -178,10 +181,17 @@ class Sphere : public DomainCreator<3> {
   using Equiangular3D =
       CoordinateMaps::ProductOf3Maps<Equiangular, Equiangular, Equiangular>;
   using BulgedCube = CoordinateMaps::BulgedCube;
+  using LogicalToSphericalMap = domain::CoordinateMaps::ProductOf3Maps<
+      domain::CoordinateMaps::Identity<1>, domain::CoordinateMaps::Identity<1>,
+      domain::CoordinateMaps::Interval>;
 
  public:
   using maps_list = tmpl::append<
       tmpl::list<
+          domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
+                                LogicalToSphericalMap,
+                                domain::CoordinateMaps::SphericalToCartesian<3>,
+                                domain::CoordinateMaps::EquatorialCompression>,
           // Inner cube
           domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
                                 BulgedCube>,
@@ -311,6 +321,22 @@ class Sphere : public DomainCreator<3> {
         "partitions."};
   };
 
+  struct ShellType {
+    using type = std::variant<domain::CoordinateMaps::ShellType,
+                              std::vector<domain::CoordinateMaps::ShellType>>;
+    static constexpr Options::String help = {
+        "Select type of each spherical shell. Specify a list of N+1 shell "
+        "types for N radial partitions, i.e., one for each shell."
+        "The shell type can be 'Cubed' to create six wedges that compose the "
+        "shell, or 'Spherical' to create a single spherical shell with "
+        "spherical harmonics. The 'Spherical' option "
+        "is not widely supported yet, so it is unlikely to work. "
+        "If the interior of the sphere is filled with a cube, the innermost "
+        "shell must have a 'Cubed' shell type. You can also specify just a "
+        "single shell type (not in a vector) which will use the same shell "
+        "type for all partitions."};
+  };
+
   struct WhichWedges {
     using type = ShellWedges;
     static constexpr Options::String help = {
@@ -340,7 +366,7 @@ class Sphere : public DomainCreator<3> {
   using basic_options =
       tmpl::list<InnerRadius, OuterRadius, Interior, InitialRefinement,
                  InitialGridPoints, UseEquiangularMap, EquatorialCompression,
-                 RadialPartitioning, RadialDistribution, WhichWedges,
+                 RadialPartitioning, RadialDistribution, ShellType, WhichWedges,
                  TimeDependentMaps>;
 
   template <typename Metavariables>
@@ -371,6 +397,8 @@ class Sphere : public DomainCreator<3> {
       std::vector<double> radial_partitioning = {},
       const typename RadialDistribution::type& radial_distribution =
           domain::CoordinateMaps::Distribution::Linear,
+      const typename ShellType::type& shell_type =
+          domain::CoordinateMaps::ShellType::Cubed,
       ShellWedges which_wedges = ShellWedges::All,
       std::optional<TimeDepOptionType> time_dependent_options = std::nullopt,
       std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
@@ -428,6 +456,7 @@ class Sphere : public DomainCreator<3> {
   std::optional<EquatorialCompressionOptions> equatorial_compression_{};
   std::vector<double> radial_partitioning_{};
   std::vector<domain::CoordinateMaps::Distribution> radial_distribution_{};
+  std::vector<domain::CoordinateMaps::ShellType> shell_types_{};
   ShellWedges which_wedges_ = ShellWedges::All;
   std::optional<TimeDepOptionType> time_dependent_options_{};
   bool use_hard_coded_maps_{false};
@@ -435,7 +464,7 @@ class Sphere : public DomainCreator<3> {
       outer_boundary_condition_;
   size_t num_shells_;
   size_t num_blocks_;
-  size_t num_blocks_per_shell_;
+  size_t num_blocks_per_cubed_shell_;
   std::vector<std::string> block_names_{};
   std::unordered_map<std::string, std::unordered_set<std::string>>
       block_groups_{};
diff --git a/src/Domain/Domain.cpp b/src/Domain/Domain.cpp
index 1c956f278835..925812b040ee 100644
--- a/src/Domain/Domain.cpp
+++ b/src/Domain/Domain.cpp
@@ -21,8 +21,14 @@ struct BlockLogical;
 }  // namespace Frame
 
 template <size_t VolumeDim>
-Domain<VolumeDim>::Domain(std::vector<Block<VolumeDim>> blocks)
-    : blocks_(std::move(blocks)) {}
+Domain<VolumeDim>::Domain(
+    std::vector<Block<VolumeDim>> blocks,
+    std::unordered_map<std::string, ExcisionSphere<VolumeDim>> excision_spheres,
+    std::unordered_map<std::string, std::unordered_set<std::string>>
+        block_groups)
+    : blocks_(std::move(blocks)),
+      excision_spheres_(std::move(excision_spheres)),
+      block_groups_(std::move(block_groups)) {}
 
 template <size_t VolumeDim>
 Domain<VolumeDim>::Domain(
diff --git a/src/Domain/Domain.hpp b/src/Domain/Domain.hpp
index b2822d62187a..50e48bfc44e5 100644
--- a/src/Domain/Domain.hpp
+++ b/src/Domain/Domain.hpp
@@ -98,7 +98,12 @@ namespace domain {}
 template <size_t VolumeDim>
 class Domain {
  public:
-  explicit Domain(std::vector<Block<VolumeDim>> blocks);
+  explicit Domain(
+      std::vector<Block<VolumeDim>> blocks,
+      std::unordered_map<std::string, ExcisionSphere<VolumeDim>>
+          excision_spheres = {},
+      std::unordered_map<std::string, std::unordered_set<std::string>>
+          block_groups = {});
 
   /*!
    * \brief Create a Domain using CoordinateMaps to encode the Orientations.
diff --git a/src/Domain/DomainHelpers.cpp b/src/Domain/DomainHelpers.cpp
index 8051cf1b281d..3273ec1cca30 100644
--- a/src/Domain/DomainHelpers.cpp
+++ b/src/Domain/DomainHelpers.cpp
@@ -550,6 +550,10 @@ void set_identified_boundaries(
   }
 }
 
+// A Block or Blocks can be wrapped in an outer layer of Blocks surrounding
+// the original Block(s). In the BBH Domain, this occurs several times, using
+// both Wedges and Frustums. The simplest example in which wrapping is used is
+// Sphere, where the central Block is wrapped with six Wedge<3>s.
 std::array<OrientationMap<3>, 6> orientations_for_sphere_wrappings() {
   return {{
       // Upper Z
@@ -592,145 +596,6 @@ size_t which_wedge_index(const ShellWedges& which_wedges) {
   }
 }
 
-std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
-    const double inner_radius, const double outer_radius,
-    const double inner_sphericity, const double outer_sphericity,
-    const bool use_equiangular_map,
-    const std::optional<std::pair<double, std::array<double, 3>>>&
-        offset_options,
-    const bool use_half_wedges, const std::vector<double>& radial_partitioning,
-    const std::vector<domain::CoordinateMaps::Distribution>&
-        radial_distribution,
-    const ShellWedges which_wedges, const double opening_angle) {
-  ASSERT(not use_half_wedges or which_wedges == ShellWedges::All,
-         "If we are using half wedges we must also be using ShellWedges::All.");
-  ASSERT(radial_distribution.size() == 1 + radial_partitioning.size(),
-         "Specify a radial distribution for every spherical shell. You "
-         "specified "
-             << radial_distribution.size() << " items, but the domain has "
-             << 1 + radial_partitioning.size() << " shells.");
-
-  const auto wedge_orientations = orientations_for_sphere_wrappings();
-
-  using Wedge3DMap = domain::CoordinateMaps::Wedge<3>;
-  using Halves = Wedge3DMap::WedgeHalves;
-  std::vector<Wedge3DMap> wedges_for_all_layers{};
-
-  const size_t number_of_layers = 1 + radial_partitioning.size();
-  double temp_inner_radius = inner_radius;
-  double temp_inner_sphericity = inner_sphericity;
-  if (offset_options.has_value()) {
-    for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) {
-      const auto& radial_distribution_this_layer =
-          radial_distribution.at(layer_i);
-      std::optional<double> optional_outer_radius{};
-      if (outer_sphericity != 0.0) {
-        optional_outer_radius = outer_radius;
-      } else {
-        optional_outer_radius = std::nullopt;
-      }
-      // Generate wedges/half-wedges a layer at a time.
-      std::vector<Wedge3DMap> wedges_for_this_layer{};
-      if (not use_half_wedges) {
-        for (size_t face_j = which_wedge_index(which_wedges); face_j < 6;
-             face_j++) {
-          wedges_for_this_layer.emplace_back(
-              temp_inner_radius, optional_outer_radius,
-              offset_options.value().first, offset_options.value().second,
-              gsl::at(wedge_orientations, face_j), use_equiangular_map,
-              Halves::Both, radial_distribution_this_layer);
-        }
-      } else {
-        for (size_t i = 0; i < 4; i++) {
-          wedges_for_this_layer.emplace_back(
-              temp_inner_radius, optional_outer_radius,
-              offset_options.value().first, offset_options.value().second,
-              gsl::at(wedge_orientations, i), use_equiangular_map,
-              Halves::LowerOnly, radial_distribution_this_layer);
-          wedges_for_this_layer.emplace_back(
-              temp_inner_radius, optional_outer_radius,
-              offset_options.value().first, offset_options.value().second,
-              gsl::at(wedge_orientations, i), use_equiangular_map,
-              Halves::UpperOnly, radial_distribution_this_layer);
-        }
-        wedges_for_this_layer.emplace_back(
-            temp_inner_radius, optional_outer_radius,
-            offset_options.value().first, offset_options.value().second,
-            gsl::at(wedge_orientations, 4), use_equiangular_map, Halves::Both,
-            radial_distribution_this_layer);
-        wedges_for_this_layer.emplace_back(
-            temp_inner_radius, optional_outer_radius,
-            offset_options.value().first, offset_options.value().second,
-            gsl::at(wedge_orientations, 5), use_equiangular_map, Halves::Both,
-            radial_distribution_this_layer);
-      }
-      for (const auto& wedge : wedges_for_this_layer) {
-        wedges_for_all_layers.push_back(wedge);
-      }
-    }
-  } else {
-    double temp_outer_radius{};
-    for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) {
-      const auto& radial_distribution_this_layer =
-          radial_distribution.at(layer_i);
-      if (layer_i != radial_partitioning.size()) {
-        temp_outer_radius = radial_partitioning.at(layer_i);
-      } else {
-        temp_outer_radius = outer_radius;
-      }
-      // Generate wedges/half-wedges a layer at a time.
-      std::vector<Wedge3DMap> wedges_for_this_layer{};
-      if (not use_half_wedges) {
-        for (size_t face_j = which_wedge_index(which_wedges); face_j < 6;
-             face_j++) {
-          wedges_for_this_layer.emplace_back(
-              temp_inner_radius, temp_outer_radius, temp_inner_sphericity,
-              outer_sphericity, gsl::at(wedge_orientations, face_j),
-              use_equiangular_map, Halves::Both, radial_distribution_this_layer,
-              std::array<double, 2>({{M_PI_2, M_PI_2}}));
-        }
-      } else {
-        for (size_t i = 0; i < 4; i++) {
-          wedges_for_this_layer.emplace_back(
-              temp_inner_radius, temp_outer_radius, temp_inner_sphericity,
-              outer_sphericity, gsl::at(wedge_orientations, i),
-              use_equiangular_map, Halves::LowerOnly,
-              radial_distribution_this_layer,
-              std::array<double, 2>({{opening_angle, M_PI_2}}));
-          wedges_for_this_layer.emplace_back(
-              temp_inner_radius, temp_outer_radius, temp_inner_sphericity,
-              outer_sphericity, gsl::at(wedge_orientations, i),
-              use_equiangular_map, Halves::UpperOnly,
-              radial_distribution_this_layer,
-              std::array<double, 2>({{opening_angle, M_PI_2}}));
-        }
-        const double endcap_opening_angle = M_PI - opening_angle;
-        const std::array<double, 2> endcap_opening_angles = {
-            {endcap_opening_angle, endcap_opening_angle}};
-        wedges_for_this_layer.emplace_back(
-            temp_inner_radius, temp_outer_radius, temp_inner_sphericity,
-            outer_sphericity, gsl::at(wedge_orientations, 4),
-            use_equiangular_map, Halves::Both, radial_distribution_this_layer,
-            endcap_opening_angles, false);
-        wedges_for_this_layer.emplace_back(
-            temp_inner_radius, temp_outer_radius, temp_inner_sphericity,
-            outer_sphericity, gsl::at(wedge_orientations, 5),
-            use_equiangular_map, Halves::Both, radial_distribution_this_layer,
-            endcap_opening_angles, false);
-      }
-      for (const auto& wedge : wedges_for_this_layer) {
-        wedges_for_all_layers.push_back(wedge);
-      }
-
-      if (layer_i != radial_partitioning.size()) {
-        temp_inner_radius = radial_partitioning.at(layer_i);
-        temp_inner_sphericity = outer_sphericity;
-      }
-    }
-  }
-  return wedges_for_all_layers;
-}
-
 std::vector<domain::CoordinateMaps::Frustum> frustum_coordinate_maps(
     const double length_inner_cube, const double length_outer_cube,
     const bool equiangular_map_at_outer, const bool equiangular_map_at_inner,
diff --git a/src/Domain/DomainHelpers.hpp b/src/Domain/DomainHelpers.hpp
index 6c9c802d28f9..40d80befd80d 100644
--- a/src/Domain/DomainHelpers.hpp
+++ b/src/Domain/DomainHelpers.hpp
@@ -16,6 +16,7 @@
 #include "DataStructures/Index.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Domain/CoordinateMaps/Distribution.hpp"
+#include "Domain/CoordinateMaps/ShellType.hpp"
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/Side.hpp"
 #include "Utilities/ConstantExpressions.hpp"
@@ -179,7 +180,10 @@ size_t which_wedge_index(const ShellWedges& which_wedges);
  * of pi minus this opening angle. This parameter only has an effect if
  * `use_half_wedges` is set to `true`.
  */
-std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
+template <typename SourceFrame, typename TargetFrame, typename... AppendMaps>
+std::vector<
+    std::unique_ptr<domain::CoordinateMapBase<SourceFrame, TargetFrame, 3>>>
+spherical_shells_coordinate_maps(
     double inner_radius, double outer_radius, double inner_sphericity,
     double outer_sphericity, bool use_equiangular_map,
     const std::optional<std::pair<double, std::array<double, 3>>>&
@@ -188,7 +192,10 @@ std::vector<domain::CoordinateMaps::Wedge<3>> sph_wedge_coordinate_maps(
     const std::vector<double>& radial_partitioning = {},
     const std::vector<domain::CoordinateMaps::Distribution>&
         radial_distribution = {domain::CoordinateMaps::Distribution::Linear},
-    ShellWedges which_wedges = ShellWedges::All, double opening_angle = M_PI_2);
+    const std::vector<domain::CoordinateMaps::ShellType>& shell_types =
+        {domain::CoordinateMaps::ShellType::Cubed},
+    ShellWedges which_wedges = ShellWedges::All, double opening_angle = M_PI_2,
+    const AppendMaps&... append_maps);
 
 /// \ingroup ComputationalDomainGroup
 /// These are the ten Frustums used in the DomainCreators for binary compact
diff --git a/src/Domain/DomainHelpers.tpp b/src/Domain/DomainHelpers.tpp
new file mode 100644
index 000000000000..8d950ee301a2
--- /dev/null
+++ b/src/Domain/DomainHelpers.tpp
@@ -0,0 +1,140 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include "Domain/DomainHelpers.hpp"
+
+#include <cstddef>
+#include <memory>
+#include <vector>
+
+#include "Domain/CoordinateMaps/CoordinateMap.tpp"
+#include "Domain/CoordinateMaps/Distribution.hpp"
+#include "Domain/CoordinateMaps/Identity.hpp"
+#include "Domain/CoordinateMaps/Interval.hpp"
+#include "Domain/CoordinateMaps/ProductMaps.tpp"
+#include "Domain/CoordinateMaps/ShellType.hpp"
+#include "Domain/CoordinateMaps/Wedge.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Domain/CoordinateMaps/SphericalToCartesian.hpp"
+#include "Domain/CoordinateMaps/Affine.hpp"
+
+template <typename SourceFrame, typename TargetFrame, typename... AppendMaps>
+std::vector<
+    std::unique_ptr<domain::CoordinateMapBase<SourceFrame, TargetFrame, 3>>>
+spherical_shells_coordinate_maps(
+    const double inner_radius, const double outer_radius,
+    const double inner_sphericity, const double outer_sphericity,
+    const bool use_equiangular_map, const bool use_half_wedges,
+    const std::vector<double>& radial_partitioning,
+    const std::vector<domain::CoordinateMaps::Distribution>&
+        radial_distribution,
+    const std::vector<domain::CoordinateMaps::ShellType>& shell_types,
+    const ShellWedges which_wedges, const double opening_angle,
+    const AppendMaps&... append_maps) {
+  ASSERT(not use_half_wedges or which_wedges == ShellWedges::All,
+         "If we are using half wedges we must also be using ShellWedges::All.");
+  ASSERT(radial_distribution.size() == 1 + radial_partitioning.size(),
+         "Specify a radial distribution for every spherical shell. You "
+         "specified "
+             << radial_distribution.size() << " items, but the domain has "
+             << 1 + radial_partitioning.size() << " shells.");
+  ASSERT(shell_types.size() == 1 + radial_partitioning.size(),
+         "Specify a type for every spherical shell. You specified "
+             << shell_types.size() << " items, but the domain has "
+             << 1 + radial_partitioning.size() << " shells.");
+
+  const auto wedge_orientations = orientations_for_sphere_wrappings();
+
+  std::vector<
+      std::unique_ptr<domain::CoordinateMapBase<SourceFrame, TargetFrame, 3>>>
+      maps{};
+
+  const size_t number_of_layers = 1 + radial_partitioning.size();
+  double temp_inner_radius = inner_radius;
+  double temp_outer_radius{};
+  double temp_inner_sphericity = inner_sphericity;
+  for (size_t layer_i = 0; layer_i < number_of_layers; layer_i++) {
+    // Determine inner and outer radius for this shell
+    if (layer_i != radial_partitioning.size()) {
+      temp_outer_radius = radial_partitioning.at(layer_i);
+    } else {
+      temp_outer_radius = outer_radius;
+    }
+    if (shell_types[layer_i] == domain::CoordinateMaps::ShellType::Cubed) {
+      // Compose the shell of wedges (deformed cubes)
+      using Wedge3DMap = domain::CoordinateMaps::Wedge<3>;
+      using Halves = Wedge3DMap::WedgeHalves;
+      std::vector<Wedge3DMap> wedges{};
+      if (not use_half_wedges) {
+        for (size_t face_j = which_wedge_index(which_wedges); face_j < 6;
+             face_j++) {
+          wedges.emplace_back(
+              temp_inner_radius, temp_outer_radius, temp_inner_sphericity,
+              outer_sphericity, gsl::at(wedge_orientations, face_j),
+              use_equiangular_map, Halves::Both, radial_distribution[layer_i],
+              std::array<double, 2>({{M_PI_2, M_PI_2}}));
+        }
+      } else {
+        for (size_t i = 0; i < 4; i++) {
+          wedges.emplace_back(temp_inner_radius, temp_outer_radius,
+                              temp_inner_sphericity, outer_sphericity,
+                              gsl::at(wedge_orientations, i),
+                              use_equiangular_map, Halves::LowerOnly,
+                              radial_distribution[layer_i],
+                              std::array<double, 2>({{opening_angle, M_PI_2}}));
+          wedges.emplace_back(temp_inner_radius, temp_outer_radius,
+                              temp_inner_sphericity, outer_sphericity,
+                              gsl::at(wedge_orientations, i),
+                              use_equiangular_map, Halves::UpperOnly,
+                              radial_distribution[layer_i],
+                              std::array<double, 2>({{opening_angle, M_PI_2}}));
+        }
+        const double endcap_opening_angle = M_PI - opening_angle;
+        const std::array<double, 2> endcap_opening_angles = {
+            {endcap_opening_angle, endcap_opening_angle}};
+        wedges.emplace_back(temp_inner_radius, temp_outer_radius,
+                            temp_inner_sphericity, outer_sphericity,
+                            gsl::at(wedge_orientations, 4), use_equiangular_map,
+                            Halves::Both, radial_distribution[layer_i],
+                            endcap_opening_angles, false);
+        wedges.emplace_back(temp_inner_radius, temp_outer_radius,
+                            temp_inner_sphericity, outer_sphericity,
+                            gsl::at(wedge_orientations, 5), use_equiangular_map,
+                            Halves::Both, radial_distribution[layer_i],
+                            endcap_opening_angles, false);
+      }
+      auto wedge_maps =
+          domain::make_vector_coordinate_map_base<SourceFrame, TargetFrame, 3>(
+              std::move(wedges), append_maps...);
+      maps.insert(maps.end(), std::make_move_iterator(wedge_maps.begin()),
+                  std::make_move_iterator(wedge_maps.end()));
+    } else if (shell_types[layer_i] ==
+               domain::CoordinateMaps::ShellType::Spherical) {
+      // Use a single spherical shell (with spherical harmonics in angular
+      // directions)
+      ASSERT(inner_sphericity == 1. and outer_sphericity == 1.,
+             "Spherical shells only support inner and outer sphericity == 1.");
+      domain::CoordinateMaps::ProductOf3Maps<
+          domain::CoordinateMaps::Identity<1>,
+          domain::CoordinateMaps::Identity<1>, domain::CoordinateMaps::Interval>
+          logical_to_spherical_map{{},
+                                   {},
+                                   {-1., 1., temp_inner_radius, outer_radius,
+                                    radial_distribution[layer_i], 0.}};
+      maps.push_back(domain::make_coordinate_map_base<SourceFrame, TargetFrame>(
+          std::move(logical_to_spherical_map),
+          domain::CoordinateMaps::SphericalToCartesian<3>{}, append_maps...));
+    } else {
+      ERROR("Unknown ShellType");
+    }
+
+    if (layer_i != radial_partitioning.size()) {
+      temp_inner_radius = radial_partitioning.at(layer_i);
+      temp_inner_sphericity = outer_sphericity;
+    }
+  }
+  return maps;
+}
diff --git a/src/Executables/ExportCoordinates/ExportCoordinates.hpp b/src/Executables/ExportCoordinates/ExportCoordinates.hpp
index d7853a1e5aae..d2e8f3e3aca3 100644
--- a/src/Executables/ExportCoordinates/ExportCoordinates.hpp
+++ b/src/Executables/ExportCoordinates/ExportCoordinates.hpp
@@ -130,13 +130,13 @@ struct ExportCoordinates {
     const double time = get<Tags::Time>(box);
 
     const auto& mesh = get<domain::Tags::Mesh<Dim>>(box);
-    const auto& inv_jacobian =
-        db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
-                                              Frame::Inertial>>(box);
+    // const auto& inv_jacobian =
+    //     db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
+    //                                           Frame::Inertial>>(box);
     const auto& inertial_coordinates =
         db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
-    const auto deriv_inertial_coordinates =
-        partial_derivative(inertial_coordinates, mesh, inv_jacobian);
+    // const auto deriv_inertial_coordinates =
+    //     partial_derivative(inertial_coordinates, mesh, inv_jacobian);
     // Collect volume data
     // Remove tensor types, only storing individual components
     std::vector<TensorComponent> components;
@@ -148,13 +148,13 @@ struct ExportCoordinates {
                               inertial_coordinates.get(d));
     }
 
-    for (size_t i = 0; i < deriv_inertial_coordinates.size(); ++i) {
-      components.emplace_back(
-          "DerivInertialCoordinates_" +
-              deriv_inertial_coordinates.component_name(
-                  deriv_inertial_coordinates.get_tensor_index(i)),
-          deriv_inertial_coordinates[i]);
-    }
+    // for (size_t i = 0; i < deriv_inertial_coordinates.size(); ++i) {
+    //   components.emplace_back(
+    //       "DerivInertialCoordinates_" +
+    //           deriv_inertial_coordinates.component_name(
+    //               deriv_inertial_coordinates.get_tensor_index(i)),
+    //       deriv_inertial_coordinates[i]);
+    // }
 
     // Also output the determinant of the inverse jacobian, which measures
     // the expansion and compression of the grid
@@ -169,27 +169,27 @@ struct ExportCoordinates {
     // Also output the jacobian diagnostic, which compares the analytic
     // Jacobian (via the CoordinateMap) to the numerical Jacobian
     // (computed via logical_partial_derivative)
-    const auto& jacobian = determinant_and_inverse(inv_jacobian).second;
-    tnsr::i<DataVector, Dim, Frame::ElementLogical> jac_diag{
-        mesh.number_of_grid_points(), 0.0};
-    domain::jacobian_diagnostic(make_not_null(&jac_diag), jacobian,
-                                inertial_coordinates, mesh);
-    for (size_t i = 0; i < Dim; ++i) {
-      components.emplace_back(
-          "JacobianDiagnostic_" +
-              jac_diag.component_name(jac_diag.get_tensor_index(i)),
-          jac_diag.get(i));
-    }
+    // const auto& jacobian = determinant_and_inverse(inv_jacobian).second;
+    // tnsr::i<DataVector, Dim, Frame::ElementLogical> jac_diag{
+    //     mesh.number_of_grid_points(), 0.0};
+    // domain::jacobian_diagnostic(make_not_null(&jac_diag), jacobian,
+    //                             inertial_coordinates, mesh);
+    // for (size_t i = 0; i < Dim; ++i) {
+    //   components.emplace_back(
+    //       "JacobianDiagnostic_" +
+    //           jac_diag.component_name(jac_diag.get_tensor_index(i)),
+    //       jac_diag.get(i));
+    // }
 
     // Also output the computation domain metric
-    const auto& flat_logical_metric =
-        db::get<domain::Tags::FlatLogicalMetric<Dim>>(box);
-    for (size_t i = 0; i < flat_logical_metric.size(); ++i) {
-      components.emplace_back(
-          db::tag_name<domain::Tags::FlatLogicalMetric<Dim>>() +
-              flat_logical_metric.component_suffix(i),
-          flat_logical_metric[i]);
-    }
+    // const auto& flat_logical_metric =
+    //     db::get<domain::Tags::FlatLogicalMetric<Dim>>(box);
+    // for (size_t i = 0; i < flat_logical_metric.size(); ++i) {
+    //   components.emplace_back(
+    //       db::tag_name<domain::Tags::FlatLogicalMetric<Dim>>() +
+    //           flat_logical_metric.component_suffix(i),
+    //       flat_logical_metric[i]);
+    // }
 
     // Send data to volume observer
     auto& local_observer = *Parallel::local_branch(
diff --git a/tests/InputFiles/ExportCoordinates/Input3D.yaml b/tests/InputFiles/ExportCoordinates/Input3D.yaml
index 919a4a658c49..0d83fa565aad 100644
--- a/tests/InputFiles/ExportCoordinates/Input3D.yaml
+++ b/tests/InputFiles/ExportCoordinates/Input3D.yaml
@@ -30,17 +30,19 @@ ResourceInfo:
   Singletons: Auto
 
 DomainCreator:
-  Brick:
-    LowerBound: [-0.5, -0.5, -0.5]
-    UpperBound: [0.5, 0.5, 0.5]
-    Distribution: [Linear, Linear, Linear]
-    IsPeriodicIn: [false, false, false]
-    InitialRefinement: [0, 0, 0]
-    InitialGridPoints: [3, 3, 3]
-    TimeDependence:
-      UniformTranslation:
-        InitialTime: 0.0
-        Velocity: [0.5, 0.0, 0.0]
+  Sphere:
+    InnerRadius: 2.
+    OuterRadius: 10.
+    Interior: Excise
+    InitialRefinement: 0
+    InitialGridPoints: 12
+    UseEquiangularMap: True
+    EquatorialCompression: None
+    WhichWedges: All
+    RadialPartitioning: [4.]
+    RadialDistribution: [Linear, Logarithmic]
+    ShellType: [Cubed, Spherical]
+    TimeDependentMaps: None
 
 SpatialDiscretization:
   ActiveGrid: Dg