From 78d7c6af4d2fd704c59fa0c21c376cbad140853f Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 4 Feb 2020 14:02:41 -0700
Subject: [PATCH] Convert more places to use GeometryInfo::face_indices().

---
 include/deal.II/base/geometry_info.h          |  4 +--
 include/deal.II/grid/grid_tools.h             |  6 ++--
 .../deal.II/grid/tria_accessor.templates.h    |  3 +-
 .../matrix_free/cuda_hanging_nodes_internal.h |  3 +-
 .../deal.II/matrix_free/face_setup_internal.h |  6 ++--
 .../matrix_free/mapping_info.templates.h      |  3 +-
 .../matrix_free/matrix_free.templates.h       |  3 +-
 .../numerics/error_estimator.templates.h      |  3 +-
 .../numerics/matrix_creator.templates.h       | 12 +++-----
 source/dofs/dof_tools_constraints.cc          |  9 ++----
 source/dofs/dof_tools_sparsity.cc             |  2 +-
 source/fe/fe_enriched.cc                      |  3 +-
 source/fe/fe_nedelec.cc                       |  3 +-
 source/fe/mapping_fe_field.cc                 |  6 ++--
 source/grid/grid_generator.cc                 |  3 +-
 source/grid/grid_out.cc                       |  6 ++--
 source/grid/grid_tools.cc                     |  5 ++--
 source/grid/tria.cc                           | 30 +++++++------------
 source/grid/tria_description.cc               |  3 +-
 source/multigrid/mg_tools.cc                  | 18 ++++-------
 20 files changed, 45 insertions(+), 86 deletions(-)

diff --git a/include/deal.II/base/geometry_info.h b/include/deal.II/base/geometry_info.h
index cb51154011..34766592b4 100644
--- a/include/deal.II/base/geometry_info.h
+++ b/include/deal.II/base/geometry_info.h
@@ -1211,7 +1211,7 @@ struct GeometryInfo<0>
    * using range-based for loops of the following kind:
    * @code
    *   for (auto &cell : triangulation.active_cell_iterators())
-   *     for (auto face_index : GeometryInfo<dim>::face_indices)
+   *     for (auto face_index : GeometryInfo<dim>::face_indices())
    *       if (cell->face(face_index)->at_boundary())
    *         ... do something ...
    * @endcode
@@ -1863,7 +1863,7 @@ struct GeometryInfo
    * using range-based for loops of the following kind:
    * @code
    *   for (auto &cell : triangulation.active_cell_iterators())
-   *     for (auto face_index : GeometryInfo<dim>::face_indices)
+   *     for (auto face_index : GeometryInfo<dim>::face_indices())
    *       if (cell->face(face_index)->at_boundary())
    *         ... do something ...
    * @endcode
diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 5820941335..a6634986ee 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -3162,8 +3162,7 @@ namespace GridTools
           cell = triangulation.begin_active(),
           endc = triangulation.end();
         for (; cell != endc; ++cell)
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (const unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->has_children() &&
                 !cell->face(face)->at_boundary())
               {
@@ -3179,8 +3178,7 @@ namespace GridTools
           cell = triangulation.begin_active(),
           endc = triangulation.end();
         for (; cell != endc; ++cell)
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (const unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->has_children() &&
                 !cell->face(face)->at_boundary())
               {
diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h
index aabb4fc29d..572cfde36a 100644
--- a/include/deal.II/grid/tria_accessor.templates.h
+++ b/include/deal.II/grid/tria_accessor.templates.h
@@ -3211,8 +3211,7 @@ inline unsigned int
 CellAccessor<dim, spacedim>::face_iterator_to_index(
   const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> &face) const
 {
-  for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell;
-       ++face_n)
+  for (const unsigned int face_n : GeometryInfo<dim>::face_indices())
     if (this->face(face_n) == face)
       return face_n;
 
diff --git a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h
index 7cc5d2c31c..a4c87d51d1 100644
--- a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h
+++ b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h
@@ -285,8 +285,7 @@ namespace CUDAWrappers
         FETools::lexicographic_to_hierarchic_numbering<dim - 1>(
           FE_Q<dim - 1>(fe_degree));
 
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const unsigned int face : GeometryInfo<dim>::face_indices())
         {
           if ((!cell->at_boundary(face)) &&
               (cell->neighbor(face)->has_children() == false))
diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index 2cfd03a943..2ad4ab8437 100644
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -210,8 +210,7 @@ namespace internal
                 continue;
               typename dealii::Triangulation<dim>::cell_iterator dcell(
                 &triangulation, cell_levels[i].first, cell_levels[i].second);
-              for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
-                   ++f)
+              for (const unsigned int f : GeometryInfo<dim>::face_indices())
                 {
                   if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
                     continue;
@@ -794,8 +793,7 @@ namespace internal
                   &triangulation,
                   cell_levels[cell].first,
                   cell_levels[cell].second);
-                for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
-                     ++f)
+                for (const unsigned int f : GeometryInfo<dim>::face_indices())
                   {
                     // boundary face
                     if (face_is_owned[dcell->face(f)->index()] ==
diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index b31ed2dcb0..f61e09d0ab 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -2014,8 +2014,7 @@ namespace internal
         fe_face_values_neigh[i].resize(face_data_by_cells[i].descriptor.size());
       for (unsigned int cell = 0; cell < cell_type.size(); ++cell)
         for (unsigned int my_q = 0; my_q < face_data_by_cells.size(); ++my_q)
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (const unsigned int face : GeometryInfo<dim>::face_indices())
             {
               if (fe_face_values[my_q][fe_index].get() == nullptr)
                 fe_face_values[my_q][fe_index].reset(
diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 674a27a9bc..d4de7d1ba1 100644
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -1622,8 +1622,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
                                          const bool,
                                          const bool)> &fu) {
               for (unsigned int c = 0; c < n_cell_batches(); ++c)
-                for (unsigned int d = 0; d < GeometryInfo<dim>::faces_per_cell;
-                     ++d)
+                for (const unsigned int d : GeometryInfo<dim>::face_indices())
                   for (unsigned int v = 0;
                        v < VectorizedArrayType::n_array_elements;
                        ++v)
diff --git a/include/deal.II/numerics/error_estimator.templates.h b/include/deal.II/numerics/error_estimator.templates.h
index d490ffa785..78594efeb6 100644
--- a/include/deal.II/numerics/error_estimator.templates.h
+++ b/include/deal.II/numerics/error_estimator.templates.h
@@ -942,8 +942,7 @@ namespace internal
     local_face_integrals.clear();
 
     // loop over all faces of this cell
-    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
-         ++face_no)
+    for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
       {
         const typename DoFHandlerType::face_iterator face = cell->face(face_no);
 
diff --git a/include/deal.II/numerics/matrix_creator.templates.h b/include/deal.II/numerics/matrix_creator.templates.h
index 4bd279065b..b24a74b452 100644
--- a/include/deal.II/numerics/matrix_creator.templates.h
+++ b/include/deal.II/numerics/matrix_creator.templates.h
@@ -1056,8 +1056,7 @@ namespace MatrixCreator
       copy_data.cell_matrix.clear();
       copy_data.cell_vector.clear();
 
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const unsigned int face : GeometryInfo<dim>::face_indices())
         // check if this face is on that part of the boundary we are
         // interested in
         if (boundary_functions.find(cell->face(face)->boundary_id()) !=
@@ -1274,8 +1273,7 @@ namespace MatrixCreator
       // inefficient, so we copy the dofs into a set, which enables binary
       // searches.
       unsigned int pos(0);
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const unsigned int face : GeometryInfo<dim>::face_indices())
         {
           // check if this face is on that part of
           // the boundary we are interested in
@@ -1499,8 +1497,7 @@ namespace MatrixCreator
       copy_data.cell_vector.clear();
 
 
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const unsigned int face : GeometryInfo<dim>::face_indices())
         // check if this face is on that part of
         // the boundary we are interested in
         if (boundary_functions.find(cell->face(face)->boundary_id()) !=
@@ -1733,8 +1730,7 @@ namespace MatrixCreator
       // inefficient, so we copy the dofs into a set, which enables binary
       // searches.
       unsigned int pos(0);
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const unsigned int face : GeometryInfo<dim>::face_indices())
         {
           // check if this face is on that part of
           // the boundary we are interested in
diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc
index 4c85b91b3e..2a0372cdee 100644
--- a/source/dofs/dof_tools_constraints.cc
+++ b/source/dofs/dof_tools_constraints.cc
@@ -663,8 +663,7 @@ namespace DoFTools
           if (cell->is_artificial())
             continue;
 
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (const unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->has_children())
               {
                 // in any case, faces can have at most two active fe indices,
@@ -812,8 +811,7 @@ namespace DoFTools
           if (cell->is_artificial())
             continue;
 
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (const unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->has_children())
               {
                 // first of all, make sure that we treat a case which is
@@ -1085,8 +1083,7 @@ namespace DoFTools
           if (cell->is_artificial())
             continue;
 
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (const unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->has_children())
               {
                 // first of all, make sure that we treat a case which is
diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc
index a4c3ddc967..62c2f5ebf5 100644
--- a/source/dofs/dof_tools_sparsity.cc
+++ b/source/dofs/dof_tools_sparsity.cc
@@ -790,7 +790,7 @@ namespace DoFTools
         Table<2, bool> support_on_face(fe.dofs_per_cell,
                                        GeometryInfo<dim>::faces_per_cell);
         for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (const unsigned int f : GeometryInfo<dim>::face_indices())
             support_on_face(i, f) = fe.has_support_on_face(i, f);
 
         // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc
index 2ff7ad0373..a41a54ade3 100644
--- a/source/fe/fe_enriched.cc
+++ b/source/fe/fe_enriched.cc
@@ -1293,8 +1293,7 @@ namespace ColorEnriched
         {
           const unsigned int           fe_index = cell->active_fe_index();
           const std::set<unsigned int> fe_set   = fe_sets.at(fe_index);
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (const unsigned int face : GeometryInfo<dim>::face_indices())
             {
               // cell shouldn't be at the boundary and
               // neighboring cell is not already visited (to avoid visiting
diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc
index 158ce585cd..b105aa9f4e 100644
--- a/source/fe/fe_nedelec.cc
+++ b/source/fe/fe_nedelec.cc
@@ -442,8 +442,7 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
       const Quadrature<dim> &face_quadrature =
         QProjector<dim>::project_to_all_faces(reference_face_quadrature);
 
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const unsigned int face : GeometryInfo<dim>::face_indices())
         for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
           {
             this->generalized_support_points[face * n_face_points + q_point +
diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc
index d88b1a6367..23de12715e 100644
--- a/source/fe/mapping_fe_field.cc
+++ b/source/fe/mapping_fe_field.cc
@@ -589,8 +589,7 @@ MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::compute_face_data(
               // ensure a counterclockwise
               // orientation of tangentials
               static const int tangential_orientation[4] = {-1, 1, 1, -1};
-              for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
-                   ++i)
+              for (const unsigned int i : GeometryInfo<dim>::face_indices())
                 {
                   Tensor<1, dim> tang;
                   tang[1 - i / 2] = tangential_orientation[i];
@@ -601,8 +600,7 @@ MappingFEField<dim, spacedim, VectorType, DoFHandlerType>::compute_face_data(
             }
           else if (dim == 3)
             {
-              for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
-                   ++i)
+              for (const unsigned int i : GeometryInfo<dim>::face_indices())
                 {
                   Tensor<1, dim> tang1, tang2;
 
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 9e19f3f724..e496d47163 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -3798,8 +3798,7 @@ namespace GridGenerator
     // The index of the first cell of the leg.
     unsigned int cell_index = 1;
     // The legs of the cross
-    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-         ++face)
+    for (const unsigned int face : GeometryInfo<dim>::face_indices())
       {
         const unsigned int oface = GeometryInfo<dim>::opposite_face[face];
         const unsigned int dir = GeometryInfo<dim>::unit_normal_direction[face];
diff --git a/source/grid/grid_out.cc b/source/grid/grid_out.cc
index 8757b4c723..51c70a635e 100644
--- a/source/grid/grid_out.cc
+++ b/source/grid/grid_out.cc
@@ -956,8 +956,7 @@ GridOut::write_dx(const Triangulation<dim, spacedim> &tria,
               << n_faces << " data follows" << '\n';
           for (const auto &cell : tria.active_cell_iterators())
             {
-              for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
-                   ++f)
+              for (const unsigned int f : GeometryInfo<dim>::face_indices())
                 out << ' ' << cell->face(f)->measure();
               out << '\n';
             }
@@ -970,8 +969,7 @@ GridOut::write_dx(const Triangulation<dim, spacedim> &tria,
               << n_faces << " data follows" << '\n';
           for (const auto &cell : tria.active_cell_iterators())
             {
-              for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
-                   ++f)
+              for (const unsigned int f : GeometryInfo<dim>::face_indices())
                 out << ' ' << cell->face(f)->diameter();
               out << '\n';
             }
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index dbf5adf3f5..e952a81502 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -103,8 +103,7 @@ namespace GridTools
     const typename Triangulation<dim, spacedim>::active_cell_iterator endc =
       tria.end();
     for (; cell != endc; ++cell)
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const unsigned int face : GeometryInfo<dim>::face_indices())
         if (cell->face(face)->at_boundary())
           for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
                ++i)
@@ -4138,7 +4137,7 @@ namespace GridTools
         continue_refinement = false;
 
         for (const auto &cell : tria.active_cell_iterators())
-          for (unsigned int j = 0; j < GeometryInfo<dim>::faces_per_cell; j++)
+          for (const unsigned int j : GeometryInfo<dim>::face_indices())
             if (cell->at_boundary(j) == false &&
                 cell->neighbor(j)->has_children())
               {
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 2333a0b342..bb39d652c0 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -10440,8 +10440,7 @@ Triangulation<dim, spacedim>::get_boundary_ids() const
       std::set<types::boundary_id> b_ids;
       active_cell_iterator         cell = begin_active();
       for (; cell != end(); ++cell)
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           if (cell->at_boundary(face))
             b_ids.insert(cell->face(face)->boundary_id());
       std::vector<types::boundary_id> boundary_ids(b_ids.begin(), b_ids.end());
@@ -10461,8 +10460,7 @@ Triangulation<dim, spacedim>::get_manifold_ids() const
     {
       m_ids.insert(cell->manifold_id());
       if (dim > 1)
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           if (cell->at_boundary(face))
             m_ids.insert(cell->face(face)->manifold_id());
     }
@@ -10641,10 +10639,8 @@ Triangulation<dim, spacedim>::create_triangulation(
           case 1:
             {
               bool values[][2] = {{false, true}, {true, false}};
-              for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
-                   ++i)
-                for (unsigned int j = 0; j < GeometryInfo<dim>::faces_per_cell;
-                     ++j)
+              for (const unsigned int i : GeometryInfo<dim>::face_indices())
+                for (const unsigned int j : GeometryInfo<dim>::face_indices())
                   correct(i, j) = (values[i][j]);
               break;
             }
@@ -10654,10 +10650,8 @@ Triangulation<dim, spacedim>::create_triangulation(
                                   {true, false, false, true},
                                   {true, false, false, true},
                                   {false, true, true, false}};
-              for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
-                   ++i)
-                for (unsigned int j = 0; j < GeometryInfo<dim>::faces_per_cell;
-                     ++j)
+              for (const unsigned int i : GeometryInfo<dim>::face_indices())
+                for (const unsigned int j : GeometryInfo<dim>::face_indices())
                   correct(i, j) = (values[i][j]);
               break;
             }
@@ -10680,8 +10674,7 @@ Triangulation<dim, spacedim>::create_triangulation(
                cell != this_round.end();
                ++cell)
             {
-              for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
-                   ++i)
+              for (const unsigned int i : GeometryInfo<dim>::face_indices())
                 {
                   if (!((*cell)->face(i)->at_boundary()))
                     {
@@ -13984,8 +13977,7 @@ namespace
         for (unsigned int c = 0; c < cell->n_children(); ++c)
           cell->child(c)->clear_coarsen_flag();
 
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           if (!cell->at_boundary(face) &&
               (!cell->neighbor(face)->is_active()) &&
               (cell_will_be_coarsened(cell->neighbor(face))))
@@ -14028,8 +14020,7 @@ namespace
       {
         // use first algorithm
         unsigned int refined_neighbors = 0, unrefined_neighbors = 0;
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           if (!cell->at_boundary(face))
             {
               if (face_will_be_refined_by_neighbor(cell, face))
@@ -14668,8 +14659,7 @@ Triangulation<dim, spacedim>::prepare_coarsening_and_refinement()
             if (cell->refine_flag_set())
               {
                 // loop over neighbors of cell
-                for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
-                     ++i)
+                for (const unsigned int i : GeometryInfo<dim>::face_indices())
                   {
                     // only do something if the face is not at the
                     // boundary and if the face will be refined with
diff --git a/source/grid/tria_description.cc b/source/grid/tria_description.cc
index 0905de3017..4f843c349d 100644
--- a/source/grid/tria_description.cc
+++ b/source/grid/tria_description.cc
@@ -242,8 +242,7 @@ namespace TriangulationDescription
                     numbers::invalid_unsigned_int;
 
                 // c) save boundary_ids of each face of this cell
-                for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
-                     ++f)
+                for (const unsigned int f : GeometryInfo<dim>::face_indices())
                   {
                     types::boundary_id boundary_ind =
                       cell->face(f)->boundary_id();
diff --git a/source/multigrid/mg_tools.cc b/source/multigrid/mg_tools.cc
index dbd96edf6c..c1d515835e 100644
--- a/source/multigrid/mg_tools.cc
+++ b/source/multigrid/mg_tools.cc
@@ -208,8 +208,7 @@ namespace MGTools
         // and add the missing
         // contribution as well as the
         // flux contributions.
-        for (unsigned int iface = 0; iface < GeometryInfo<dim>::faces_per_cell;
-             ++iface)
+        for (const unsigned int iface : GeometryInfo<dim>::face_indices())
           {
             bool level_boundary = cell->at_boundary(iface);
             typename DoFHandler<dim, spacedim>::cell_iterator neighbor;
@@ -460,8 +459,7 @@ namespace MGTools
         // and add the missing
         // contribution as well as the
         // flux contributions.
-        for (unsigned int iface = 0; iface < GeometryInfo<dim>::faces_per_cell;
-             ++iface)
+        for (const unsigned int iface : GeometryInfo<dim>::face_indices())
           {
             bool level_boundary = cell->at_boundary(iface);
             typename DoFHandler<dim, spacedim>::cell_iterator neighbor;
@@ -627,8 +625,7 @@ namespace MGTools
             sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
 
         // Loop over all interior neighbors
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           {
             bool use_face = false;
             if ((!cell->at_boundary(face)) &&
@@ -705,8 +702,7 @@ namespace MGTools
 
         cell->get_mg_dof_indices(dofs_on_this_cell);
         // Loop over all interior neighbors
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           {
             // Neighbor is coarser
             bool use_face = false;
@@ -814,8 +810,7 @@ namespace MGTools
               sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
 
         // Loop over all interior neighbors
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           {
             typename DoFHandler<dim, spacedim>::face_iterator cell_face =
               cell->face(face);
@@ -984,8 +979,7 @@ namespace MGTools
 
         cell->get_mg_dof_indices(dofs_on_this_cell);
         // Loop over all interior neighbors
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
           {
             // Neighbor is coarser
             bool use_face = false;
-- 
2.39.5