Add Topology::B1Radial

Author: michaeldmurphy1

src/Domain/BlockLogicalCoordinates.cpp

@@ -114,6 +114,7 @@ block_logical_coordinates_single_point(
   for (size_t d = 0; d < Dim; ++d) {
     const auto topology = gsl::at(block.topologies(), d);
     if (topology == domain::Topology::I1 or
+        topology == domain::Topology::B1Radial or
         topology == domain::Topology::B2Radial or
         topology == domain::Topology::B3Radial) {
       // Map inverses may report logical coordinates outside [-1, 1] due to

src/Domain/Structure/CreateInitialMesh.cpp

@@ -11,6 +11,7 @@
 #include "Domain/Structure/Element.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Structure/OrientationMap.hpp"
+#include "Domain/Structure/Topology.hpp"
 #include "NumericalAlgorithms/Spectral/Basis.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
@@ -52,6 +53,8 @@ std::array<Spectral::Basis, Dim> make_basis(
         [[fallthrough]];
       case (domain::Topology::CartoonCylinder):
         return Spectral::Basis::Cartoon;
+      case (domain::Topology::B1Radial):
+        return Spectral::Basis::ZernikeB1;
       default:
         ERROR("Invalid topology");
     }
@@ -93,6 +96,8 @@ std::array<Spectral::Quadrature, Dim> make_quadrature(
             return Spectral::Quadrature::SphericalSymmetry;
           case (domain::Topology::CartoonCylinder):
             return Spectral::Quadrature::AxialSymmetry;
+          case (domain::Topology::B1Radial):
+            return Spectral::Quadrature::GaussRadauUpper;
           default:
             ERROR("Invalid topology");
         }

src/Domain/Structure/HasBoundary.cpp

@@ -8,7 +8,7 @@
 
 namespace domain {
 bool has_boundary(const domain::Topology topology, const Side side) {
-  return topology == Topology::I1 or
+  return topology == Topology::I1 or topology == Topology::B1Radial or
          (side == Side::Upper and
           (topology == Topology::B2Radial or topology == Topology::B3Radial));
 }

src/Domain/Structure/NeighborIsConforming.cpp

@@ -69,9 +69,28 @@ bool neighbor_is_conforming(
           return false;
         }
       }
+      // permute_from_neighbor reorders by dimension index only, never flipping
+      // sides. So the correct direction into the permuted array has the
+      // dimension from searching which self index maps to the neighbor shared
+      // face dimension, and the side directly from the orientation-mapped
+      // shared face direction.
+      const auto neighbor_shared_face =
+          orientation_of_neighbor(direction_to_neighbor.opposite());
+      std::optional<size_t> permuted_dim{};
+      for (size_t i = 0; i < VolumeDim; ++i) {
+        if (orientation_of_neighbor(i) == neighbor_shared_face.dimension()) {
+          permuted_dim = i;
+          break;
+        }
+      }
+      ASSERT(permuted_dim.has_value(),
+             "No dimension found matching neighbor shared face dimension "
+                 << neighbor_shared_face.dimension());
+      const Direction<VolumeDim> permuted_neighbor_direction{
+          permuted_dim.value(), neighbor_shared_face.side()};
       const auto neighbor_boundary_topologies = boundary_topologies(
           orientation_of_neighbor.permute_from_neighbor(neighbor_topologies),
-          direction_to_neighbor);
+          permuted_neighbor_direction);
       return self_boundary_topologies == neighbor_boundary_topologies;
     }
   }

src/Domain/Structure/Topology.cpp

@@ -20,6 +20,8 @@ std::ostream& operator<<(std::ostream& os, const Topology topology) {
       return os << "S2Colatitude";
     case Topology::S2Longitude:
       return os << "S2Longitude";
+    case Topology::B1Radial:
+      return os << "B1Radial";
     case Topology::B2Radial:
       return os << "B2Radial";
     case Topology::B2Angular:

src/Domain/Structure/Topology.hpp

@@ -47,13 +47,14 @@ enum class Topology : uint8_t {
   S1 = 2,
   S2Colatitude = 3,
   S2Longitude = 4,
-  B2Radial = 5,
-  B2Angular = 6,
-  B3Radial = 7,
-  B3Colatitude = 8,
-  B3Longitude = 9,
-  CartoonSphere = 10,
-  CartoonCylinder = 11
+  B1Radial = 5,
+  B2Radial = 6,
+  B2Angular = 7,
+  B3Radial = 8,
+  B3Colatitude = 9,
+  B3Longitude = 10,
+  CartoonSphere = 11,
+  CartoonCylinder = 12
 };
 
 /// Output operator for a Topology.
@@ -83,8 +84,14 @@ static constexpr auto full_sphere = std::array{
 static constexpr auto cartoon_sphere =
     std::array{Topology::I1, Topology::CartoonSphere, Topology::CartoonSphere};
 
+static constexpr auto cartoon_sphere_inner = std::array{
+    Topology::B1Radial, Topology::CartoonSphere, Topology::CartoonSphere};
+
 static constexpr auto cartoon_cylinder =
     std::array{Topology::I1, Topology::I1, Topology::CartoonCylinder};
+
+static constexpr auto cartoon_cylinder_inner =
+    std::array{Topology::B1Radial, Topology::I1, Topology::CartoonCylinder};
 }  // namespace topologies
 
 }  // namespace domain

tests/Unit/Domain/Structure/Test_HasBoundary.cpp

@@ -11,6 +11,7 @@ namespace {
 void test() {
   for (const auto side : std::array{Side::Lower, Side::Upper}) {
     CHECK(domain::has_boundary(domain::Topology::I1, side));
+    CHECK(domain::has_boundary(domain::Topology::B1Radial, side));
     CHECK_FALSE(domain::has_boundary(domain::Topology::S1, side));
     CHECK_FALSE(domain::has_boundary(domain::Topology::S2Colatitude, side));
     CHECK_FALSE(domain::has_boundary(domain::Topology::S2Longitude, side));

tests/Unit/Domain/Structure/Test_NeighborIsConforming.cpp

@@ -224,6 +224,35 @@ void test_3d() {
   CHECK(neighbor_is_conforming(domain::topologies::cylindrical_shell,
                                domain::topologies::full_cylinder,
                                Direction<3>::upper_zeta(), shell_to_full));
+
+  // Configurations of blocks in CartoonSphere2D domain
+  const OrientationMap<3> cartoonsphere2d_map(std::array<Direction<3>, 3>{
+      Direction<3>::lower_eta(), Direction<3>::upper_xi(),
+      Direction<3>::upper_zeta()});
+  const OrientationMap<3> cartoonsphere2d_map2(std::array<Direction<3>, 3>{
+      Direction<3>::upper_eta(), Direction<3>::lower_xi(),
+      Direction<3>::upper_zeta()});
+  const OrientationMap<3> cartoonsphere2d_map3(std::array<Direction<3>, 3>{
+      Direction<3>::lower_xi(), Direction<3>::lower_eta(),
+      Direction<3>::upper_zeta()});
+  CHECK(neighbor_is_conforming(domain::topologies::cartoon_cylinder_inner,
+                               domain::topologies::cartoon_cylinder,
+                               Direction<3>::upper_xi(), aligned));
+  CHECK(neighbor_is_conforming(domain::topologies::cartoon_cylinder,
+                               domain::topologies::cartoon_cylinder_inner,
+                               Direction<3>::lower_xi(), aligned));
+  CHECK(neighbor_is_conforming(domain::topologies::cartoon_cylinder_inner,
+                               domain::topologies::cartoon_cylinder,
+                               Direction<3>::upper_xi(), cartoonsphere2d_map2));
+  CHECK(neighbor_is_conforming(domain::topologies::cartoon_cylinder,
+                               domain::topologies::cartoon_cylinder_inner,
+                               Direction<3>::lower_eta(), cartoonsphere2d_map));
+  CHECK(neighbor_is_conforming(domain::topologies::cartoon_cylinder_inner,
+                               domain::topologies::cartoon_cylinder,
+                               Direction<3>::upper_xi(), cartoonsphere2d_map3));
+  CHECK(neighbor_is_conforming(domain::topologies::cartoon_cylinder,
+                               domain::topologies::cartoon_cylinder_inner,
+                               Direction<3>::upper_xi(), cartoonsphere2d_map3));
 }
 }  // namespace
 

tests/Unit/Domain/Structure/Test_Topology.cpp

@@ -12,6 +12,7 @@ SPECTRE_TEST_CASE("Unit.Domain.Structure.Topology", "[Domain][Unit]") {
   CHECK(get_output(domain::Topology::S1) == "S1");
   CHECK(get_output(domain::Topology::S2Colatitude) == "S2Colatitude");
   CHECK(get_output(domain::Topology::S2Longitude) == "S2Longitude");
+  CHECK(get_output(domain::Topology::B1Radial) == "B1Radial");
   CHECK(get_output(domain::Topology::B2Radial) == "B2Radial");
   CHECK(get_output(domain::Topology::B2Angular) == "B2Angular");
   CHECK(get_output(domain::Topology::B3Radial) == "B3Radial");