From d3aeb482fe34679cf448f3c3fb99e12e0f56d1ce Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Fri, 10 Jul 2020 08:35:35 +0200
Subject: [PATCH] Remove many GeometryInfos from grid/grid_tools

---
 source/grid/grid_tools.cc              | 169 +++++++++++--------------
 source/grid/grid_tools_dof_handlers.cc |  51 ++++----
 2 files changed, 101 insertions(+), 119 deletions(-)

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 766a65177e..ad2d211ef0 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -105,10 +105,9 @@ namespace GridTools
     const typename Triangulation<dim, spacedim>::active_cell_iterator endc =
       tria.end();
     for (; cell != endc; ++cell)
-      for (const unsigned int face : GeometryInfo<dim>::face_indices())
+      for (const unsigned int face : cell->face_indices())
         if (cell->face(face)->at_boundary())
-          for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
-               ++i)
+          for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
             boundary_vertices[cell->face(face)->vertex_index(i)] = true;
 
     // now traverse the list of boundary vertices and check distances.
@@ -361,8 +360,7 @@ namespace GridTools
       insert_face_data(const FaceIteratorType &face)
       {
         CellData<dim - 1> face_cell_data;
-        for (unsigned int vertex_n = 0;
-             vertex_n < GeometryInfo<dim>::vertices_per_face;
+        for (unsigned int vertex_n = 0; vertex_n < face->n_vertices();
              ++vertex_n)
           face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
         face_cell_data.boundary_id = face->boundary_id();
@@ -424,8 +422,7 @@ namespace GridTools
 
     unsigned int max_level_0_vertex_n = 0;
     for (const auto &cell : tria.cell_iterators_on_level(0))
-      for (const unsigned int cell_vertex_n :
-           GeometryInfo<dim>::vertex_indices())
+      for (const unsigned int cell_vertex_n : cell->vertex_indices())
         max_level_0_vertex_n =
           std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
     vertices.resize(max_level_0_vertex_n + 1);
@@ -438,8 +435,7 @@ namespace GridTools
       {
         // Save cell data
         CellData<dim> cell_data;
-        for (const unsigned int cell_vertex_n :
-             GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int cell_vertex_n : cell->vertex_indices())
           {
             Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
                    ExcInternalError());
@@ -455,20 +451,17 @@ namespace GridTools
         // Save face data
         if (dim > 1)
           {
-            for (const unsigned int face_n : GeometryInfo<dim>::face_indices())
+            for (const unsigned int face_n : cell->face_indices())
               face_data.insert_face_data(cell->face(face_n));
           }
         // Save line data
         if (dim == 3)
           {
-            for (unsigned int line_n = 0;
-                 line_n < GeometryInfo<dim>::lines_per_cell;
-                 ++line_n)
+            for (unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
               {
                 const auto  line = cell->line(line_n);
                 CellData<1> line_cell_data;
-                for (unsigned int vertex_n = 0;
-                     vertex_n < GeometryInfo<2>::vertices_per_face;
+                for (unsigned int vertex_n = 0; vertex_n < line->n_vertices();
                      ++vertex_n)
                   line_cell_data.vertices[vertex_n] =
                     line->vertex_index(vertex_n);
@@ -485,9 +478,8 @@ namespace GridTools
     {
       std::vector<bool> used_vertices(vertices.size());
       for (const CellData<dim> &cell_data : cells)
-        for (const unsigned int cell_vertex_n :
-             GeometryInfo<dim>::vertex_indices())
-          used_vertices[cell_data.vertices[cell_vertex_n]] = true;
+        for (const auto v : cell_data.vertices)
+          used_vertices[v] = true;
       Assert(std::find(used_vertices.begin(), used_vertices.end(), false) ==
                used_vertices.end(),
              ExcMessage("The level zero vertices should form a contiguous "
@@ -525,7 +517,7 @@ namespace GridTools
     // first check which vertices are actually used
     std::vector<bool> vertex_used(vertices.size(), false);
     for (unsigned int c = 0; c < cells.size(); ++c)
-      for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+      for (unsigned int v = 0; v < cells[c].vertices.size(); ++v)
         {
           Assert(cells[c].vertices[v] < vertices.size(),
                  ExcMessage("Invalid vertex index encountered! cells[" +
@@ -554,13 +546,15 @@ namespace GridTools
 
     // next replace old vertex numbers by the new ones
     for (unsigned int c = 0; c < cells.size(); ++c)
-      for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
-        cells[c].vertices[v] = new_vertex_numbers[cells[c].vertices[v]];
+      for (auto &v : cells[c].vertices)
+        v = new_vertex_numbers[v];
 
     // same for boundary data
     for (unsigned int c = 0; c < subcelldata.boundary_lines.size(); // NOLINT
          ++c)
-      for (const unsigned int v : GeometryInfo<1>::vertex_indices())
+      for (unsigned int v = 0;
+           v < subcelldata.boundary_lines[c].vertices.size();
+           ++v)
         {
           Assert(subcelldata.boundary_lines[c].vertices[v] <
                    new_vertex_numbers.size(),
@@ -580,7 +574,9 @@ namespace GridTools
 
     for (unsigned int c = 0; c < subcelldata.boundary_quads.size(); // NOLINT
          ++c)
-      for (const unsigned int v : GeometryInfo<2>::vertex_indices())
+      for (unsigned int v = 0;
+           v < subcelldata.boundary_quads[c].vertices.size();
+           ++v)
         {
           Assert(subcelldata.boundary_quads[c].vertices[v] <
                    new_vertex_numbers.size(),
@@ -920,7 +916,7 @@ namespace GridTools
       {
         // loop over all vertices of the cell and see if it is listed in the map
         // given as first argument of the function
-        for (const unsigned int vertex_no : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int vertex_no : cell->vertex_indices())
           {
             const unsigned int vertex_index = cell->vertex_index(vertex_no);
             const Point<dim> & vertex_point = cell->vertex(vertex_no);
@@ -960,7 +956,7 @@ namespace GridTools
     // according to the computed values
     std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
     for (const auto &cell : dof_handler.active_cell_iterators())
-      for (const unsigned int vertex_no : GeometryInfo<dim>::vertex_indices())
+      for (const unsigned int vertex_no : cell->vertex_indices())
         if (vertex_touched[cell->vertex_index(vertex_no)] == false)
           {
             Point<dim> &v = cell->vertex(vertex_no);
@@ -987,14 +983,13 @@ namespace GridTools
       endc = tria.end();
     for (; cell != endc; ++cell)
       {
-        for (unsigned int i : GeometryInfo<dim>::face_indices())
+        for (unsigned int i : cell->face_indices())
           {
             const typename Triangulation<dim, spacedim>::face_iterator &face =
               cell->face(i);
             if (face->at_boundary())
               {
-                for (unsigned j = 0; j < GeometryInfo<dim>::vertices_per_face;
-                     ++j)
+                for (unsigned j = 0; j < face->n_vertices(); ++j)
                   {
                     const Point<spacedim> &vertex       = face->vertex(j);
                     const unsigned int     vertex_index = face->vertex_index(j);
@@ -1055,8 +1050,7 @@ namespace GridTools
         {
           if (dim > 1)
             {
-              for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
-                   ++i)
+              for (unsigned int i = 0; i < cell->n_lines(); ++i)
                 {
                   const typename Triangulation<dim, spacedim>::line_iterator
                     line = cell->line(i);
@@ -1116,8 +1110,7 @@ namespace GridTools
         for (const auto &cell : triangulation.active_cell_iterators())
           if (cell->is_locally_owned())
             {
-              for (const unsigned int vertex_no :
-                   GeometryInfo<dim>::vertex_indices())
+              for (const unsigned int vertex_no : cell->vertex_indices())
                 {
                   const unsigned global_vertex_no =
                     cell->vertex_index(vertex_no);
@@ -1188,8 +1181,7 @@ namespace GridTools
 
         // now do the actual move of the vertices
         for (const auto &cell : triangulation.active_cell_iterators())
-          for (const unsigned int vertex_no :
-               GeometryInfo<dim>::vertex_indices())
+          for (const unsigned int vertex_no : cell->vertex_indices())
             cell->vertex(vertex_no) =
               new_vertex_locations[cell->vertex_index(vertex_no)];
       }
@@ -1209,7 +1201,7 @@ namespace GridTools
           endc = triangulation.end();
         for (; cell != endc; ++cell)
           if (!cell->is_artificial())
-            for (const unsigned int face : GeometryInfo<dim>::face_indices())
+            for (const unsigned int face : cell->face_indices())
               if (cell->face(face)->has_children() &&
                   !cell->face(face)->at_boundary())
                 {
@@ -1429,7 +1421,7 @@ namespace GridTools
     // list, but std::set throws out those cells already entered
     for (const auto &cell : mesh.active_cell_iterators())
       {
-        for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           if (cell->vertex_index(v) == vertex)
             {
               // OK, we found a cell that contains
@@ -1445,7 +1437,7 @@ namespace GridTools
                 for (unsigned int vface = 0; vface < dim; vface++)
                   {
                     const unsigned int face =
-                      GeometryInfo<dim>::vertex_to_face[v][vface];
+                      GeometryInfo<dim>::vertex_to_face[v][vface]; // TODO
 
                     if (!cell->at_boundary(face) &&
                         cell->neighbor(face)->is_active())
@@ -1471,7 +1463,7 @@ namespace GridTools
         // in 3d also loop over the edges
         if (dim >= 3)
           {
-            for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+            for (unsigned int e = 0; e < cell->n_lines(); ++e)
               if (cell->line(e)->has_children())
                 // the only place where this vertex could have been
                 // hiding is on the mid-edge point of the edge we
@@ -1756,7 +1748,7 @@ namespace GridTools
     double       minimum_distance = position.distance_square(vertices[0]);
     unsigned int closest_vertex   = 0;
 
-    for (unsigned int v = 1; v < GeometryInfo<dim>::vertices_per_cell; ++v)
+    for (unsigned int v = 1; v < cell->n_vertices(); ++v)
       {
         const double vertex_distance = position.distance_square(vertices[v]);
         if (vertex_distance < minimum_distance)
@@ -1811,8 +1803,7 @@ namespace GridTools
 
         for (; i < active_cells.size(); ++i)
           if (predicate(active_cells[i]))
-            for (const unsigned int v :
-                 GeometryInfo<MeshType::dimension>::vertex_indices())
+            for (const unsigned int v : active_cells[i]->vertex_indices())
               for (unsigned int d = 0; d < spacedim; ++d)
                 {
                   minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
@@ -2084,23 +2075,21 @@ namespace GridTools
       cell = triangulation.begin_active(),
       endc = triangulation.end();
     for (; cell != endc; ++cell)
-      for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+      for (const unsigned int i : cell->vertex_indices())
         vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
 
     // Take care of hanging nodes
     cell = triangulation.begin_active();
     for (; cell != endc; ++cell)
       {
-        for (unsigned int i : GeometryInfo<dim>::face_indices())
+        for (unsigned int i : cell->face_indices())
           {
             if ((cell->at_boundary(i) == false) &&
                 (cell->neighbor(i)->is_active()))
               {
                 typename Triangulation<dim, spacedim>::active_cell_iterator
                   adjacent_cell = cell->neighbor(i);
-                for (unsigned int j = 0;
-                     j < GeometryInfo<dim>::vertices_per_face;
-                     ++j)
+                for (unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
                   vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(
                     adjacent_cell);
               }
@@ -2109,7 +2098,7 @@ namespace GridTools
         // in 3d also loop over the edges
         if (dim == 3)
           {
-            for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell; ++i)
+            for (unsigned int i = 0; i < cell->n_lines(); ++i)
               if (cell->line(i)->has_children())
                 // the only place where this vertex could have been
                 // hiding is on the mid-edge point of the edge we
@@ -2169,7 +2158,7 @@ namespace GridTools
       {
         if (cell->is_locally_owned())
           {
-            for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+            for (const unsigned int i : cell->vertex_indices())
               {
                 types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
                 typename std::set<active_cell_iterator>::iterator
@@ -2237,7 +2226,7 @@ namespace GridTools
         // received.
         if (cell->is_ghost())
           {
-            for (unsigned int i : GeometryInfo<dim>::face_indices())
+            for (unsigned int i : cell->face_indices())
               {
                 if (cell->at_boundary(i) == false)
                   {
@@ -2252,7 +2241,7 @@ namespace GridTools
                               adjacent_cell->subdomain_id();
                             if (cell->subdomain_id() < adj_subdomain_id)
                               for (unsigned int j = 0;
-                                   j < GeometryInfo<dim>::vertices_per_face;
+                                   j < cell->face(i)->n_vertices();
                                    ++j)
                                 {
                                   vertices_to_recv[cell->subdomain_id()].insert(
@@ -2508,7 +2497,7 @@ namespace GridTools
       {
         const unsigned int index = cell->active_cell_index();
         cell_connectivity.add(index, index);
-        for (auto f : GeometryInfo<dim>::face_indices())
+        for (auto f : cell->face_indices())
           if ((cell->at_boundary(f) == false) &&
               (cell->neighbor(f)->has_children() == false))
             {
@@ -2535,7 +2524,7 @@ namespace GridTools
       triangulation.n_vertices());
     for (const auto &cell : triangulation.active_cell_iterators())
       {
-        for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           vertex_to_cell[cell->vertex_index(v)].push_back(
             cell->active_cell_index());
       }
@@ -2544,7 +2533,7 @@ namespace GridTools
                              triangulation.n_active_cells());
     for (const auto &cell : triangulation.active_cell_iterators())
       {
-        for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           for (unsigned int n = 0;
                n < vertex_to_cell[cell->vertex_index(v)].size();
                ++n)
@@ -2568,7 +2557,7 @@ namespace GridTools
          cell != triangulation.end(level);
          ++cell)
       {
-        for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
       }
 
@@ -2579,7 +2568,7 @@ namespace GridTools
          cell != triangulation.end(level);
          ++cell)
       {
-        for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           for (unsigned int n = 0;
                n < vertex_to_cell[cell->vertex_index(v)].size();
                ++n)
@@ -2998,7 +2987,7 @@ namespace GridTools
         if (cell->is_artificial() ||
             (cell->is_ghost() &&
              (cell->subdomain_id() < tr->locally_owned_subdomain())))
-          for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+          for (const unsigned int v : cell->vertex_indices())
             locally_owned_vertices[cell->vertex_index(v)] = false;
 
     return locally_owned_vertices;
@@ -3140,7 +3129,7 @@ namespace GridTools
         Tensor<spacedim - structdim, spacedim>
           parent_alternating_forms[GeometryInfo<structdim>::vertices_per_cell];
 
-        for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
+        for (const unsigned int i : object->vertex_indices())
           parent_vertices[i] = object->vertex(i);
 
         GeometryInfo<structdim>::alternating_form_at_vertices(
@@ -3168,7 +3157,7 @@ namespace GridTools
           [GeometryInfo<structdim>::vertices_per_cell];
 
         for (unsigned int c = 0; c < object->n_children(); ++c)
-          for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
+          for (const unsigned int i : object->child(c)->vertex_indices())
             child_vertices[c][i] = object->child(c)->vertex(i);
 
         // replace mid-object
@@ -3197,7 +3186,7 @@ namespace GridTools
         // objective function
         double objective = 0;
         for (unsigned int c = 0; c < object->n_children(); ++c)
-          for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
+          for (const unsigned int i : object->child(c)->vertex_indices())
             objective +=
               (child_alternating_forms[c][i] -
                average_parent_alternating_form / std::pow(2., 1. * structdim))
@@ -3285,10 +3274,8 @@ namespace GridTools
           Iterator::AccessorType::structure_dimension;
 
         double diameter = object->diameter();
-        for (const unsigned int f : GeometryInfo<structdim>::face_indices())
-          for (unsigned int e = f + 1;
-               e < GeometryInfo<structdim>::faces_per_cell;
-               ++e)
+        for (const unsigned int f : object->face_indices())
+          for (unsigned int e = f + 1; e < object->n_faces(); ++e)
             diameter = std::min(
               diameter,
               get_face_midpoint(object,
@@ -3425,7 +3412,7 @@ namespace GridTools
                         [GeometryInfo<structdim>::vertices_per_cell];
 
         for (unsigned int c = 0; c < object->n_children(); ++c)
-          for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
+          for (const unsigned int i : object->child(c)->vertex_indices())
             child_vertices[c][i] = object->child(c)->vertex(i);
 
         Tensor<spacedim - structdim, spacedim> child_alternating_forms
@@ -3439,9 +3426,8 @@ namespace GridTools
         old_min_product =
           child_alternating_forms[0][0] * parent_alternating_forms[0];
         for (unsigned int c = 0; c < object->n_children(); ++c)
-          for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
-            for (const unsigned int j :
-                 GeometryInfo<structdim>::vertex_indices())
+          for (const unsigned int i : object->child(c)->vertex_indices())
+            for (const unsigned int j : object->vertex_indices())
               old_min_product = std::min<double>(old_min_product,
                                                  child_alternating_forms[c][i] *
                                                    parent_alternating_forms[j]);
@@ -3463,9 +3449,8 @@ namespace GridTools
         new_min_product =
           child_alternating_forms[0][0] * parent_alternating_forms[0];
         for (unsigned int c = 0; c < object->n_children(); ++c)
-          for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
-            for (const unsigned int j :
-                 GeometryInfo<structdim>::vertex_indices())
+          for (const unsigned int i : object->child(c)->vertex_indices())
+            for (const unsigned int j : object->vertex_indices())
               new_min_product = std::min<double>(new_min_product,
                                                  child_alternating_forms[c][i] *
                                                    parent_alternating_forms[j]);
@@ -3506,7 +3491,7 @@ namespace GridTools
         // distorted but the neighbor is even more refined, then the face had
         // been deformed before already, and had been ignored at the time; we
         // should then also be able to ignore it this time as well
-        for (auto f : GeometryInfo<dim>::face_indices())
+        for (auto f : cell->face_indices())
           {
             Assert(cell->face(f)->has_children(), ExcInternalError());
             Assert(cell->face(f)->refinement_case() ==
@@ -3620,9 +3605,9 @@ namespace GridTools
     // and reset later
     if (dim >= 3)
       for (const auto &cell : tria.active_cell_iterators())
-        for (auto f : GeometryInfo<dim>::face_indices())
+        for (auto f : cell->face_indices())
           if (cell->face(f)->at_boundary())
-            for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_face; ++e)
+            for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
               {
                 const auto         bid = cell->face(f)->line(e)->boundary_id();
                 const unsigned int ind = std::find(src_boundary_ids.begin(),
@@ -3636,7 +3621,7 @@ namespace GridTools
 
     // now do cells
     for (const auto &cell : tria.active_cell_iterators())
-      for (auto f : GeometryInfo<dim>::face_indices())
+      for (auto f : cell->face_indices())
         if (cell->face(f)->at_boundary())
           {
             const auto         bid = cell->face(f)->boundary_id();
@@ -3653,8 +3638,7 @@ namespace GridTools
               }
 
             if (dim >= 3)
-              for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_face;
-                   ++e)
+              for (unsigned int e = 0; e < cell->face(f)->n_lines(); ++e)
                 {
                   const auto bid = cell->face(f)->line(e)->boundary_id();
                   const unsigned int ind = std::find(src_boundary_ids.begin(),
@@ -3683,7 +3667,7 @@ namespace GridTools
         cell->set_manifold_id(cell->material_id());
         if (compute_face_ids == true)
           {
-            for (auto f : GeometryInfo<dim>::face_indices())
+            for (auto f : cell->face_indices())
               {
                 if (cell->at_boundary(f) == false)
                   cell->face(f)->set_manifold_id(
@@ -3721,7 +3705,7 @@ namespace GridTools
     for (auto &cell : tria.active_cell_iterators())
       {
         if (dim > 1)
-          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          for (unsigned int l = 0; l < cell->n_lines(); ++l)
             {
               if (cell->line(l)->user_index() == 0)
                 {
@@ -3734,7 +3718,7 @@ namespace GridTools
                   cell->manifold_id());
             }
         if (dim > 2)
-          for (unsigned int l = 0; l < GeometryInfo<dim>::quads_per_cell; ++l)
+          for (unsigned int l = 0; l < cell->n_faces(); ++l)
             {
               if (cell->quad(l)->user_index() == 0)
                 {
@@ -3750,7 +3734,7 @@ namespace GridTools
     for (auto &cell : tria.active_cell_iterators())
       {
         if (dim > 1)
-          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          for (unsigned int l = 0; l < cell->n_lines(); ++l)
             {
               const auto id = cell->line(l)->user_index();
               // Make sure we change the manifold indicator only once
@@ -3765,7 +3749,7 @@ namespace GridTools
                 }
             }
         if (dim > 2)
-          for (unsigned int l = 0; l < GeometryInfo<dim>::quads_per_cell; ++l)
+          for (unsigned int l = 0; l < cell->n_faces(); ++l)
             {
               const auto id = cell->quad(l)->user_index();
               // Make sure we change the manifold indicator only once
@@ -3833,7 +3817,7 @@ namespace GridTools
         continue_refinement = false;
 
         for (const auto &cell : tria.active_cell_iterators())
-          for (const unsigned int j : GeometryInfo<dim>::face_indices())
+          for (const unsigned int j : cell->face_indices())
             if (cell->at_boundary(j) == false &&
                 cell->neighbor(j)->has_children())
               {
@@ -3905,7 +3889,7 @@ namespace GridTools
         for (const auto &cell : tria.active_cell_iterators())
           {
             unsigned int boundary_face_counter = 0;
-            for (auto f : GeometryInfo<dim>::face_indices())
+            for (auto f : cell->face_indices())
               if (cell->face(f)->at_boundary())
                 boundary_face_counter++;
             if (boundary_face_counter > dim)
@@ -4115,7 +4099,7 @@ namespace GridTools
         if (cells_to_remove[cell->active_cell_index()] == false)
           {
             CellData<dim> c;
-            for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+            for (const unsigned int v : cell->vertex_indices())
               c.vertices[v] = cell->vertex_index(v);
             c.manifold_id = cell->manifold_id();
             c.material_id = cell->material_id();
@@ -4132,19 +4116,19 @@ namespace GridTools
            face->manifold_id() != numbers::flat_manifold_id) &&
           faces_to_remove[face->index()] == false)
         {
-          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_face; ++l)
+          for (unsigned int l = 0; l < face->n_lines(); ++l)
             {
               CellData<1> line;
               if (dim == 2)
                 {
-                  for (const unsigned int v : GeometryInfo<1>::vertex_indices())
+                  for (const unsigned int v : face->vertex_indices())
                     line.vertices[v] = face->vertex_index(v);
                   line.boundary_id = face->boundary_id();
                   line.manifold_id = face->manifold_id();
                 }
               else
                 {
-                  for (const unsigned int v : GeometryInfo<1>::vertex_indices())
+                  for (const unsigned int v : face->line(l)->vertex_indices())
                     line.vertices[v] = face->line(l)->vertex_index(v);
                   line.boundary_id = face->line(l)->boundary_id();
                   line.manifold_id = face->line(l)->manifold_id();
@@ -4154,7 +4138,7 @@ namespace GridTools
           if (dim == 3)
             {
               CellData<2> quad;
-              for (const unsigned int v : GeometryInfo<2>::vertex_indices())
+              for (const unsigned int v : face->vertex_indices())
                 quad.vertices[v] = face->vertex_index(v);
               quad.boundary_id = face->boundary_id();
               quad.manifold_id = face->manifold_id();
@@ -5489,9 +5473,10 @@ namespace GridTools
             pair.second.first.first->face(pair.second.first.second);
           const auto mask = pair.second.second;
 
+          AssertDimension(face_a->n_vertices(), face_b->n_vertices());
+
           // loop over all vertices on face
-          for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
-               ++i)
+          for (unsigned int i = 0; i < face_a->n_vertices(); ++i)
             {
               const bool face_orientation = mask[0];
               const bool face_flip        = mask[1];
@@ -5586,7 +5571,7 @@ namespace GridTools
     std::vector<bool> vertex_of_own_cell(tria.n_vertices(), false);
     for (const auto &cell : tria.active_cell_iterators())
       if (cell->is_locally_owned())
-        for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           vertex_of_own_cell[cell->vertex_index(v)] = true;
 
     // 3) for each vertex belonging to a locally owned cell all ghost
@@ -5600,7 +5585,7 @@ namespace GridTools
           const types::subdomain_id owner = cell->subdomain_id();
 
           // loop over all its vertices
-          for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+          for (const unsigned int v : cell->vertex_indices())
             {
               // set owner if vertex belongs to a local cell
               if (vertex_of_own_cell[cell->vertex_index(v)])
diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc
index 8f782e61bf..84f8b1c9e0 100644
--- a/source/grid/grid_tools_dof_handlers.cc
+++ b/source/grid/grid_tools_dof_handlers.cc
@@ -253,7 +253,7 @@ namespace GridTools
     // list, but std::set throws out those cells already entered
     for (; cell != endc; ++cell)
       {
-        for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           if (cell->vertex_index(v) == vertex)
             {
               // OK, we found a cell that contains
@@ -295,7 +295,7 @@ namespace GridTools
         // in 3d also loop over the edges
         if (dim >= 3)
           {
-            for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
+            for (unsigned int e = 0; e < cell->n_lines(); ++e)
               if (cell->line(e)->has_children())
                 // the only place where this vertex could have been
                 // hiding is on the mid-edge point of the edge we
@@ -797,8 +797,7 @@ namespace GridTools
          cell != mesh.end(level);
          ++cell)
       if (predicate(cell)) // True predicate --> Part of subdomain
-        for (const unsigned int v :
-             GeometryInfo<MeshType::dimension>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
             true;
 
@@ -809,8 +808,7 @@ namespace GridTools
          cell != mesh.end(level);
          ++cell)
       if (!predicate(cell)) // False predicate --> Potential halo cell
-        for (const unsigned int v :
-             GeometryInfo<MeshType::dimension>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
                 v)] == true)
             {
@@ -909,8 +907,7 @@ namespace GridTools
     for (const auto &cell : mesh.active_cell_iterators())
       if (!predicate(cell)) // Negation of predicate --> Not Part of subdomain
         {
-          for (const unsigned int v :
-               GeometryInfo<MeshType::dimension>::vertex_indices())
+          for (const unsigned int v : cell->vertex_indices())
             vertices_outside_subdomain[cell->vertex_index(v)] = true;
           n_non_predicate_cells++;
         }
@@ -928,8 +925,7 @@ namespace GridTools
     for (const auto &cell : mesh.active_cell_iterators())
       if (predicate(cell)) // True predicate --> Potential boundary cell of the
                            // subdomain
-        for (const unsigned int v :
-             GeometryInfo<MeshType::dimension>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           if (vertices_outside_subdomain[cell->vertex_index(v)] == true)
             {
               subdomain_boundary_cells.push_back(cell);
@@ -1094,8 +1090,7 @@ namespace GridTools
     // if it belongs to the predicate domain, extend the bounding box.
     for (const auto &cell : mesh.active_cell_iterators())
       if (predicate(cell)) // True predicate --> Part of subdomain
-        for (const unsigned int v :
-             GeometryInfo<MeshType::dimension>::vertex_indices())
+        for (const unsigned int v : cell->vertex_indices())
           if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
               false)
             {
@@ -1238,9 +1233,13 @@ namespace GridTools
                                                            mesh_2.begin(0),
                                                          endc = mesh_1.end(0);
     for (; cell_1 != endc; ++cell_1, ++cell_2)
-      for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
-        if (cell_1->vertex(v) != cell_2->vertex(v))
+      {
+        if (cell_1->n_vertices() != cell_2->n_vertices())
           return false;
+        for (const unsigned int v : cell_1->vertex_indices())
+          if (cell_1->vertex(v) != cell_2->vertex(v))
+            return false;
+      }
 
     // if we've gotten through all
     // this, then the meshes really
@@ -1407,8 +1406,7 @@ namespace GridTools
 
     std::vector<typename MeshType::active_cell_iterator> patch;
     patch.push_back(cell);
-    for (const unsigned int face_number :
-         GeometryInfo<MeshType::dimension>::face_indices())
+    for (const unsigned int face_number : cell->face_indices())
       if (cell->face(face_number)->at_boundary() == false)
         {
           if (cell->neighbor(face_number)->has_children() == false)
@@ -1517,8 +1515,7 @@ namespace GridTools
          uniform_cell != uniform_cells.end();
          ++uniform_cell)
       {
-        for (const unsigned int v :
-             GeometryInfo<Container::dimension>::vertex_indices())
+        for (const unsigned int v : (*uniform_cell)->vertex_indices())
           {
             Point<Container::space_dimension> position =
               (*uniform_cell)->vertex(v);
@@ -1594,7 +1591,7 @@ namespace GridTools
                       // adjust the cell vertices of the local_triangulation to
                       // match cell vertices of the global triangulation
                       for (const unsigned int v :
-                           GeometryInfo<Container::dimension>::vertex_indices())
+                           active_tria_cell->vertex_indices())
                         active_tria_cell->vertex(v) = patch[i]->vertex(v);
 
                       Assert(active_tria_cell->center().distance(
@@ -1760,8 +1757,7 @@ namespace GridTools
             // few and we don't have to use them.
             if (cell->is_artificial() == false)
               {
-                for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
-                     ++l)
+                for (unsigned int l = 0; l < cell->n_lines(); ++l)
                   if (cell->line(l)->has_children())
                     for (unsigned int c = 0; c < cell->line(l)->n_children();
                          ++c)
@@ -1810,7 +1806,7 @@ namespace GridTools
             // face (or line).
 
             // Take care of dofs on neighbor faces
-            for (const unsigned int f : GeometryInfo<dim>::face_indices())
+            for (const unsigned int f : cell->face_indices())
               {
                 if (cell->face(f)->has_children())
                   {
@@ -1917,8 +1913,7 @@ namespace GridTools
             // and dofs on all children of parent line.
             if (dim == 3)
               {
-                for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
-                     ++l)
+                for (unsigned int l = 0; l < cell->n_lines(); ++l)
                   {
                     if (cell->line(l)->has_children())
                       {
@@ -2227,7 +2222,7 @@ namespace GridTools
          cell != mesh.end(0);
          ++cell)
       {
-        for (unsigned int i : GeometryInfo<dim>::face_indices())
+        for (unsigned int i : cell->face_indices())
           {
             const typename MeshType::face_iterator face = cell->face(i);
             if (face->at_boundary() && face->boundary_id() == b_id1)
@@ -2442,11 +2437,13 @@ namespace GridTools
 
     std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
 
+    AssertDimension(face1->n_vertices(), face2->n_vertices());
+
     std::set<unsigned int> face2_vertices;
-    for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
+    for (unsigned int i = 0; i < face1->n_vertices(); ++i)
       face2_vertices.insert(i);
 
-    for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
+    for (unsigned int i = 0; i < face1->n_vertices(); ++i)
       {
         for (std::set<unsigned int>::iterator it = face2_vertices.begin();
              it != face2_vertices.end();
-- 
2.39.5