Merge pull request #4675 from kidder/amr_new_mesh

src/ParallelAlgorithms/Amr/CMakeLists.txt

@@ -16,3 +16,4 @@ target_link_libraries(
 
 add_subdirectory(Actions)
 add_subdirectory(Criteria)
+add_subdirectory(Projectors)

src/ParallelAlgorithms/Amr/Projectors/CMakeLists.txt

@@ -0,0 +1,27 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(LIBRARY AmrProjectors)
+
+add_spectre_library(${LIBRARY})
+
+spectre_target_sources(
+  ${LIBRARY}
+  PRIVATE
+  Mesh.cpp
+  )
+
+spectre_target_headers(
+  ${LIBRARY}
+  INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
+  HEADERS
+  Mesh.hpp
+  )
+
+target_link_libraries(
+  ${LIBRARY}
+  PRIVATE
+  Amr
+  DataStructures
+  Spectral
+  )

src/ParallelAlgorithms/Amr/Projectors/Mesh.cpp

@@ -0,0 +1,75 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "ParallelAlgorithms/Amr/Projectors/Mesh.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <vector>
+
+#include "DataStructures/Index.hpp"
+#include "Domain/Amr/Flag.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "Utilities/Algorithm.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/MakeArray.hpp"
+
+namespace amr::projectors {
+template <size_t Dim>
+Mesh<Dim> mesh(const Mesh<Dim>& old_mesh,
+               const std::array<amr::domain::Flag, Dim>& flags) {
+  std::array<size_t, Dim> new_extents = old_mesh.extents().indices();
+  for (size_t d = 0; d < Dim; ++d) {
+    if (gsl::at(flags, d) == amr::domain::Flag::IncreaseResolution) {
+      ++gsl::at(new_extents, d);
+    } else if (gsl::at(flags, d) == amr::domain::Flag::DecreaseResolution) {
+      --gsl::at(new_extents, d);
+    }
+  }
+  return {new_extents, old_mesh.basis(), old_mesh.quadrature()};
+}
+
+template <size_t Dim>
+Mesh<Dim> parent_mesh(const std::vector<Mesh<Dim>>& children_meshes) {
+  const auto& parent_quadrature = children_meshes.front().quadrature();
+  const auto& parent_basis = children_meshes.front().basis();
+
+  ASSERT(alg::all_of(
+             children_meshes,
+             [&parent_quadrature, &parent_basis](const Mesh<Dim>& child_mesh) {
+               return (child_mesh.quadrature() == parent_quadrature and
+                       child_mesh.basis() == parent_basis);
+             }),
+         "AMR does not currently support joining elements with different "
+         "quadratures or bases");
+
+  // loop over each mesh, returning an array containing the max extent in each
+  // dimension
+  auto parent_extents = std::accumulate(
+      std::next(children_meshes.begin()), children_meshes.end(),
+      children_meshes.front().extents().indices(),
+      [](std::array<size_t, Dim>& extents, const Mesh<Dim>& mesh) {
+        alg::transform(extents, mesh.extents().indices(), extents.begin(),
+                       [](size_t a, size_t b) { return std::max(a, b); });
+        return extents;
+      });
+
+  return {parent_extents, parent_basis, parent_quadrature};
+}
+
+#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data)                                  \
+  template Mesh<DIM(data)> mesh(                              \
+      const Mesh<DIM(data)>& old_mesh,                        \
+      const std::array<amr::domain::Flag, DIM(data)>& flags); \
+  template Mesh<DIM(data)> parent_mesh(                       \
+      const std::vector<Mesh<DIM(data)>>& children_meshes);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
+
+#undef DIM
+#undef INSTANTIATE
+}  // namespace amr::projectors

src/ParallelAlgorithms/Amr/Projectors/Mesh.hpp

@@ -0,0 +1,30 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <vector>
+
+#include "Domain/Amr/Flag.hpp"
+
+/// \cond
+template <size_t>
+class Mesh;
+/// \endcond
+
+namespace amr::projectors {
+/// Given `old_mesh` returns a new Mesh based on the refinement `flags`
+template <size_t Dim>
+Mesh<Dim> mesh(const Mesh<Dim>& old_mesh,
+               const std::array<amr::domain::Flag, Dim>& flags);
+
+/// Given the Mesh%es for a set of joining Element%s, returns the Mesh
+/// for the newly created parent Element.
+///
+/// \details In each dimension the extent of the parent Mesh will be the
+/// maximum over the extents of each child Mesh
+template <size_t Dim>
+Mesh<Dim> parent_mesh(const std::vector<Mesh<Dim>>& children_meshes);
+}  // namespace amr::projectors

tests/Unit/ParallelAlgorithms/Amr/CMakeLists.txt

@@ -8,6 +8,7 @@ set(LIBRARY_SOURCES
   Actions/Test_Initialize.cpp
   Actions/Test_UpdateAmrDecision.cpp
   Criteria/Test_Random.cpp
+  Projectors/Test_Mesh.cpp
   )
 
 add_test_library(
@@ -22,6 +23,7 @@ target_link_libraries(
   PRIVATE
   Amr
   AmrCriteria
+  AmrProjectors
   DomainStructure
   Utilities
   )

tests/Unit/ParallelAlgorithms/Amr/Projectors/Test_Mesh.cpp

@@ -0,0 +1,185 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <array>
+#include <cstddef>
+#include <vector>
+
+#include "Domain/Amr/Flag.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "NumericalAlgorithms/Spectral/Spectral.hpp"
+#include "ParallelAlgorithms/Amr/Projectors/Mesh.hpp"
+#include "Utilities/Literals.hpp"
+
+namespace {
+void test_mesh_1d() {
+  const auto legendre = Spectral::Basis::Legendre;
+  const auto gauss_lobatto = Spectral::Quadrature::GaussLobatto;
+  const auto refine = std::array{amr::domain::Flag::IncreaseResolution};
+  const auto coarsen = std::array{amr::domain::Flag::DecreaseResolution};
+  const auto split = std::array{amr::domain::Flag::Split};
+  const auto join = std::array{amr::domain::Flag::Join};
+  const auto stay = std::array{amr::domain::Flag::DoNothing};
+  Mesh<1> mesh_3{std::array{3_st}, legendre, gauss_lobatto};
+  Mesh<1> mesh_4{std::array{4_st}, legendre, gauss_lobatto};
+  CHECK(amr::projectors::mesh(mesh_3, refine) == mesh_4);
+  CHECK(amr::projectors::mesh(mesh_4, coarsen) == mesh_3);
+  CHECK(amr::projectors::mesh(mesh_3, join) == mesh_3);
+  CHECK(amr::projectors::mesh(mesh_3, split) == mesh_3);
+  CHECK(amr::projectors::mesh(mesh_3, stay) == mesh_3);
+
+  CHECK(amr::projectors::parent_mesh(std::vector{mesh_3, mesh_4}) == mesh_4);
+}
+
+void test_mesh_2d() {
+  const auto legendre = Spectral::Basis::Legendre;
+  const auto gauss_lobatto = Spectral::Quadrature::GaussLobatto;
+  const auto refine_refine = std::array{amr::domain::Flag::IncreaseResolution,
+                                        amr::domain::Flag::IncreaseResolution};
+  const auto refine_stay = std::array{amr::domain::Flag::IncreaseResolution,
+                                      amr::domain::Flag::DoNothing};
+  const auto refine_coarsen = std::array{amr::domain::Flag::IncreaseResolution,
+                                         amr::domain::Flag::DecreaseResolution};
+  const auto coarsen_refine = std::array{amr::domain::Flag::DecreaseResolution,
+                                         amr::domain::Flag::IncreaseResolution};
+  const auto coarsen_stay = std::array{amr::domain::Flag::DecreaseResolution,
+                                       amr::domain::Flag::DoNothing};
+  const auto coarsen_coarsen =
+      std::array{amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::DecreaseResolution};
+  Mesh<2> mesh_3_5{std::array{3_st, 5_st}, legendre, gauss_lobatto};
+  Mesh<2> mesh_3_6{std::array{3_st, 6_st}, legendre, gauss_lobatto};
+  Mesh<2> mesh_4_5{std::array{4_st, 5_st}, legendre, gauss_lobatto};
+  Mesh<2> mesh_4_6{std::array{4_st, 6_st}, legendre, gauss_lobatto};
+  CHECK(amr::projectors::mesh(mesh_3_5, refine_refine) == mesh_4_6);
+  CHECK(amr::projectors::mesh(mesh_3_5, refine_stay) == mesh_4_5);
+  CHECK(amr::projectors::mesh(mesh_3_6, refine_coarsen) == mesh_4_5);
+  CHECK(amr::projectors::mesh(mesh_4_5, coarsen_refine) == mesh_3_6);
+  CHECK(amr::projectors::mesh(mesh_4_6, coarsen_stay) == mesh_3_6);
+  CHECK(amr::projectors::mesh(mesh_4_6, coarsen_coarsen) == mesh_3_5);
+
+  CHECK(amr::projectors::parent_mesh(
+            std::vector{mesh_3_5, mesh_4_5, mesh_3_6}) == mesh_4_6);
+#ifdef SPECTRE_DEBUG
+  Mesh<2> mesh_mismatch_basis{std::array{2_st, 4_st},
+                              std::array{legendre, Spectral::Basis::Chebyshev},
+                              std::array{gauss_lobatto, gauss_lobatto}};
+  CHECK_THROWS_WITH(
+      amr::projectors::parent_mesh(std::vector{mesh_3_5, mesh_mismatch_basis}),
+      Catch::Contains("AMR does not currently support joining elements with "
+                      "different quadratures or bases"));
+  Mesh<2> mesh_mismatch_quadrature{
+      std::array{2_st, 4_st}, std::array{legendre, legendre},
+      std::array{Spectral::Quadrature::Gauss, gauss_lobatto}};
+  CHECK_THROWS_WITH(
+      amr::projectors::parent_mesh(
+          std::vector{mesh_3_5, mesh_mismatch_quadrature}),
+      Catch::Contains("AMR does not currently support joining elements with "
+                      "different quadratures or bases"));
+#endif
+}
+
+void test_mesh_3d() {
+  const auto legendre = Spectral::Basis::Legendre;
+  const auto gauss_lobatto = Spectral::Quadrature::GaussLobatto;
+  const auto refine_refine_refine =
+      std::array{amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::IncreaseResolution};
+  const auto refine_stay_refine = std::array{
+      amr::domain::Flag::IncreaseResolution, amr::domain::Flag::DoNothing,
+      amr::domain::Flag::IncreaseResolution};
+  const auto refine_coarsen_refine =
+      std::array{amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::IncreaseResolution};
+  const auto coarsen_refine_refine =
+      std::array{amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::IncreaseResolution};
+  const auto coarsen_stay_refine = std::array{
+      amr::domain::Flag::DecreaseResolution, amr::domain::Flag::DoNothing,
+      amr::domain::Flag::IncreaseResolution};
+  const auto coarsen_coarsen_refine =
+      std::array{amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::IncreaseResolution};
+  const auto refine_refine_coarsen =
+      std::array{amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::DecreaseResolution};
+  const auto refine_stay_coarsen = std::array{
+      amr::domain::Flag::IncreaseResolution, amr::domain::Flag::DoNothing,
+      amr::domain::Flag::DecreaseResolution};
+  const auto refine_coarsen_coarsen =
+      std::array{amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::DecreaseResolution};
+  const auto coarsen_refine_coarsen =
+      std::array{amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::IncreaseResolution,
+                 amr::domain::Flag::DecreaseResolution};
+  const auto coarsen_stay_coarsen = std::array{
+      amr::domain::Flag::DecreaseResolution, amr::domain::Flag::DoNothing,
+      amr::domain::Flag::DecreaseResolution};
+  const auto coarsen_coarsen_coarsen =
+      std::array{amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::DecreaseResolution,
+                 amr::domain::Flag::DecreaseResolution};
+  Mesh<3> mesh_3_5_7{std::array{3_st, 5_st, 7_st}, legendre, gauss_lobatto};
+  Mesh<3> mesh_3_6_7{std::array{3_st, 6_st, 7_st}, legendre, gauss_lobatto};
+  Mesh<3> mesh_4_5_7{std::array{4_st, 5_st, 7_st}, legendre, gauss_lobatto};
+  Mesh<3> mesh_4_6_7{std::array{4_st, 6_st, 7_st}, legendre, gauss_lobatto};
+  Mesh<3> mesh_3_5_8{std::array{3_st, 5_st, 8_st}, legendre, gauss_lobatto};
+  Mesh<3> mesh_3_6_8{std::array{3_st, 6_st, 8_st}, legendre, gauss_lobatto};
+  Mesh<3> mesh_4_5_8{std::array{4_st, 5_st, 8_st}, legendre, gauss_lobatto};
+  Mesh<3> mesh_4_6_8{std::array{4_st, 6_st, 8_st}, legendre, gauss_lobatto};
+  CHECK(amr::projectors::mesh(mesh_3_5_7, refine_refine_refine) == mesh_4_6_8);
+  CHECK(amr::projectors::mesh(mesh_3_5_7, refine_stay_refine) == mesh_4_5_8);
+  CHECK(amr::projectors::mesh(mesh_3_6_7, refine_coarsen_refine) == mesh_4_5_8);
+  CHECK(amr::projectors::mesh(mesh_4_5_7, coarsen_refine_refine) == mesh_3_6_8);
+  CHECK(amr::projectors::mesh(mesh_4_6_7, coarsen_stay_refine) == mesh_3_6_8);
+  CHECK(amr::projectors::mesh(mesh_4_6_7, coarsen_coarsen_refine) ==
+        mesh_3_5_8);
+  CHECK(amr::projectors::mesh(mesh_3_5_8, refine_refine_coarsen) == mesh_4_6_7);
+  CHECK(amr::projectors::mesh(mesh_3_5_8, refine_stay_coarsen) == mesh_4_5_7);
+  CHECK(amr::projectors::mesh(mesh_3_6_8, refine_coarsen_coarsen) ==
+        mesh_4_5_7);
+  CHECK(amr::projectors::mesh(mesh_4_5_8, coarsen_refine_coarsen) ==
+        mesh_3_6_7);
+  CHECK(amr::projectors::mesh(mesh_4_6_8, coarsen_stay_coarsen) == mesh_3_6_7);
+  CHECK(amr::projectors::mesh(mesh_4_6_8, coarsen_coarsen_coarsen) ==
+        mesh_3_5_7);
+
+  CHECK(amr::projectors::parent_mesh(
+            std::vector{mesh_3_5_7, mesh_4_5_8, mesh_3_6_7}) == mesh_4_6_8);
+#ifdef SPECTRE_DEBUG
+  Mesh<3> mesh_mismatch_basis{
+      std::array{2_st, 4_st, 3_st},
+      std::array{legendre, Spectral::Basis::Chebyshev, legendre},
+      std::array{gauss_lobatto, gauss_lobatto, gauss_lobatto}};
+  CHECK_THROWS_WITH(
+      amr::projectors::parent_mesh(
+          std::vector{mesh_3_5_7, mesh_mismatch_basis}),
+      Catch::Contains("AMR does not currently support joining elements with "
+                      "different quadratures or bases"));
+  Mesh<3> mesh_mismatch_quadrature{
+      std::array{2_st, 4_st, 3_st}, std::array{legendre, legendre, legendre},
+      std::array{Spectral::Quadrature::Gauss, gauss_lobatto, gauss_lobatto}};
+  CHECK_THROWS_WITH(
+      amr::projectors::parent_mesh(
+          std::vector{mesh_3_5_7, mesh_mismatch_quadrature}),
+      Catch::Contains("AMR does not currently support joining elements with "
+                      "different quadratures or bases"));
+#endif
+}
+
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.Amr.Projectors.Mesh", "[ParallelAlgorithms][Unit]") {
+  test_mesh_1d();
+  test_mesh_2d();
+  test_mesh_3d();
+}