From 3c9bae0954d8559a21002382a7e56b38a1ad2648 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@rub.de>
Date: Fri, 4 Oct 2024 12:55:42 +0200
Subject: [PATCH] Increase performance of
 GridTools::get_vertex_connectivity_of_cells

---
 source/grid/grid_tools_topology.cc | 46 ++++++++++++++++++------------
 1 file changed, 28 insertions(+), 18 deletions(-)

diff --git a/source/grid/grid_tools_topology.cc b/source/grid/grid_tools_topology.cc
index a98e3f9f86..97b88c422e 100644
--- a/source/grid/grid_tools_topology.cc
+++ b/source/grid/grid_tools_topology.cc
@@ -1855,8 +1855,9 @@ namespace GridTools
     const Triangulation<dim, spacedim> &triangulation,
     DynamicSparsityPattern             &cell_connectivity)
   {
-    std::vector<std::vector<unsigned int>> vertex_to_cell(
-      triangulation.n_vertices());
+    // Assume at most 16 neighbors per cell (a structured hex mesh has 8)
+    std::vector<boost::container::small_vector<unsigned int, 16>>
+      vertex_to_cell(triangulation.n_vertices());
     for (const auto &cell : triangulation.active_cell_iterators())
       {
         for (const unsigned int v : cell->vertex_indices())
@@ -1866,14 +1867,20 @@ namespace GridTools
 
     cell_connectivity.reinit(triangulation.n_active_cells(),
                              triangulation.n_active_cells());
+    std::vector<types::global_dof_index> neighbors;
     for (const auto &cell : triangulation.active_cell_iterators())
       {
+        neighbors.clear();
         for (const unsigned int v : cell->vertex_indices())
-          for (unsigned int n = 0;
-               n < vertex_to_cell[cell->vertex_index(v)].size();
-               ++n)
-            cell_connectivity.add(cell->active_cell_index(),
-                                  vertex_to_cell[cell->vertex_index(v)][n]);
+          neighbors.insert(neighbors.end(),
+                           vertex_to_cell[cell->vertex_index(v)].begin(),
+                           vertex_to_cell[cell->vertex_index(v)].end());
+        std::sort(neighbors.begin(), neighbors.end());
+        cell_connectivity.add_entries(cell->active_cell_index(),
+                                      neighbors.begin(),
+                                      std::unique(neighbors.begin(),
+                                                  neighbors.end()),
+                                      true);
       }
   }
 
@@ -1885,8 +1892,8 @@ namespace GridTools
     const unsigned int                  level,
     DynamicSparsityPattern             &cell_connectivity)
   {
-    std::vector<std::vector<unsigned int>> vertex_to_cell(
-      triangulation.n_vertices());
+    std::vector<boost::container::small_vector<unsigned int, 16>>
+      vertex_to_cell(triangulation.n_vertices());
     for (typename Triangulation<dim, spacedim>::cell_iterator cell =
            triangulation.begin(level);
          cell != triangulation.end(level);
@@ -1898,17 +1905,20 @@ namespace GridTools
 
     cell_connectivity.reinit(triangulation.n_cells(level),
                              triangulation.n_cells(level));
-    for (typename Triangulation<dim, spacedim>::cell_iterator cell =
-           triangulation.begin(level);
-         cell != triangulation.end(level);
-         ++cell)
+    std::vector<types::global_dof_index> neighbors;
+    for (const auto &cell : triangulation.cell_iterators_on_level(level))
       {
+        neighbors.clear();
         for (const unsigned int v : cell->vertex_indices())
-          for (unsigned int n = 0;
-               n < vertex_to_cell[cell->vertex_index(v)].size();
-               ++n)
-            cell_connectivity.add(cell->index(),
-                                  vertex_to_cell[cell->vertex_index(v)][n]);
+          neighbors.insert(neighbors.end(),
+                           vertex_to_cell[cell->vertex_index(v)].begin(),
+                           vertex_to_cell[cell->vertex_index(v)].end());
+        std::sort(neighbors.begin(), neighbors.end());
+        cell_connectivity.add_entries(cell->index(),
+                                      neighbors.begin(),
+                                      std::unique(neighbors.begin(),
+                                                  neighbors.end()),
+                                      true);
       }
   }
 } /* namespace GridTools */
-- 
2.39.5