Merge pull request #2602 from nilsdeppe/quadrature_to_create_initial_…

…mesh Add quadrature to create_initial_mesh function
sxs-collaboration · Nov 25, 2020 · ebec864 · ebec864
2 parents fd88f47 + d8e03a3
commit ebec864
Show file tree

Hide file tree

Showing 9 changed files with 59 additions and 55 deletions.
diff --git a/src/Domain/Structure/CreateInitialMesh.cpp b/src/Domain/Structure/CreateInitialMesh.cpp
@@ -19,26 +19,26 @@ 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 ElementId<Dim>& element_id, const Spectral::Quadrature quadrature,
     const OrientationMap<Dim>& orientation) noexcept {
   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,
-          Spectral::Quadrature::GaussLobatto};
+  return {extents.indices(), Spectral::Basis::Legendre, quadrature};
 }
 }  // namespace domain::Initialization
 
 /// \cond
 #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 OrientationMap<DIM(data)>&) noexcept;
+#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)>&) noexcept;
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
 

diff --git a/src/Domain/Structure/CreateInitialMesh.hpp b/src/Domain/Structure/CreateInitialMesh.hpp
@@ -16,8 +16,7 @@ template <size_t Dim>
 struct OrientationMap;
 /// \endcond
 
-namespace domain {
-namespace Initialization {
+namespace domain::Initialization {
 /// \ingroup InitializationGroup
 /// \brief Construct the initial Mesh of an Element.
 ///
@@ -28,11 +27,11 @@ namespace Initialization {
 ///
 /// \param initial_extents the initial extents of each Block in the Domain
 /// \param element_id id of an Element or its neighbor
+/// \param quadrature the quadrature rule/grid point distribution
 /// \param orientation OrientationMap of (neighboring) `element_id`
 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 ElementId<Dim>& element_id, Spectral::Quadrature quadrature,
     const OrientationMap<Dim>& orientation = {}) noexcept;
-}  // namespace Initialization
-}  // namespace domain
+}  // namespace domain::Initialization
diff --git a/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp b/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.cpp
@@ -24,7 +24,8 @@ namespace evolution::dg::Initialization {
 template <size_t Dim, bool AddFluxBoundaryConditionMortars>
 auto Mortars<Dim, AddFluxBoundaryConditionMortars>::apply_impl(
     const std::vector<std::array<size_t, Dim>>& initial_extents,
-    const Element<Dim>& element, const TimeStepId& next_temporal_id,
+    const Spectral::Quadrature quadrature, const Element<Dim>& element,
+    const TimeStepId& next_temporal_id,
     const std::unordered_map<Direction<Dim>, Mesh<Dim - 1>>& interface_meshes,
     const std::unordered_map<Direction<Dim>, Mesh<Dim - 1>>&
         boundary_meshes) noexcept
@@ -41,11 +42,12 @@ auto Mortars<Dim, AddFluxBoundaryConditionMortars>::apply_impl(
       const auto mortar_id = std::make_pair(direction, neighbor);
       mortar_data[mortar_id];  // Default initialize data
       mortar_meshes.emplace(
-          mortar_id, ::dg::mortar_mesh(
-                         interface_meshes.at(direction),
-                         ::domain::Initialization::create_initial_mesh(
-                             initial_extents, neighbor, neighbors.orientation())
-                             .slice_away(direction.dimension())));
+          mortar_id,
+          ::dg::mortar_mesh(interface_meshes.at(direction),
+                            ::domain::Initialization::create_initial_mesh(
+                                initial_extents, neighbor, quadrature,
+                                neighbors.orientation())
+                                .slice_away(direction.dimension())));
       mortar_sizes.emplace(
           mortar_id,
           ::dg::mortar_size(element.id(), neighbor, direction.dimension(),

diff --git a/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp b/src/Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp
@@ -94,6 +94,7 @@ struct Mortars {
       auto [mortar_data, mortar_meshes, mortar_sizes,
             mortar_next_temporal_ids] =
           apply_impl(db::get<::domain::Tags::InitialExtents<Dim>>(box),
+                     Spectral::Quadrature::GaussLobatto,
                      db::get<::domain::Tags::Element<Dim>>(box),
                      db::get<::Tags::TimeStepId>(box),
                      db::get<::domain::Tags::Interface<
@@ -122,7 +123,8 @@ struct Mortars {
                     MortarMap<TimeStepId>>
   apply_impl(
       const std::vector<std::array<size_t, Dim>>& initial_extents,
-      const Element<Dim>& element, const TimeStepId& next_temporal_id,
+      Spectral::Quadrature quadrature, const Element<Dim>& element,
+      const TimeStepId& next_temporal_id,
       const std::unordered_map<Direction<Dim>, Mesh<Dim - 1>>& interface_meshes,
       const std::unordered_map<Direction<Dim>, Mesh<Dim - 1>>&
           boundary_meshes) noexcept;

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

diff --git a/src/ParallelAlgorithms/DiscontinuousGalerkin/InitializeDomain.hpp b/src/ParallelAlgorithms/DiscontinuousGalerkin/InitializeDomain.hpp
@@ -101,7 +101,7 @@ struct InitializeDomain {
     const ElementId<Dim> element_id{array_index};
     const auto& my_block = domain.blocks()[element_id.block_id()];
     Mesh<Dim> mesh = domain::Initialization::create_initial_mesh(
-        initial_extents, element_id);
+        initial_extents, element_id, Spectral::Quadrature::GaussLobatto);
     Element<Dim> element = domain::Initialization::create_initial_element(
         element_id, my_block, initial_refinement);
     if (my_block.is_time_dependent()) {

diff --git a/src/ParallelAlgorithms/DiscontinuousGalerkin/InitializeMortars.hpp b/src/ParallelAlgorithms/DiscontinuousGalerkin/InitializeMortars.hpp
@@ -129,7 +129,8 @@ struct InitializeMortars {
             dg::mortar_mesh(
                 interface_meshes.at(direction),
                 domain::Initialization::create_initial_mesh(
-                    initial_extents, neighbor, neighbors.orientation())
+                    initial_extents, neighbor,
+                    Spectral::Quadrature::GaussLobatto, neighbors.orientation())
                     .slice_away(direction.dimension())));
         mortar_sizes.emplace(
             mortar_id,

diff --git a/tests/Unit/Domain/Structure/Test_CreateInitialMesh.cpp b/tests/Unit/Domain/Structure/Test_CreateInitialMesh.cpp
@@ -17,44 +17,43 @@ namespace domain::Initialization {
 SPECTRE_TEST_CASE("Unit.Domain.Structure.CreateInitialMesh", "[Domain][Unit]") {
   {
     INFO("Single element");
-    CHECK(create_initial_mesh({{{3}}}, ElementId<1>{0}) ==
-          Mesh<1>{3, Spectral::Basis::Legendre,
-                  Spectral::Quadrature::GaussLobatto});
-    CHECK(create_initial_mesh({{{3, 2}}}, ElementId<2>{0}) ==
-          Mesh<2>{{{3, 2}},
-                  Spectral::Basis::Legendre,
-                  Spectral::Quadrature::GaussLobatto});
-    CHECK(create_initial_mesh({{{3, 2, 4}}}, ElementId<3>{0}) ==
-          Mesh<3>{{{3, 2, 4}},
-                  Spectral::Basis::Legendre,
-                  Spectral::Quadrature::GaussLobatto});
+    for (const auto& quadrature :
+         {Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss}) {
+      CHECK(create_initial_mesh({{{3}}}, ElementId<1>{0}, quadrature) ==
+            Mesh<1>{3, Spectral::Basis::Legendre, quadrature});
+      CHECK(create_initial_mesh({{{3, 2}}}, ElementId<2>{0}, quadrature) ==
+            Mesh<2>{{{3, 2}}, Spectral::Basis::Legendre, quadrature});
+      CHECK(create_initial_mesh({{{3, 2, 4}}}, ElementId<3>{0}, quadrature) ==
+            Mesh<3>{{{3, 2, 4}}, Spectral::Basis::Legendre, quadrature});
+    }
   }
   {
     INFO("Another element");
-    CHECK(create_initial_mesh({{{3}}, {{2}}}, ElementId<1>{1}) ==
-          Mesh<1>{2, Spectral::Basis::Legendre,
-                  Spectral::Quadrature::GaussLobatto});
-    CHECK(create_initial_mesh({{{3, 3}}, {{2, 2}}}, ElementId<2>{1}) ==
-          Mesh<2>{2, Spectral::Basis::Legendre,
-                  Spectral::Quadrature::GaussLobatto});
-    CHECK(create_initial_mesh({{{3, 3, 3}}, {{2, 2, 2}}}, ElementId<3>{1}) ==
-          Mesh<3>{2, Spectral::Basis::Legendre,
-                  Spectral::Quadrature::GaussLobatto});
+    for (const auto& quadrature :
+         {Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss}) {
+      CHECK(create_initial_mesh({{{3}}, {{2}}}, ElementId<1>{1}, quadrature) ==
+            Mesh<1>{2, Spectral::Basis::Legendre, quadrature});
+      CHECK(create_initial_mesh({{{3, 3}}, {{2, 2}}}, ElementId<2>{1},
+                                quadrature) ==
+            Mesh<2>{2, Spectral::Basis::Legendre, quadrature});
+      CHECK(create_initial_mesh({{{3, 3, 3}}, {{2, 2, 2}}}, ElementId<3>{1},
+                                quadrature) ==
+            Mesh<3>{2, Spectral::Basis::Legendre, quadrature});
+    }
   }
   {
     INFO("Unaligned orientation");
     OrientationMap<2> unaligned(
         make_array(Direction<2>::lower_eta(), Direction<2>::upper_xi()));
-    CHECK(
-        create_initial_mesh({{{2, 3}}, {{4, 5}}}, ElementId<2>{0}, unaligned) ==
-        Mesh<2>{{{3, 2}},
-                Spectral::Basis::Legendre,
-                Spectral::Quadrature::GaussLobatto});
-    CHECK(
-        create_initial_mesh({{{2, 3}}, {{4, 5}}}, ElementId<2>{1}, unaligned) ==
-        Mesh<2>{{{5, 4}},
-                Spectral::Basis::Legendre,
-                Spectral::Quadrature::GaussLobatto});
+    for (const auto& quadrature :
+         {Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss}) {
+      CHECK(create_initial_mesh({{{2, 3}}, {{4, 5}}}, ElementId<2>{0},
+                                quadrature, unaligned) ==
+            Mesh<2>{{{3, 2}}, Spectral::Basis::Legendre, quadrature});
+      CHECK(create_initial_mesh({{{2, 3}}, {{4, 5}}}, ElementId<2>{1},
+                                quadrature, unaligned) ==
+            Mesh<2>{{{5, 4}}, Spectral::Basis::Legendre, quadrature});
+    }
   }
 }
 

diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Initialization/Test_Mortars.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Initialization/Test_Mortars.cpp
@@ -102,8 +102,9 @@ void test_impl(const std::vector<std::array<size_t, Dim>>& initial_extents,
   ActionTesting::emplace_component_and_initialize<component<metavars>>(
       &runner, element.id(),
       {time_step_id, element,
-       domain::Initialization::create_initial_mesh(initial_extents,
-                                                   element.id(), {})});
+       domain::Initialization::create_initial_mesh(
+           initial_extents, element.id(), Spectral::Quadrature::GaussLobatto,
+           {})});
 
   ActionTesting::set_phase(make_not_null(&runner), metavars::Phase::Testing);