Implement cache to speed up TriaAccessor::vertex_index()

diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h
index 9fbe2107da83d7ccd695ef8e3ce062dc4e9fccf6..deea012397158680928d81a7091c7bfb8f494526 100644 (file)
--- a/include/deal.II/grid/tria.h
+++ b/include/deal.II/grid/tria.h
@@ -3897,6 +3897,12 @@ private:
   void
   reset_global_cell_indices();
 
+  /**
+   * Reset cache for the cells' vertex indices.
+   */
+  void
+  reset_cell_vertex_indices_cache();
+
   /**
    * Refine all cells on all levels which were previously flagged for
    * refinement.
@@ -4255,6 +4261,7 @@ Triangulation<dim, spacedim>::load(Archive &ar, const unsigned int)
         level->global_active_cell_indices.resize(level->refine_flags.size());
         level->global_level_cell_indices.resize(level->refine_flags.size());
       }
+    reset_cell_vertex_indices_cache();
     reset_active_cell_indices();
     reset_global_cell_indices();
   }

diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h
index 1524e05e9f920e9faa1e8d0bcafefd3c9a8579d9..5e249501a28776724bdb73b95d3d1f60fa5bdc53 100644 (file)
--- a/include/deal.II/grid/tria_accessor.templates.h
+++ b/include/deal.II/grid/tria_accessor.templates.h
@@ -1112,8 +1112,23 @@ TriaAccessor<structdim, dim, spacedim>::vertex_index(
 {
   AssertIndexRange(corner, this->n_vertices());
 
-  return dealii::internal::TriaAccessorImplementation::Implementation::
-    vertex_index(*this, corner);
+  if (structdim == dim)
+    {
+      constexpr unsigned int max_vertices_per_cell = 1 << dim;
+      const std::size_t      my_index =
+        static_cast<std::size_t>(this->present_index) * max_vertices_per_cell;
+      AssertIndexRange(my_index + corner,
+                       this->tria->levels[this->present_level]
+                         ->cell_vertex_indices_cache.size());
+      const unsigned int vertex_index =
+        this->tria->levels[this->present_level]
+          ->cell_vertex_indices_cache[my_index + corner];
+      Assert(vertex_index != numbers::invalid_unsigned_int, ExcInternalError());
+      return vertex_index;
+    }
+  else
+    return dealii::internal::TriaAccessorImplementation::Implementation::
+      vertex_index(*this, corner);
 }
 
 
@@ -3206,7 +3221,8 @@ namespace internal
            ((i == 1) && cell.at_boundary(1) ?
               dealii::TriaAccessor<0, 1, spacedim>::right_vertex :
               dealii::TriaAccessor<0, 1, spacedim>::interior_vertex)),
-        cell.vertex_index(i));
+        dealii::internal::TriaAccessorImplementation::Implementation::
+          vertex_index(cell, i));
       return dealii::TriaIterator<dealii::TriaAccessor<0, 1, spacedim>>(a);
     }
 

diff --git a/include/deal.II/grid/tria_levels.h b/include/deal.II/grid/tria_levels.h
index 39969f954edcd0d75a4e077045a0eb7861de01d1..64f774477f19619f228a627bf0b2e1d05cf43a84 100644 (file)
--- a/include/deal.II/grid/tria_levels.h
+++ b/include/deal.II/grid/tria_levels.h
@@ -224,6 +224,14 @@ namespace internal
        */
       std::vector<dealii::ReferenceCell> reference_cell;
 
+      /**
+       * A cache for the vertex indices of the cells (`structdim == dim`), in
+       * order to more quickly retrieve these frequently accessed quantities.
+       * For simplified addressing, the information is indexed by the maximum
+       * number of vertices possible for a cell (quadrilateral/hexahedron).
+       */
+      std::vector<unsigned int> cell_vertex_indices_cache;
+
       /**
        * Determine an estimate for the memory consumption (in bytes) of this
        * object.
@@ -250,9 +258,9 @@ namespace internal
 
       ar &refine_flags &coarsen_flags;
 
-      // do not serialize 'active_cell_indices' here. instead of storing them
-      // to the stream and re-reading them again later, we just rebuild them
-      // in Triangulation::load()
+      // do not serialize `active_cell_indices` and `vertex_indices_cache`
+      // here. instead of storing them to the stream and re-reading them again
+      // later, we just rebuild them in Triangulation::load()
 
       ar &neighbors;
       ar &subdomain_ids;

diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index a5a42d27c104bee58260396f79b5ae0fdf2b4801..4f8c7d1b83ddf7fdc60f643cf4647f8705b52958 100644 (file)
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -10539,6 +10539,7 @@ Triangulation<dim, spacedim>::create_triangulation(
 
   // update our counts of the various elements of a triangulation, and set
   // active_cell_indices of all cells
+  reset_cell_vertex_indices_cache();
   internal::TriangulationImplementation::Implementation::compute_number_cache(
     *this, levels.size(), number_cache);
   reset_active_cell_indices();
@@ -13390,6 +13391,8 @@ Triangulation<dim, spacedim>::execute_coarsening_and_refinement()
 
   const DistortedCellList cells_with_distorted_children = execute_refinement();
 
+  reset_cell_vertex_indices_cache();
+
   // verify a case with which we have had
   // some difficulty in the past (see the
   // deal.II/coarsening_* tests)
@@ -13453,6 +13456,29 @@ Triangulation<dim, spacedim>::reset_global_cell_indices()
 
 
 
+template <int dim, int spacedim>
+void
+Triangulation<dim, spacedim>::reset_cell_vertex_indices_cache()
+{
+  for (unsigned int l = 0; l < levels.size(); ++l)
+    {
+      constexpr unsigned int     max_vertices_per_cell = 1 << dim;
+      std::vector<unsigned int> &cache = levels[l]->cell_vertex_indices_cache;
+      cache.clear();
+      cache.resize(levels[l]->refine_flags.size() * max_vertices_per_cell,
+                   numbers::invalid_unsigned_int);
+      for (const auto &cell : cell_iterators_on_level(l))
+        {
+          const unsigned int my_index = cell->index() * max_vertices_per_cell;
+          for (const unsigned int i : cell->vertex_indices())
+            cache[my_index + i] = internal::TriaAccessorImplementation::
+              Implementation::vertex_index(*cell, i);
+        }
+    }
+}
+
+
+
 template <int dim, int spacedim>
 void
 Triangulation<dim, spacedim>::update_periodic_face_map()