Merge pull request sxs-collaboration#5643 from kidder/neighbor_mesh

Mutate NeighborMesh at beginning of AMR

src/Domain/Amr/NeighborsOfChild.cpp

@@ -14,23 +14,28 @@
 #include "Domain/Amr/NewNeighborIds.hpp"
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/DirectionMap.hpp"
+#include "Domain/Structure/DirectionalIdMap.hpp"
 #include "Domain/Structure/Element.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Structure/Neighbors.hpp"
 #include "Domain/Structure/OrientationMap.hpp"
 #include "Domain/Structure/SegmentId.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 
 namespace amr {
 template <size_t VolumeDim>
-DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_child(
-    const Element<VolumeDim>& parent,
-    const std::array<Flag, VolumeDim>& parent_flags,
+std::pair<DirectionMap<VolumeDim, Neighbors<VolumeDim>>,
+          DirectionalIdMap<VolumeDim, Mesh<VolumeDim>>>
+neighbors_of_child(
+    const Element<VolumeDim>& parent, const Info<VolumeDim>& parent_info,
     const std::unordered_map<ElementId<VolumeDim>, Info<VolumeDim>>&
         parent_neighbor_info,
     const ElementId<VolumeDim>& child_id) {
-  DirectionMap<VolumeDim, Neighbors<VolumeDim>> result;
+  std::pair<DirectionMap<VolumeDim, Neighbors<VolumeDim>>,
+            DirectionalIdMap<VolumeDim, Mesh<VolumeDim>>>
+      result;
 
   const auto sibling_id = [&child_id](const size_t dim) {
     auto sibling_segment_ids = child_id.segment_ids();
@@ -42,23 +47,33 @@ DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_child(
 
   for (const auto& [direction, old_neighbors] : parent.neighbors()) {
     const auto dim = direction.dimension();
-    if (gsl::at(parent_flags, dim) == Flag::Split and
+    if (gsl::at(parent_info.flags, dim) == Flag::Split and
         has_potential_sibling(child_id, direction)) {
-      result.emplace(direction, Neighbors<VolumeDim>{{sibling_id(dim)}, {}});
+      const auto id = sibling_id(dim);
+      result.first.emplace(direction, Neighbors<VolumeDim>{{id}, {}});
+      result.second.insert({{direction, id}, parent_info.new_mesh});
     } else {
-      result.emplace(direction,
-                     Neighbors<VolumeDim>{
-                         new_neighbor_ids(child_id, direction, old_neighbors,
-                                          parent_neighbor_info),
-                         old_neighbors.orientation()});
+      const auto new_neighbor_ids_and_meshes = amr::new_neighbor_ids(
+          child_id, direction, old_neighbors, parent_neighbor_info);
+      std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids;
+      for (const auto& [id, mesh] : new_neighbor_ids_and_meshes) {
+        result.second.insert_or_assign({direction, id}, mesh);
+        new_neighbor_ids.insert(id);
+      }
+
+      result.first.emplace(
+          direction,
+          Neighbors<VolumeDim>{new_neighbor_ids, old_neighbors.orientation()});
     }
   }
 
   for (const auto& direction : parent.external_boundaries()) {
     const auto dim = direction.dimension();
-    if (gsl::at(parent_flags, dim) == Flag::Split and
+    if (gsl::at(parent_info.flags, dim) == Flag::Split and
         has_potential_sibling(child_id, direction)) {
-      result.emplace(direction, Neighbors<VolumeDim>{{sibling_id(dim)}, {}});
+      const auto id = sibling_id(dim);
+      result.first.emplace(direction, Neighbors<VolumeDim>{{id}, {}});
+      result.second.insert({{direction, id}, parent_info.new_mesh});
     }
   }
 
@@ -67,12 +82,13 @@ DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_child(
 
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
-#define INSTANTIATE(_, data)                                                 \
-  template DirectionMap<DIM(data), Neighbors<DIM(data)>> neighbors_of_child( \
-      const Element<DIM(data)>& parent,                                      \
-      const std::array<Flag, DIM(data)>& parent_flags,                       \
-      const std::unordered_map<ElementId<DIM(data)>, Info<DIM(data)>>&       \
-          parent_neighbor_info,                                              \
+#define INSTANTIATE(_, data)                                                \
+  template std::pair<DirectionMap<DIM(data), Neighbors<DIM(data)>>,         \
+                     DirectionalIdMap<DIM(data), Mesh<DIM(data)>>>          \
+  neighbors_of_child(                                                       \
+      const Element<DIM(data)>& parent, const Info<DIM(data)>& parent_info, \
+      const std::unordered_map<ElementId<DIM(data)>, Info<DIM(data)>>&      \
+          parent_neighbor_info,                                             \
       const ElementId<DIM(data)>& child_id);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))

src/Domain/Amr/NeighborsOfChild.hpp

@@ -15,24 +15,29 @@ template <size_t VolumeDim>
 struct Info;
 }  // namespace amr
 template <size_t VolumeDim, typename T>
+class DirectionalIdMap;
+template <size_t VolumeDim, typename T>
 class DirectionMap;
 template <size_t VolumeDim>
 class Element;
 template <size_t VolumeDim>
 class ElementId;
 template <size_t VolumeDim>
+class Mesh;
+template <size_t VolumeDim>
 class Neighbors;
 /// \endcond
 
 namespace amr {
 /// \ingroup AmrGroup
-/// \brief returns the neighbors of the Element with ElementId `child_id`,
-/// whose parent Element is `parent` which has refinement flags `parent_flags`
+/// \brief returns the neighbors and their Mesh%es of the Element with ElementId
+/// `child_id`, whose parent Element is `parent` which has Info `parent_info`
 /// and neighbor Info `parent_neighbor_info`
 template <size_t VolumeDim>
-DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_child(
-    const Element<VolumeDim>& parent,
-    const std::array<Flag, VolumeDim>& parent_flags,
+std::pair<DirectionMap<VolumeDim, Neighbors<VolumeDim>>,
+          DirectionalIdMap<VolumeDim, Mesh<VolumeDim>>>
+neighbors_of_child(
+    const Element<VolumeDim>& parent, const Info<VolumeDim>& parent_info,
     const std::unordered_map<ElementId<VolumeDim>, Info<VolumeDim>>&
         parent_neighbor_info,
     const ElementId<VolumeDim>& child_id);

src/Domain/Amr/NeighborsOfParent.cpp

@@ -17,22 +17,28 @@
 #include "Domain/Amr/NewNeighborIds.hpp"
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/DirectionMap.hpp"
+#include "Domain/Structure/DirectionalIdMap.hpp"
 #include "Domain/Structure/Element.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Structure/Neighbors.hpp"
 #include "Domain/Structure/OrientationMap.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "Utilities/Algorithm.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 
 namespace amr {
 template <size_t VolumeDim>
-DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_parent(
+std::pair<DirectionMap<VolumeDim, Neighbors<VolumeDim>>,
+          DirectionalIdMap<VolumeDim, Mesh<VolumeDim>>>
+neighbors_of_parent(
     const ElementId<VolumeDim>& parent_id,
     const std::vector<std::tuple<
         const Element<VolumeDim>&,
         const std::unordered_map<ElementId<VolumeDim>, Info<VolumeDim>>&>>&
         children_elements_and_neighbor_info) {
-  DirectionMap<VolumeDim, Neighbors<VolumeDim>> result;
+  std::pair<DirectionMap<VolumeDim, Neighbors<VolumeDim>>,
+            DirectionalIdMap<VolumeDim, Mesh<VolumeDim>>>
+      result;
 
   std::vector<ElementId<VolumeDim>> children_ids;
   children_ids.reserve(children_elements_and_neighbor_info.size());
@@ -51,15 +57,19 @@ DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_parent(
           is_child(*(child_neighbors.ids().begin()))) {
         continue;  // neighbor in this direction was joined sibling
       }
-      if (0 == result.count(direction)) {
-        result.emplace(direction, Neighbors<VolumeDim>{
-                                      amr::new_neighbor_ids(
-                                          parent_id, direction, child_neighbors,
-                                          child_neighbor_info),
-                                      child_neighbors.orientation()});
+      const auto new_neighbor_ids_and_meshes = amr::new_neighbor_ids(
+          parent_id, direction, child_neighbors, child_neighbor_info);
+      std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids;
+      for (const auto& [id, mesh] : new_neighbor_ids_and_meshes) {
+        result.second.insert_or_assign({direction, id}, mesh);
+        new_neighbor_ids.insert(id);
+      }
+      if (0 == result.first.count(direction)) {
+        result.first.emplace(
+            direction, Neighbors<VolumeDim>{new_neighbor_ids,
+                                            child_neighbors.orientation()});
       } else {
-        result.at(direction).add_ids(amr::new_neighbor_ids(
-            parent_id, direction, child_neighbors, child_neighbor_info));
+        result.first.at(direction).add_ids(new_neighbor_ids);
       }
     }
   }
@@ -69,7 +79,9 @@ DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_parent(
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
 #define INSTANTIATE(_, data)                                                  \
-  template DirectionMap<DIM(data), Neighbors<DIM(data)>> neighbors_of_parent( \
+  template std::pair<DirectionMap<DIM(data), Neighbors<DIM(data)>>,           \
+                     DirectionalIdMap<DIM(data), Mesh<DIM(data)>>>            \
+  neighbors_of_parent(                                                        \
       const ElementId<DIM(data)>& parent_id,                                  \
       const std::vector<std::tuple<                                           \
           const Element<DIM(data)>&,                                          \

src/Domain/Amr/NeighborsOfParent.hpp

@@ -15,21 +15,27 @@ template <size_t VolumeDim>
 struct Info;
 }  // namespace amr
 template <size_t VolumeDim, typename T>
+class DirectionalIdMap;
+template <size_t VolumeDim, typename T>
 class DirectionMap;
 template <size_t VolumeDim>
 class Element;
 template <size_t VolumeDim>
 class ElementId;
 template <size_t VolumeDim>
+class Mesh;
+template <size_t VolumeDim>
 class Neighbors;
 /// \endcond
 
 namespace amr {
 /// \ingroup AmrGroup
-/// \brief returns the neighbors of the Element with ElementId `parent_id`,
-/// that is created from its `children_elements_and_neighbor_info`
+/// \brief returns the neighbors and their Mesh%es of the Element with ElementId
+/// `parent_id`, that is created from its `children_elements_and_neighbor_info`
 template <size_t VolumeDim>
-DirectionMap<VolumeDim, Neighbors<VolumeDim>> neighbors_of_parent(
+std::pair<DirectionMap<VolumeDim, Neighbors<VolumeDim>>,
+          DirectionalIdMap<VolumeDim, Mesh<VolumeDim>>>
+neighbors_of_parent(
     const ElementId<VolumeDim>& parent_id,
     const std::vector<std::tuple<
         const Element<VolumeDim>&,

src/Domain/Amr/NewNeighborIds.cpp

@@ -18,6 +18,7 @@
 #include "Domain/Structure/Neighbors.hpp"
 #include "Domain/Structure/OrientationMap.hpp"
 #include "Domain/Structure/SegmentId.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/StdHelpers.hpp"
@@ -39,12 +40,13 @@ bool overlapping_within_one_level(const SegmentId& segment1,
 
 namespace amr {
 template <size_t VolumeDim>
-std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids(
+std::unordered_map<ElementId<VolumeDim>, Mesh<VolumeDim>> new_neighbor_ids(
     const ElementId<VolumeDim>& my_id, const Direction<VolumeDim>& direction,
     const Neighbors<VolumeDim>& previous_neighbors_in_direction,
     const std::unordered_map<ElementId<VolumeDim>, Info<VolumeDim>>&
         previous_neighbors_amr_info) {
-  std::unordered_set<ElementId<VolumeDim>> new_neighbors_in_direction;
+  std::unordered_map<ElementId<VolumeDim>, Mesh<VolumeDim>>
+      new_neighbors_in_direction;
 
   const OrientationMap<VolumeDim>& orientation_map_from_me_to_neighbors =
       previous_neighbors_in_direction.orientation();
@@ -65,9 +67,10 @@ std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids(
                     ? previous_segment_id.id_of_child(
                           direction_to_me_in_neighbor_frame.side())
                     : previous_segment_id));
-    new_neighbors_in_direction.emplace(
-        previous_neighbor_id.block_id(),
-        std::array<SegmentId, 1>({{new_segment_id}}));
+    new_neighbors_in_direction.try_emplace(
+        {previous_neighbor_id.block_id(),
+         std::array<SegmentId, 1>({{new_segment_id}})},
+        previous_neighbors_amr_info.at(previous_neighbor_id).new_mesh);
     return new_neighbors_in_direction;
   }
 
@@ -87,6 +90,8 @@ std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids(
     // or children if I am the result of h-refinement)
     bool there_is_no_neighbor = false;
     const auto neighbor_segment_ids = previous_neighbor_id.segment_ids();
+    const auto& new_neighbor_mesh =
+        previous_neighbors_amr_info.at(previous_neighbor_id).new_mesh;
     for (size_t d = 0; d < VolumeDim; ++d) {
       const amr::Flag neighbor_flag =
           previous_neighbors_amr_info.at(previous_neighbor_id).flags.at(d);
@@ -155,16 +160,22 @@ std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids(
       if constexpr (VolumeDim > 1) {
         for (const auto segment_id_eta : valid_neighbor_segment_ids[1]) {
           if constexpr (VolumeDim == 2) {
-            new_neighbors_in_direction.emplace(
-                previous_neighbor_id.block_id(),
-                std::array{segment_id_xi, segment_id_eta},
-                previous_neighbor_id.grid_index());
+            // multiple previous neighbors may have joined to form a new
+            // neighbor; it is assumed they have the same new_mesh
+            new_neighbors_in_direction.try_emplace(
+                {previous_neighbor_id.block_id(),
+                 std::array{segment_id_xi, segment_id_eta},
+                 previous_neighbor_id.grid_index()},
+                new_neighbor_mesh);
           } else if constexpr (VolumeDim == 3) {
             for (const auto segment_id_zeta : valid_neighbor_segment_ids[2]) {
-              new_neighbors_in_direction.emplace(
-                  previous_neighbor_id.block_id(),
-                  std::array{segment_id_xi, segment_id_eta, segment_id_zeta},
-                  previous_neighbor_id.grid_index());
+              // multiple previous neighbors may have joined to form a new
+              // neighbor; it is assumed they have the same new_mesh
+              new_neighbors_in_direction.try_emplace(
+                  {previous_neighbor_id.block_id(),
+                   std::array{segment_id_xi, segment_id_eta, segment_id_zeta},
+                   previous_neighbor_id.grid_index()},
+                  new_neighbor_mesh);
             }
           }
         }
@@ -178,7 +189,8 @@ std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids(
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
 #define INSTANTIATE(_, data)                                           \
-  template std::unordered_set<ElementId<DIM(data)>> new_neighbor_ids(  \
+  template std::unordered_map<ElementId<DIM(data)>, Mesh<DIM(data)>>   \
+  new_neighbor_ids(                                                    \
       const ElementId<DIM(data)>& my_id,                               \
       const Direction<DIM(data)>& direction,                           \
       const Neighbors<DIM(data)>& previous_neighbors_in_direction,     \

src/Domain/Amr/NewNeighborIds.hpp

@@ -6,7 +6,6 @@
 #include <array>
 #include <cstddef>
 #include <unordered_map>
-#include <unordered_set>
 
 /// \cond
 namespace amr {
@@ -18,14 +17,16 @@ class Direction;
 template <size_t VolumeDim>
 class ElementId;
 template <size_t VolumeDim>
+class Mesh;
+template <size_t VolumeDim>
 class Neighbors;
 /// \endcond
 
 namespace amr {
 /// \ingroup AmrGroup
-/// \brief returns the ElementId%s of the neighbors in the given `direction` of
-/// the Element whose ElementId is `my_id` given the
-/// `previous_neighbors_in_direction` and their amr::Flag%s.
+/// \brief returns the ElementId and Mesh of the new neighbors in the given
+/// `direction` of the Element whose ElementId is `my_id` given the
+/// `previous_neighbors_in_direction` and their amr::Info\.
 ///
 /// \note `previous_neighbors_in_direction` should be from the parent (or a
 /// child) of the Element with `my_id` if `my_id` corresponds to a newly created
@@ -34,7 +35,7 @@ namespace amr {
 /// \note `previous_neighbors_amr_info` may contain info from neighbors in
 /// directions other than `direction`
 template <size_t VolumeDim>
-std::unordered_set<ElementId<VolumeDim>> new_neighbor_ids(
+std::unordered_map<ElementId<VolumeDim>, Mesh<VolumeDim>> new_neighbor_ids(
     const ElementId<VolumeDim>& my_id, const Direction<VolumeDim>& direction,
     const Neighbors<VolumeDim>& previous_neighbors_in_direction,
     const std::unordered_map<ElementId<VolumeDim>, Info<VolumeDim>>&

src/Domain/Tags/CMakeLists.txt

@@ -8,5 +8,6 @@ spectre_target_headers(
   ElementDistribution.hpp
   FaceNormal.hpp
   Faces.hpp
+  NeighborMesh.hpp
   SurfaceJacobian.hpp
   )

src/Evolution/DiscontinuousGalerkin/Tags/NeighborMesh.hpp → src/Domain/Tags/NeighborMesh.hpp

@@ -11,7 +11,7 @@
 #include "Domain/Structure/DirectionalIdMap.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 
-namespace evolution::dg::Tags {
+namespace domain::Tags {
 /*!
  * \brief Holds the mesh of each neighboring element.
  *
@@ -21,6 +21,6 @@ namespace evolution::dg::Tags {
  */
 template <size_t Dim>
 struct NeighborMesh : db::SimpleTag {
-  using type = DirectionalIdMap<Dim, Mesh<Dim>>;
+  using type = DirectionalIdMap<Dim, ::Mesh<Dim>>;
 };
-}  // namespace evolution::dg::Tags
+}  // namespace domain::Tags