From 183b788d08738d536d6ed784c4d1f61b83a9026a Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 29 Jun 2020 19:23:20 +0200
Subject: [PATCH] Add edge ref list transfer to diamond mesh

---
 src/mesh/DiamondDualMeshBuilder.cpp | 113 +++++++++++++++++++++++++---
 1 file changed, 103 insertions(+), 10 deletions(-)

diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 52857868c..d113bcc23 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -46,6 +46,17 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
     diamond_descriptor.cell_number_vector[i_primal_face] = primal_face_number[i_primal_face];
   }
 
+  if constexpr (Dimension == 3) {
+    const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
+    diamond_descriptor.edge_number_vector.resize(number_of_edges);
+    for (size_t i_edge = 0; i_edge < number_of_edges; ++i_edge) {
+      diamond_descriptor.edge_number_vector[i_edge] = i_edge;
+    }
+    if (parallel::size() > 1) {
+      throw NotImplementedError("parallel edge numbering is undefined");
+    }
+  }
+
   diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
 
   diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);
@@ -169,6 +180,45 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
   }
 
   {
+    const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] {
+      std::unordered_map<unsigned int, NodeId> node_to_id_map;
+      for (size_t i_node = 0; i_node < diamond_descriptor.node_number_vector.size(); ++i_node) {
+        node_to_id_map[diamond_descriptor.node_number_vector[i_node]] = i_node;
+      }
+      return node_to_id_map;
+    }();
+
+    for (size_t i_node_list = 0; i_node_list < primal_connectivity.template numberOfRefItemList<ItemType::node>();
+         ++i_node_list) {
+      const auto& primal_ref_node_list = primal_connectivity.template refItemList<ItemType::node>(i_node_list);
+      std::cout << "treating " << primal_ref_node_list.refId() << '\n';
+      const auto& primal_node_list = primal_ref_node_list.list();
+
+      const std::vector<NodeId> diamond_node_list = [&]() {
+        std::vector<NodeId> diamond_node_list;
+
+        for (size_t i_primal_node = 0; i_primal_node < primal_node_list.size(); ++i_primal_node) {
+          auto primal_node_id       = primal_node_list[i_primal_node];
+          const auto i_diamond_node = node_to_id_map.find(primal_connectivity.nodeNumber()[primal_node_id]);
+          if (i_diamond_node != node_to_id_map.end()) {
+            diamond_node_list.push_back(i_diamond_node->second);
+          }
+        }
+
+        return diamond_node_list;
+      }();
+
+      if (parallel::allReduceOr(diamond_node_list.size() > 0)) {
+        Array<NodeId> node_array(diamond_node_list.size());
+        for (size_t i = 0; i < diamond_node_list.size(); ++i) {
+          node_array[i] = diamond_node_list[i];
+        }
+        diamond_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), node_array});
+      }
+    }
+  }
+
+  if constexpr (Dimension > 1) {
     const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
 
     using Face = ConnectivityFace<Dimension>;
@@ -222,6 +272,58 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
     }
   }
 
+  if constexpr (Dimension > 2) {
+    const auto& primal_edge_to_node_matrix = primal_connectivity.edgeToNodeMatrix();
+    using Edge                             = ConnectivityFace<2>;
+
+    const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map = [&] {
+      std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
+      for (EdgeId l = 0; l < diamond_descriptor.edge_to_node_vector.size(); ++l) {
+        const auto& node_vector = diamond_descriptor.edge_to_node_vector[l];
+        edge_to_id_map[Edge(node_vector, diamond_descriptor.node_number_vector)] = l;
+      }
+      return edge_to_id_map;
+    }();
+
+    for (size_t i_edge_list = 0; i_edge_list < primal_connectivity.template numberOfRefItemList<ItemType::edge>();
+         ++i_edge_list) {
+      const auto& primal_ref_edge_list = primal_connectivity.template refItemList<ItemType::edge>(i_edge_list);
+      std::cout << "treating " << primal_ref_edge_list.refId() << '\n';
+      const auto& primal_edge_list = primal_ref_edge_list.list();
+
+      const std::vector<EdgeId> diamond_edge_list = [&]() {
+        std::vector<EdgeId> diamond_edge_list;
+        diamond_edge_list.reserve(primal_edge_list.size());
+        for (size_t i = 0; i < primal_edge_list.size(); ++i) {
+          EdgeId primal_edge_id = primal_edge_list[i];
+
+          const auto& primal_edge_node_list = primal_edge_to_node_matrix[primal_edge_id];
+
+          const auto i_diamond_edge = [&]() {
+            std::vector<unsigned int> node_list(primal_edge_node_list.size());
+            for (size_t i = 0; i < primal_edge_node_list.size(); ++i) {
+              node_list[i] = primal_edge_node_list[i];
+            }
+            return edge_to_id_map.find(Edge(node_list, diamond_descriptor.node_number_vector));
+          }();
+
+          if (i_diamond_edge != edge_to_id_map.end()) {
+            diamond_edge_list.push_back(i_diamond_edge->second);
+          }
+        }
+        return diamond_edge_list;
+      }();
+
+      if (parallel::allReduceOr(diamond_edge_list.size() > 0)) {
+        Array<EdgeId> edge_array(diamond_edge_list.size());
+        for (size_t i = 0; i < diamond_edge_list.size(); ++i) {
+          edge_array[i] = diamond_edge_list[i];
+        }
+        diamond_descriptor.addRefItemList(RefEdgeList{primal_ref_edge_list.refId(), edge_array});
+      }
+    }
+  }
+
   const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
   const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
 
@@ -261,16 +363,7 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
   }
 
   if constexpr (Dimension == 3) {
-    const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
-    diamond_descriptor.edge_number_vector.resize(number_of_edges);
-    for (size_t i_edge = 0; i_edge < number_of_edges; ++i_edge) {
-      diamond_descriptor.edge_number_vector[i_edge] = i_edge;
-    }
-    if (parallel::size() > 1) {
-      throw NotImplementedError("parallel edge numbering is undefined");
-    }
-
-    std::vector<int> edge_cell_owner(number_of_edges);
+    std::vector<int> edge_cell_owner(diamond_descriptor.edge_number_vector.size());
     std::fill(std::begin(edge_cell_owner), std::end(edge_cell_owner), parallel::size());
 
     for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
-- 
GitLab