Use Cache::get_vertex_to_cell_map() in GT::find_all_locally_owned_act…

include/deal.II/grid/grid_tools.h

@@ -1810,6 +1810,11 @@ namespace GridTools
    * corresponding neighboring cells with points in unit coordinates are also
    * identified.
    *
+   * The parameter @p vertex_to_cells allows to accelerate the process of
+   * identifying the neighbors of a cell, by first precomputing a map from the
+   * vertex indices to the cells. Such data structure is, e.g., provided by
+   * GridTools::Cache::get_vertex_to_cell_map().
+   *
    * This function is useful e.g. for discontinuous function spaces where, for
    * the case the given point `p` lies on a vertex, edge or face, several
    * cells might hold independent values of the solution that get combined in
@@ -1843,7 +1848,10 @@ namespace GridTools
       const Point<spacedim> &        p,
       const double                   tolerance,
       const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
-                      Point<dim>> &  first_cell);
+                      Point<dim>> &  first_cell,
+      const std::vector<
+        std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
+        *vertex_to_cells = nullptr);
 
   /**
    * A variant of the previous function that internally calls one of the

source/grid/grid_tools.cc

@@ -6073,7 +6073,8 @@ namespace GridTools
               cache.get_triangulation(),
               point,
               tolerance,
-              first_cell);
+              first_cell,
+              &cache.get_vertex_to_cell_map());
 
           if (enforce_unique_mapping)
             {

source/grid/grid_tools_dof_handlers.cc

@@ -615,7 +615,10 @@ namespace GridTools
       const Point<spacedim> &        p,
       const double                   tolerance,
       const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
-                      Point<dim>> &  first_cell)
+                      Point<dim>> &  first_cell,
+      const std::vector<
+        std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
+        *vertex_to_cells)
   {
     std::vector<
       std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
@@ -671,14 +674,19 @@ namespace GridTools
         unsigned int local_vertex_index = 0;
         for (unsigned int d = 0; d < dim; ++d)
           local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
-        std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
-          cells = find_cells_adjacent_to_vertex(
-            mesh, my_cell->vertex_index(local_vertex_index));
-        for (const auto &cell : cells)
-          {
+
+        const auto fu = [&](const auto &tentative_cells) {
+          for (const auto &cell : tentative_cells)
             if (cell != my_cell)
               cells_to_add.push_back(cell);
-          }
+        };
+
+        const auto vertex_index = my_cell->vertex_index(local_vertex_index);
+
+        if (vertex_to_cells != nullptr)
+          fu((*vertex_to_cells)[vertex_index]);
+        else
+          fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
       }
     // point on line in 3d: We cannot simply take the intersection between
     // the two vertices of cells because of hanging nodes. So instead we
@@ -712,21 +720,28 @@ namespace GridTools
           (vertex_indices[1].second << vertex_indices[1].first);
         for (unsigned int d = 0; d < 2; ++d)
           {
-            auto tentative_cells = find_cells_adjacent_to_vertex(
-              mesh,
-              my_cell->vertex_index(first_vertex + (d << free_direction)));
-            for (const auto &cell : tentative_cells)
-              {
-                bool cell_not_yet_present = true;
-                for (const auto &other_cell : cells_to_add)
-                  if (cell == other_cell)
-                    {
-                      cell_not_yet_present = false;
-                      break;
-                    }
-                if (cell_not_yet_present)
-                  cells_to_add.push_back(cell);
-              }
+            const auto fu = [&](const auto &tentative_cells) {
+              for (const auto &cell : tentative_cells)
+                {
+                  bool cell_not_yet_present = true;
+                  for (const auto &other_cell : cells_to_add)
+                    if (cell == other_cell)
+                      {
+                        cell_not_yet_present = false;
+                        break;
+                      }
+                  if (cell_not_yet_present)
+                    cells_to_add.push_back(cell);
+                }
+            };
+
+            const auto vertex_index =
+              my_cell->vertex_index(first_vertex + (d << free_direction));
+
+            if (vertex_to_cells != nullptr)
+              fu((*vertex_to_cells)[vertex_index]);
+            else
+              fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
           }
       }
 

source/grid/grid_tools_dof_handlers.inst.in

@@ -69,7 +69,8 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS;
         const Point<deal_II_space_dimension> &,
         const double,
         const std::pair<typename X::active_cell_iterator,
-                        Point<deal_II_dimension>> &);
+                        Point<deal_II_dimension>> &,
+        const std::vector<std::set<typename X::active_cell_iterator>> *);
 
       template std::vector<
         std::pair<dealii::internal::ActiveCellIterator<deal_II_dimension,