From 89952e7de968b0b2bbf6f9d125c5ff308e2cd6b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 3 Feb 2025 19:09:01 +0100
Subject: [PATCH] Remove m_cell_global_index from Connectivity

- now this information is computed when needed (load balance)
- closes #11
---
 src/mesh/Connectivity.cpp           | 86 ++++++++++++++++++++++++++---
 src/mesh/Connectivity.hpp           | 51 +----------------
 src/mesh/ConnectivityDispatcher.cpp |  2 +-
 3 files changed, 80 insertions(+), 59 deletions(-)

diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 65a880495..44c07b628 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -46,14 +46,6 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
   m_cell_number = WeakCellValue<const int>(*this, descriptor.cellNumberVector());
   m_node_number = WeakNodeValue<const int>(*this, descriptor.nodeNumberVector());
 
-  {
-    WeakCellValue<int> cell_global_index(*this);
-    int first_index = 0;
-    parallel_for(
-      this->numberOfCells(), PUGS_LAMBDA(CellId j) { cell_global_index[j] = first_index + j; });
-    m_cell_global_index = cell_global_index;
-  }
-
   m_cell_owner = WeakCellValue<const int>(*this, descriptor.cellOwnerVector());
 
   {
@@ -310,6 +302,80 @@ Connectivity<Dimension>::_write(std::ostream& os) const
   return os;
 }
 
+template <size_t Dimension>
+CRSGraph
+Connectivity<Dimension>::ownCellToCellGraph() const
+{
+  std::vector<std::vector<int>> cell_cells(this->numberOfCells());
+  const auto& face_to_cell_matrix = this->faceToCellMatrix();
+
+  for (FaceId l = 0; l < this->numberOfFaces(); ++l) {
+    const auto& face_to_cell = face_to_cell_matrix[l];
+    if (face_to_cell.size() > 1) {
+      const CellId cell_0 = face_to_cell[0];
+      const CellId cell_1 = face_to_cell[1];
+
+      if (m_cell_is_owned[cell_0]) {
+        cell_cells[cell_0].push_back(cell_1);
+      }
+      if (m_cell_is_owned[cell_1]) {
+        cell_cells[cell_1].push_back(cell_0);
+      }
+    }
+  }
+
+  std::vector<std::vector<int>> owned_cell_cells;
+  owned_cell_cells.reserve(this->numberOfCells());
+
+  size_t number_of_owned_cells = 0;
+  for (CellId cell_id = 0; cell_id < this->numberOfCells(); ++cell_id) {
+    if (m_cell_is_owned[cell_id]) {
+      ++number_of_owned_cells;
+    }
+  }
+
+  for (size_t i_cell = 0; i_cell < this->numberOfCells(); ++i_cell) {
+    if (cell_cells[i_cell].size() > 0) {
+      owned_cell_cells.emplace_back(std::move(cell_cells[i_cell]));
+    }
+  }
+
+  Array number_of_owned_cells_by_rank = parallel::allGather(number_of_owned_cells);
+
+  // index is defined locally. First we compute the first index on the
+  // processor
+  size_t global_index = 0;
+  for (size_t i = 0; i < parallel::rank(); ++i) {
+    global_index += number_of_owned_cells_by_rank[i];
+  }
+
+  CellValue<int> global_cell_index{*this};
+  for (CellId cell_id = 0; cell_id < this->numberOfCells(); ++cell_id) {
+    if (m_cell_is_owned[cell_id]) {
+      global_cell_index[cell_id] = global_index++;
+    }
+  }
+  synchronize(global_cell_index);
+
+  Array<int> entries(owned_cell_cells.size() + 1);
+  entries[0] = 0;
+  for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) {
+    entries[i_cell + 1] = entries[i_cell] + owned_cell_cells[i_cell].size();
+  }
+  Array<int> neighbors(entries[owned_cell_cells.size()]);
+  {
+    size_t i = 0;
+    for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) {
+      for (CellId cell_id : owned_cell_cells[i_cell]) {
+        neighbors[i] = global_cell_index[cell_id];
+        ++i;
+      }
+    }
+  }
+
+  return CRSGraph(entries, neighbors);
+}
+
 template std::ostream& Connectivity<1>::_write(std::ostream&) const;
 template std::ostream& Connectivity<2>::_write(std::ostream&) const;
 template std::ostream& Connectivity<3>::_write(std::ostream&) const;
@@ -333,3 +399,7 @@ template Connectivity<3>::Connectivity();
 template void Connectivity<1>::_buildFrom(const ConnectivityDescriptor&);
 template void Connectivity<2>::_buildFrom(const ConnectivityDescriptor&);
 template void Connectivity<3>::_buildFrom(const ConnectivityDescriptor&);
+
+template CRSGraph Connectivity<1>::ownCellToCellGraph() const;
+template CRSGraph Connectivity<2>::ownCellToCellGraph() const;
+template CRSGraph Connectivity<3>::ownCellToCellGraph() const;
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index faee4c28c..6d89a820f 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -59,8 +59,6 @@ class Connectivity final : public IConnectivity
   ConnectivityMatrix m_item_to_item_matrix[Dimension + 1][Dimension + 1];
   WeakCellValue<const CellType> m_cell_type;
 
-  WeakCellValue<const int> m_cell_global_index;
-
   WeakCellValue<const int> m_cell_number;
   WeakFaceValue<const int> m_face_number;
   WeakEdgeValue<const int> m_edge_number;
@@ -655,54 +653,7 @@ class Connectivity final : public IConnectivity
     }
   }
 
-  PUGS_INLINE
-  CRSGraph
-  cellToCellGraph() const
-  {
-    std::vector<std::vector<int>> cell_cells(this->numberOfCells());
-    const auto& face_to_cell_matrix = this->faceToCellMatrix();
-
-    for (FaceId l = 0; l < this->numberOfFaces(); ++l) {
-      const auto& face_to_cell = face_to_cell_matrix[l];
-      if (face_to_cell.size() > 1) {
-        const CellId cell_0 = face_to_cell[0];
-        const CellId cell_1 = face_to_cell[1];
-
-        if (m_cell_is_owned[cell_0]) {
-          cell_cells[cell_0].push_back(cell_1);
-        }
-        if (m_cell_is_owned[cell_1]) {
-          cell_cells[cell_1].push_back(cell_0);
-        }
-      }
-    }
-
-    std::vector<std::vector<int>> owned_cell_cells;
-    owned_cell_cells.reserve(this->numberOfCells());
-
-    for (size_t i_cell = 0; i_cell < this->numberOfCells(); ++i_cell) {
-      if (cell_cells[i_cell].size() > 0) {
-        owned_cell_cells.emplace_back(std::move(cell_cells[i_cell]));
-      }
-    }
-
-    Array<int> entries(owned_cell_cells.size() + 1);
-    entries[0] = 0;
-    for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) {
-      entries[i_cell + 1] = entries[i_cell] + owned_cell_cells[i_cell].size();
-    }
-    Array<int> neighbors(entries[owned_cell_cells.size()]);
-    {
-      size_t i = 0;
-      for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) {
-        for (CellId cell_id : owned_cell_cells[i_cell]) {
-          neighbors[i] = m_cell_global_index[cell_id];
-          ++i;
-        }
-      }
-    }
-    return CRSGraph(entries, neighbors);
-  }
+  CRSGraph ownCellToCellGraph() const;
 
   PUGS_INLINE
   size_t
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index 65133fda2..8c205a871 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -13,7 +13,7 @@ void
 ConnectivityDispatcher<Dimension>::_buildNewOwner()
 {
   if constexpr (item_type == ItemType::cell) {
-    CRSGraph connectivity_graph = m_connectivity.cellToCellGraph();
+    CRSGraph connectivity_graph = m_connectivity.ownCellToCellGraph();
     Partitioner P;
     Array new_owner_array   = P.partition(connectivity_graph);
     CellValue cell_is_owned = m_connectivity.cellIsOwned();
-- 
GitLab