Create Ylm meshes in spherical shells

sxs-collaboration · Dec 14, 2024 · e4da1bd · e4da1bd
1 parent b89743b
commit e4da1bd
Show file tree

Hide file tree

Showing 8 changed files with 40 additions and 21 deletions.
diff --git a/src/Domain/ElementDistribution.cpp b/src/Domain/ElementDistribution.cpp
@@ -63,7 +63,7 @@ double get_num_points_and_grid_spacing_cost(
     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);
+      initial_extents, element_id, block.geometry(), quadrature);
   Element<Dim> element = ::domain::Initialization::create_initial_element(
       element_id, blocks, initial_refinement_levels);
   ElementMap<Dim, Frame::Grid> element_map{element_id, block};

diff --git a/src/Domain/Structure/CreateInitialMesh.cpp b/src/Domain/Structure/CreateInitialMesh.cpp
@@ -20,25 +20,38 @@ namespace domain::Initialization {
 template <size_t Dim>
 Mesh<Dim> create_initial_mesh(
     const std::vector<std::array<size_t, Dim>>& initial_extents,
-    const ElementId<Dim>& element_id, const Spectral::Quadrature quadrature,
+    const ElementId<Dim>& element_id, const domain::BlockGeometry geometry,
+    const Spectral::Quadrature quadrature,
     const OrientationMap<Dim>& orientation) {
   const auto& unoriented_extents = initial_extents[element_id.block_id()];
   Index<Dim> extents;
   for (size_t i = 0; i < Dim; ++i) {
     extents[i] = gsl::at(unoriented_extents, orientation(i));
   }
-  return {extents.indices(), Spectral::Basis::Legendre, quadrature};
+  if (geometry == domain::BlockGeometry::Cube) {
+    return {extents.indices(), Spectral::Basis::Legendre, quadrature};
+  } else {
+    if constexpr (Dim == 3) {
+      return {extents.indices(),
+              {{Spectral::Basis::Legendre, Spectral::Basis::SphericalHarmonic,
+                Spectral::Basis::SphericalHarmonic}},
+              {{quadrature, Spectral::Quadrature::Gauss,
+                Spectral::Quadrature::Equiangular}}};
+    } else {
+      ERROR("not implemented");
+    }
+  }
 }
 }  // namespace domain::Initialization
 
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
-#define INSTANTIATE(_, data)                                              \
-  template Mesh<DIM(data)>                                                \
-  domain::Initialization::create_initial_mesh<DIM(data)>(                 \
-      const std::vector<std::array<size_t, DIM(data)>>&,                  \
-      const ElementId<DIM(data)>&, const Spectral::Quadrature quadrature, \
-      const OrientationMap<DIM(data)>&);
+#define INSTANTIATE(_, data)                                       \
+  template Mesh<DIM(data)>                                         \
+  domain::Initialization::create_initial_mesh<DIM(data)>(          \
+      const std::vector<std::array<size_t, DIM(data)>>&,           \
+      const ElementId<DIM(data)>&, domain::BlockGeometry geometry, \
+      Spectral::Quadrature quadrature, const OrientationMap<DIM(data)>&);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
 

diff --git a/src/Domain/Structure/CreateInitialMesh.hpp b/src/Domain/Structure/CreateInitialMesh.hpp
@@ -8,12 +8,14 @@
 #include <cstdint>
 #include <vector>
 
+#include "Domain/Structure/BlockGeometry.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "NumericalAlgorithms/Spectral/Quadrature.hpp"
+
 /// \cond
 template <size_t Dim>
 struct ElementId;
 template <size_t Dim>
-class Mesh;
-template <size_t Dim>
 struct OrientationMap;
 namespace Spectral {
 enum class Quadrature : uint8_t;
@@ -36,7 +38,9 @@ namespace domain::Initialization {
 template <size_t Dim>
 Mesh<Dim> create_initial_mesh(
     const std::vector<std::array<size_t, Dim>>& initial_extents,
-    const ElementId<Dim>& element_id, Spectral::Quadrature quadrature,
+    const ElementId<Dim>& element_id,
+    domain::BlockGeometry geometry = domain::BlockGeometry::Cube,
+    Spectral::Quadrature quadrature = Spectral::Quadrature::GaussLobatto,
     const OrientationMap<Dim>& orientation =
         OrientationMap<Dim>::create_aligned());
 }  // namespace domain::Initialization
diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp
@@ -97,14 +97,15 @@ void InitializeGeometry<Dim>::apply(
     const Domain<Dim>& domain,
     const domain::FunctionsOfTimeMap& functions_of_time,
     const Spectral::Quadrature quadrature, const ElementId<Dim>& element_id) {
+  const auto& block = domain.blocks()[element_id.block_id()];
   // Mesh
   ASSERT(quadrature == Spectral::Quadrature::GaussLobatto or
              quadrature == Spectral::Quadrature::Gauss,
          "The elliptic DG scheme supports Gauss and Gauss-Lobatto "
          "grids, but the chosen quadrature is: "
              << quadrature);
-  *mesh = domain::Initialization::create_initial_mesh(initial_extents,
-                                                      element_id, quadrature);
+  *mesh = domain::Initialization::create_initial_mesh(
+      initial_extents, element_id, block.geometry(), quadrature);
   // Element
   *element = domain::Initialization::create_initial_element(
       element_id, domain.blocks(), initial_refinement);
@@ -117,7 +118,6 @@ void InitializeGeometry<Dim>::apply(
     }
   }
   // Element map
-  const auto& block = domain.blocks()[element_id.block_id()];
   *element_map = ElementMap<Dim, Frame::Inertial>{element_id, block};
   // Coordinates and Jacobians
   detail::initialize_coords_and_jacobians(

diff --git a/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp b/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp
@@ -56,8 +56,8 @@ mortars_apply_impl(const std::vector<std::array<size_t, Dim>>& initial_extents,
           mortar_id,
           ::dg::mortar_mesh(volume_mesh.slice_away(direction.dimension()),
                             ::domain::Initialization::create_initial_mesh(
-                                initial_extents, neighbor, quadrature,
-                                neighbors.orientation())
+                                initial_extents, neighbor, neighbors.geometry(),
+                                quadrature, neighbors.orientation())
                                 .slice_away(direction.dimension())));
       mortar_sizes.emplace(
           mortar_id,

diff --git a/src/Evolution/Initialization/DgDomain.hpp b/src/Evolution/Initialization/DgDomain.hpp
@@ -145,7 +145,7 @@ struct Domain {
       const ElementId<Dim>& element_id) {
     const auto& my_block = domain.blocks()[element_id.block_id()];
     *mesh = ::domain::Initialization::create_initial_mesh(
-        initial_extents, element_id, quadrature);
+        initial_extents, element_id, my_block.geometry(), quadrature);
     *element = ::domain::Initialization::create_initial_element(
         element_id, domain.blocks(), initial_refinement);
     *element_map = ElementMap<Dim, Frame::Grid>{element_id, my_block};

diff --git a/...ems/CurvedScalarWave/Worldtube/SingletonActions/InitializeElementFacesGridCoordinates.hpp b/...ems/CurvedScalarWave/Worldtube/SingletonActions/InitializeElementFacesGridCoordinates.hpp
@@ -82,10 +82,11 @@ struct InitializeElementFacesGridCoordinates {
       for (const auto element_id : element_ids) {
         const auto direction = excision_sphere.abutting_direction(element_id);
         if (direction.has_value()) {
+          const auto& current_block = blocks.at(block_id);
           const auto mesh = ::domain::Initialization::create_initial_mesh(
-              initial_extents, element_id, quadrature);
+              initial_extents, element_id, current_block.geometry(),
+              quadrature);
           const auto face_mesh = mesh.slice_away(direction.value().dimension());
-          const auto& current_block = blocks.at(block_id);
           const ElementMap<Dim, Frame::Grid> element_map{
               element_id,
               current_block.is_time_dependent()

diff --git a/src/Executables/Examples/RandomAmr/InitializeDomain.hpp b/src/Executables/Examples/RandomAmr/InitializeDomain.hpp
@@ -63,8 +63,9 @@ struct Domain {
       const std::vector<std::array<size_t, Dim>>& initial_refinement,
       const Spectral::Quadrature& quadrature,
       const ElementId<Dim>& element_id) {
+    const auto& block = domain.blocks()[element_id.block_id()];
     *mesh = ::domain::Initialization::create_initial_mesh(
-        initial_extents, element_id, quadrature);
+        initial_extents, element_id, block.geometry(), quadrature);
     *element = ::domain::Initialization::create_initial_element(
         element_id, domain.blocks(), initial_refinement);
   }