From 8aa28dd2f179d419694111e09b056ca26b782a7f Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 21 Jan 2020 20:00:08 -0700
Subject: [PATCH] Use GeometryInfo::face_indices().

---
 include/deal.II/base/geometry_info.h          | 31 ++++++++++++++-
 include/deal.II/dofs/dof_accessor.templates.h |  2 +-
 .../deal.II/grid/tria_accessor.templates.h    |  2 +-
 .../deal.II/matrix_free/face_setup_internal.h |  2 +-
 .../matrix_free/matrix_free.templates.h       |  2 +-
 .../matrix_free/shape_info.templates.h        |  4 +-
 include/deal.II/meshworker/dof_info.h         |  8 ++--
 include/deal.II/meshworker/loop.h             |  2 +-
 .../deal.II/multigrid/mg_constrained_dofs.h   |  2 +-
 .../deal.II/numerics/vector_tools.templates.h | 28 +++++---------
 source/base/data_out_base.cc                  | 10 ++---
 source/base/polynomials_bernardi_raugel.cc    |  2 +-
 source/base/quadrature_lib.cc                 |  2 +-
 source/distributed/tria.cc                    |  4 +-
 source/dofs/dof_tools.cc                      |  2 +-
 source/fe/fe_abf.cc                           |  6 +--
 source/fe/fe_bdm.cc                           |  2 +-
 source/fe/fe_bernardi_raugel.cc               |  2 +-
 source/fe/fe_face.cc                          |  2 +-
 source/fe/fe_raviart_thomas.cc                |  4 +-
 source/fe/mapping_q_generic.cc                |  2 +-
 source/grid/grid_generator.cc                 | 38 +++++++++----------
 source/grid/grid_in.cc                        |  2 +-
 source/grid/grid_out.cc                       |  6 +--
 source/grid/grid_tools.cc                     | 18 ++++-----
 source/grid/grid_tools_dof_handlers.cc        |  2 +-
 source/grid/tria.cc                           | 20 +++++-----
 source/hp/dof_handler.cc                      |  4 +-
 source/multigrid/mg_tools.cc                  |  4 +-
 source/opencascade/utilities.cc               |  2 +-
 30 files changed, 119 insertions(+), 98 deletions(-)

diff --git a/include/deal.II/base/geometry_info.h b/include/deal.II/base/geometry_info.h
index 76e752e8e2..48993f08de 100644
--- a/include/deal.II/base/geometry_info.h
+++ b/include/deal.II/base/geometry_info.h
@@ -24,6 +24,7 @@
 
 #include <boost/range/irange.hpp>
 
+#include <array>
 #include <cstdint>
 
 
@@ -1204,6 +1205,25 @@ struct GeometryInfo<0>
    */
   static constexpr unsigned int faces_per_cell = 0;
 
+  /**
+   * Return an object that can be thought of as an array containing all
+   * indices from zero to `faces_per_cell`. This allows to write code
+   * 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)
+   *       if (cell->face(face_index)->at_boundary())
+   *         ... do something ...
+   * @endcode
+   * Here, we are looping over all faces of all cells, with `face_index`
+   * taking on all valid indices.
+   *
+   * Of course, since this class is for the case `dim==0`, the
+   * returned object is actually an empty array.
+   */
+  static std::array<unsigned int, 0>
+  face_indices();
+
   /**
    * Maximum number of children of a refined face, i.e. the number of children
    * of an isotropically refined face.
@@ -2686,14 +2706,23 @@ GeometryInfo<3>::unit_cell_vertex(const unsigned int vertex)
 
 
 
+inline std::array<unsigned int, 0>
+GeometryInfo<0>::face_indices()
+{
+  return {};
+}
+
+
+
 template <int dim>
-boost::integer_range<unsigned int>
+inline boost::integer_range<unsigned int>
 GeometryInfo<dim>::face_indices()
 {
   return boost::irange(0U, faces_per_cell);
 }
 
 
+
 template <int dim>
 inline Point<dim>
 GeometryInfo<dim>::unit_cell_vertex(const unsigned int)
diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h
index 4766a25089..269f454547 100644
--- a/include/deal.II/dofs/dof_accessor.templates.h
+++ b/include/deal.II/dofs/dof_accessor.templates.h
@@ -3795,7 +3795,7 @@ DoFCellAccessor<DoFHandlerType, level_dof_access>::face_iterators() const
     face_iterators;
 
   const unsigned int dim = DoFHandlerType::dimension;
-  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+  for (unsigned int i : GeometryInfo<dim>::face_indices())
     face_iterators[i] =
       dealii::internal::DoFCellAccessorImplementation::get_face(
         *this, i, std::integral_constant<int, dim>());
diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h
index a9fd164534..aabb4fc29d 100644
--- a/include/deal.II/grid/tria_accessor.templates.h
+++ b/include/deal.II/grid/tria_accessor.templates.h
@@ -3232,7 +3232,7 @@ CellAccessor<dim, spacedim>::face_iterators() const
              GeometryInfo<dim>::faces_per_cell>
     face_iterators;
 
-  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+  for (unsigned int i : GeometryInfo<dim>::face_indices())
     face_iterators[i] =
       dealii::internal::CellAccessorImplementation::get_face(*this, i);
 
diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index f475a4d44a..2cfd03a943 100644
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -579,7 +579,7 @@ namespace internal
             &triangulation, cell_levels[i].first, cell_levels[i].second);
           if (use_active_cells)
             Assert(dcell->is_active(), ExcNotImplemented());
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             {
               if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
                 face_is_owned[dcell->face(f)->index()] =
diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index dcfa56aa04..5297c499a3 100644
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -1855,7 +1855,7 @@ namespace internal
         new_indices.clear();
         typename dealii::Triangulation<dim>::cell_iterator dcell(
           &tria, cell_level_index[cell].first, cell_level_index[cell].second);
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           {
             // Only inner faces couple different cells
             if (dcell->at_boundary(f) == false &&
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 3be0eea9b9..c89ba34f22 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -426,7 +426,7 @@ namespace internal
         {
           face_to_cell_index_nodal.reinit(GeometryInfo<dim>::faces_per_cell,
                                           dofs_per_component_on_face);
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             {
               const unsigned int direction = f / 2;
               const unsigned int stride =
@@ -459,7 +459,7 @@ namespace internal
         {
           face_to_cell_index_hermite.reinit(GeometryInfo<dim>::faces_per_cell,
                                             2 * dofs_per_component_on_face);
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             {
               const unsigned int direction = f / 2;
               const unsigned int stride =
diff --git a/include/deal.II/meshworker/dof_info.h b/include/deal.II/meshworker/dof_info.h
index 2f0222ee54..177c7a20e0 100644
--- a/include/deal.II/meshworker/dof_info.h
+++ b/include/deal.II/meshworker/dof_info.h
@@ -427,7 +427,7 @@ namespace MeshWorker
     : cell(seed)
     , cell_valid(true)
   {
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       {
         exterior[i]                = seed;
         interior[i]                = seed;
@@ -443,7 +443,7 @@ namespace MeshWorker
     : cell(other.cell)
     , cell_valid(other.cell_valid)
   {
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       {
         exterior[i]                = other.exterior[i];
         interior[i]                = other.interior[i];
@@ -458,7 +458,7 @@ namespace MeshWorker
   DoFInfoBox<dim, DOFINFO>::reset()
   {
     cell_valid = false;
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       {
         interior_face_available[i] = false;
         exterior_face_available[i] = false;
@@ -475,7 +475,7 @@ namespace MeshWorker
       return;
 
     assembler.assemble(cell);
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       {
         // Only do something if data available
         if (interior_face_available[i])
diff --git a/include/deal.II/meshworker/loop.h b/include/deal.II/meshworker/loop.h
index 69eff70e64..390a01552d 100644
--- a/include/deal.II/meshworker/loop.h
+++ b/include/deal.II/meshworker/loop.h
@@ -461,7 +461,7 @@ namespace MeshWorker
     DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
 
     assembler.initialize_info(dof_info.cell, false);
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       {
         assembler.initialize_info(dof_info.interior[i], true);
         assembler.initialize_info(dof_info.exterior[i], true);
diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h
index 14348002ef..d6641039a1 100644
--- a/include/deal.II/multigrid/mg_constrained_dofs.h
+++ b/include/deal.II/multigrid/mg_constrained_dofs.h
@@ -277,7 +277,7 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim, spacedim> &dof)
       for (; cell != endc; ++cell)
         if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
           {
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               if (cell->has_periodic_neighbor(f) &&
                   cell->periodic_neighbor(f)->level() == cell->level())
                 {
diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 3a97314fa8..3520aadf0f 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -745,7 +745,7 @@ namespace VectorTools
         endc = dof_handler.end();
       std::vector<types::global_dof_index> face_dof_indices;
       for (; cell != endc; ++cell)
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->at_boundary(f))
             {
               face_dof_indices.resize(cell->get_fe().dofs_per_face);
@@ -2504,8 +2504,7 @@ namespace VectorTools
         std::vector<double> rhs_values(n_q_points);
 
         for (; cell != endc; ++cell)
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->at_boundary() &&
                 (boundary_ids.empty() ||
                  (boundary_ids.find(cell->face(face)->boundary_id()) !=
@@ -2536,8 +2535,7 @@ namespace VectorTools
                                                Vector<double>(n_components));
 
         for (; cell != endc; ++cell)
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->at_boundary() &&
                 (boundary_ids.empty() ||
                  (boundary_ids.find(cell->face(face)->boundary_id()) !=
@@ -2651,8 +2649,7 @@ namespace VectorTools
         std::vector<double> rhs_values;
 
         for (; cell != endc; ++cell)
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->at_boundary() &&
                 (boundary_ids.empty() ||
                  (boundary_ids.find(cell->face(face)->boundary_id()) !=
@@ -2690,8 +2687,7 @@ namespace VectorTools
         std::vector<Vector<double>> rhs_values;
 
         for (; cell != endc; ++cell)
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
+          for (unsigned int face : GeometryInfo<dim>::face_indices())
             if (cell->face(face)->at_boundary() &&
                 (boundary_ids.empty() ||
                  (boundary_ids.find(cell->face(face)->boundary_id()) !=
@@ -3542,7 +3538,7 @@ namespace VectorTools
                  cell = dof.begin_active();
                cell != dof.end();
                ++cell)
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               {
                 if (cell->at_boundary(f))
                   {
@@ -4895,8 +4891,7 @@ namespace VectorTools
     const hp::MappingCollection<dim> mapping_collection(mapping);
     hp::QCollection<dim>             face_quadrature_collection;
 
-    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-         ++face)
+    for (unsigned int face : GeometryInfo<dim>::face_indices())
       face_quadrature_collection.push_back(
         QProjector<dim>::project_to_face(reference_face_quadrature, face));
 
@@ -5131,8 +5126,7 @@ namespace VectorTools
         const QGauss<dim - 1> reference_face_quadrature(
           2 * fe_collection[i].degree);
 
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (unsigned int face : GeometryInfo<dim>::face_indices())
           face_quadrature_collection.push_back(
             QProjector<dim>::project_to_face(reference_face_quadrature, face));
       }
@@ -6656,8 +6650,7 @@ namespace VectorTools
     const hp::MappingCollection<dim> mapping_collection(mapping);
     hp::QCollection<dim>             quadrature_collection;
 
-    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-         ++face)
+    for (unsigned int face : GeometryInfo<dim>::face_indices())
       quadrature_collection.push_back(
         QProjector<dim>::project_to_face(face_quadrature, face));
 
@@ -6828,8 +6821,7 @@ namespace VectorTools
 
         face_quadrature_collection.push_back(quadrature);
 
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (unsigned int face : GeometryInfo<dim>::face_indices())
           quadrature_collection.push_back(
             QProjector<dim>::project_to_face(quadrature, face));
       }
diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc
index 32f832dfb4..8612c0ea37 100644
--- a/source/base/data_out_base.cc
+++ b/source/base/data_out_base.cc
@@ -1757,7 +1757,7 @@ namespace DataOutBase
   // all the other data has a constructor of its own, except for the "neighbors"
   // field, which we set to invalid values.
   {
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       neighbors[i] = no_neighbor;
 
     AssertIndexRange(dim, spacedim + 1);
@@ -1776,7 +1776,7 @@ namespace DataOutBase
       if (vertices[i].distance(patch.vertices[i]) > epsilon)
         return false;
 
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       if (neighbors[i] != patch.neighbors[i])
         return false;
 
@@ -8856,7 +8856,7 @@ DataOutReader<dim, spacedim>::merge(const DataOutReader<dim, spacedim> &source)
 
   // adjust patch neighbors
   for (unsigned int i = old_n_patches; i < patches.size(); ++i)
-    for (unsigned int n = 0; n < GeometryInfo<dim>::faces_per_cell; ++n)
+    for (unsigned int n : GeometryInfo<dim>::face_indices())
       if (patches[i].neighbors[n] !=
           dealii::DataOutBase::Patch<dim, spacedim>::no_neighbor)
         patches[i].neighbors[n] += old_n_patches;
@@ -9071,7 +9071,7 @@ namespace DataOutBase
       out << patch.vertices[GeometryInfo<dim>::ucd_to_deal[i]] << ' ';
     out << '\n';
 
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       out << patch.neighbors[i] << ' ';
     out << '\n';
 
@@ -9119,7 +9119,7 @@ namespace DataOutBase
     for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
       in >> patch.vertices[GeometryInfo<dim>::ucd_to_deal[i]];
 
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       in >> patch.neighbors[i];
 
     in >> patch.patch_index >> patch.n_subdivisions;
diff --git a/source/base/polynomials_bernardi_raugel.cc b/source/base/polynomials_bernardi_raugel.cc
index a42ac870e4..38e822bc7a 100644
--- a/source/base/polynomials_bernardi_raugel.cc
+++ b/source/base/polynomials_bernardi_raugel.cc
@@ -120,7 +120,7 @@ PolynomialsBernardiRaugel<dim>::evaluate(
   // Normal vectors point in the +x, +y, and +z directions for
   // consistent orientation across edges
   std::vector<Tensor<1, dim>> normals;
-  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+  for (unsigned int i : GeometryInfo<dim>::face_indices())
     {
       Tensor<1, dim> normal;
       normal[i / 2] = 1;
diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc
index 4f5223fdf6..875f519e65 100644
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -1328,7 +1328,7 @@ QSplit<dim>::QSplit(const QSimplex<dim> &base, const Point<dim> &split_point)
   // face. In dimension three, we need to split the face in two triangles, so
   // we use once the first dim vertices of each face, and the second time the
   // the dim vertices of each face starting from 1.
-  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+  for (auto f : GeometryInfo<dim>::face_indices())
     for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
       {
         for (unsigned int i = 0; i < dim; ++i)
diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc
index 02d22bc7db..2477b09b0c 100644
--- a/source/distributed/tria.cc
+++ b/source/distributed/tria.cc
@@ -160,7 +160,7 @@ namespace
 
         // neighborship information. if a cell is at a boundary, then enter
         // the index of the cell itself here
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->face(f)->at_boundary() == false)
             connectivity
               ->tree_to_tree[index * GeometryInfo<dim>::faces_per_cell + f] =
@@ -172,7 +172,7 @@ namespace
 
         // fill tree_to_face, which is essentially neighbor_to_neighbor;
         // however, we have to remap the resulting face number as well
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->face(f)->at_boundary() == false)
             {
               switch (dim)
diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc
index 6a1d7c84dd..134ef0719a 100644
--- a/source/dofs/dof_tools.cc
+++ b/source/dofs/dof_tools.cc
@@ -1038,7 +1038,7 @@ namespace DoFTools
 
         for (const auto &cell : dof_handler.active_cell_iterators())
           if (!cell->is_artificial())
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               {
                 const typename dealii::DoFHandler<dim, spacedim>::face_iterator
                   face = cell->face(f);
diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc
index f6e5c01382..c343605fe0 100644
--- a/source/fe/fe_abf.cc
+++ b/source/fe/fe_abf.cc
@@ -333,7 +333,7 @@ FE_ABF<dim>::initialize_restriction()
   const unsigned int n_face_points = q_base.size();
   // First, compute interpolation on
   // subfaces
-  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+  for (unsigned int face : GeometryInfo<dim>::face_indices())
     {
       // The shape functions of the
       // child cell are evaluated
@@ -546,7 +546,7 @@ FE_ABF<dim>::convert_generalized_support_point_values_to_dof_values(
   std::fill(nodal_values.begin(), nodal_values.end(), 0.);
 
   const unsigned int n_face_points = boundary_weights.size(0);
-  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+  for (unsigned int face : GeometryInfo<dim>::face_indices())
     for (unsigned int k = 0; k < n_face_points; ++k)
       for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
         {
@@ -580,7 +580,7 @@ FE_ABF<dim>::convert_generalized_support_point_values_to_dof_values(
           support_point_values[k + start_cell_points][d];
 
   // Face integral of ABF terms
-  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+  for (unsigned int face : GeometryInfo<dim>::face_indices())
     {
       const double n_orient = GeometryInfo<dim>::unit_normal_orientation[face];
       for (unsigned int fp = 0; fp < n_face_points; ++fp)
diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc
index 26b747877b..eb66ffec30 100644
--- a/source/fe/fe_bdm.cc
+++ b/source/fe/fe_bdm.cc
@@ -149,7 +149,7 @@ FE_BDM<dim>::convert_generalized_support_point_values_to_dof_values(
   unsigned int dbase = 0;
   // The index of the first generalized support point on this face or the cell
   unsigned int pbase = 0;
-  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+  for (auto f : GeometryInfo<dim>::face_indices())
     {
       // Old version with no moments in 2D. See comment below in
       // initialize_support_points()
diff --git a/source/fe/fe_bernardi_raugel.cc b/source/fe/fe_bernardi_raugel.cc
index 4408442b88..97113e5d36 100644
--- a/source/fe/fe_bernardi_raugel.cc
+++ b/source/fe/fe_bernardi_raugel.cc
@@ -101,7 +101,7 @@ FE_BernardiRaugel<dim>::convert_generalized_support_point_values_to_dof_values(
          ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell));
 
   std::vector<Tensor<1, dim>> normals;
-  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+  for (unsigned int i : GeometryInfo<dim>::face_indices())
     {
       Tensor<1, dim> normal;
       normal[i / 2] = 1;
diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc
index 02db2fb446..e548282a32 100644
--- a/source/fe/fe_face.cc
+++ b/source/fe/fe_face.cc
@@ -979,7 +979,7 @@ std::pair<Table<2, bool>, std::vector<unsigned int>>
 FE_FaceP<dim, spacedim>::get_constant_modes() const
 {
   Table<2, bool> constant_modes(1, this->dofs_per_cell);
-  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+  for (unsigned int face : GeometryInfo<dim>::face_indices())
     constant_modes(0, face * this->dofs_per_face) = true;
   return std::pair<Table<2, bool>, std::vector<unsigned int>>(
     constant_modes, std::vector<unsigned int>(1, 0));
diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc
index 7fe49b26f9..3702b2ce91 100644
--- a/source/fe/fe_raviart_thomas.cc
+++ b/source/fe/fe_raviart_thomas.cc
@@ -271,7 +271,7 @@ FE_RaviartThomas<dim>::initialize_restriction()
   const unsigned int n_face_points = q_base.size();
   // First, compute interpolation on
   // subfaces
-  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+  for (unsigned int face : GeometryInfo<dim>::face_indices())
     {
       // The shape functions of the
       // child cell are evaluated
@@ -497,7 +497,7 @@ FE_RaviartThomas<dim>::convert_generalized_support_point_values_to_dof_values(
   std::fill(nodal_values.begin(), nodal_values.end(), 0.);
 
   const unsigned int n_face_points = boundary_weights.size(0);
-  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+  for (unsigned int face : GeometryInfo<dim>::face_indices())
     for (unsigned int k = 0; k < n_face_points; ++k)
       for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
         {
diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index d5d911c8b3..ddc4115f7b 100644
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -4087,7 +4087,7 @@ MappingQGeneric<dim, spacedim>::compute_mapping_support_points(
           dynamic_cast<const TransfiniteInterpolationManifold<dim, spacedim> *>(
             &cell->get_manifold()) == nullptr)
         {
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             if (&cell->face(f)->get_manifold() != &cell->get_manifold())
               all_manifold_ids_are_equal = false;
 
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 1afe467000..a09b1c086c 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -896,7 +896,7 @@ namespace GridGenerator
         static void set_boundary_ids(Triangulation<2> &tria)
         {
           for (auto cell : tria.active_cell_iterators())
-            for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
+            for (unsigned int f : GeometryInfo<2>::face_indices())
               {
                 if (cell->face(f)->at_boundary() == false)
                   continue;
@@ -1235,7 +1235,7 @@ namespace GridGenerator
           // simple task
           const typename Triangulation<dim, spacedim>::cell_iterator cell =
             tria.begin();
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             cell->face(f)->set_boundary_id(f);
         }
     }
@@ -1443,7 +1443,7 @@ namespace GridGenerator
       Triangulation<3>::cell_iterator cell = tria.begin();
 
       for (; cell != tria.end(); ++cell)
-        for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f)
+        for (unsigned int f : GeometryInfo<3>::face_indices())
           {
             if (!cell->face(f)->at_boundary())
               continue;
@@ -2020,7 +2020,7 @@ namespace GridGenerator
     for (auto &cell : tria.cell_iterators())
       {
         // identify faces on torus surface and set manifold to 1
-        for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f)
+        for (unsigned int f : GeometryInfo<3>::face_indices())
           {
             // faces 4 and 5 are those with normal vector aligned with torus
             // centerline
@@ -2944,7 +2944,7 @@ namespace GridGenerator
         for (; cell != endc; ++cell)
           {
             Point<2> cell_center = cell->center();
-            for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
+            for (unsigned int f : GeometryInfo<2>::face_indices())
               if (cell->face(f)->boundary_id() == 0)
                 {
                   Point<2> face_center = cell->face(f)->center();
@@ -3055,7 +3055,7 @@ namespace GridGenerator
         for (; cell != endc; ++cell)
           {
             Point<dim> cell_center = cell->center();
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               if (cell->face(f)->boundary_id() == 0)
                 {
                   Point<dim> face_center = cell->face(f)->center();
@@ -3778,7 +3778,7 @@ namespace GridGenerator
 
     std::vector<Point<spacedim>> points;
     unsigned int                 n_cells = 1;
-    for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    for (unsigned int i : GeometryInfo<dim>::face_indices())
       n_cells += sizes[i];
 
     std::vector<CellData<dim>> cells(n_cells);
@@ -4456,7 +4456,7 @@ namespace GridGenerator
 
     while (cell != end)
       {
-        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+        for (unsigned int i : GeometryInfo<dim>::face_indices())
           {
             if (cell->face(i)->boundary_id() ==
                 numbers::internal_face_boundary_id)
@@ -4522,7 +4522,7 @@ namespace GridGenerator
 
     while (cell != end)
       {
-        for (unsigned int i = 0; i < GeometryInfo<2>::faces_per_cell; ++i)
+        for (unsigned int i : GeometryInfo<2>::face_indices())
           {
             if (cell->face(i)->boundary_id() ==
                 numbers::internal_face_boundary_id)
@@ -5100,7 +5100,7 @@ namespace GridGenerator
     tria.set_all_manifold_ids_on_boundary(0);
 
     for (; cell != end; ++cell)
-      for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
+      for (unsigned int i : GeometryInfo<3>::face_indices())
         if (cell->at_boundary(i))
           {
             if (cell->face(i)->center()(0) > half_length - 1.e-5)
@@ -5193,7 +5193,7 @@ namespace GridGenerator
     tria.set_all_manifold_ids_on_boundary(0);
     while (cell != end)
       {
-        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+        for (unsigned int i : GeometryInfo<dim>::face_indices())
           {
             if (cell->face(i)->boundary_id() ==
                 numbers::internal_face_boundary_id)
@@ -5301,7 +5301,7 @@ namespace GridGenerator
     // set to one
     while (cell != end)
       {
-        for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
+        for (unsigned int i : GeometryInfo<3>::face_indices())
           {
             if (!cell->at_boundary(i))
               continue;
@@ -5551,7 +5551,7 @@ namespace GridGenerator
         // to all faces in a first step.
         Triangulation<3>::cell_iterator cell = tria.begin();
         for (; cell != tria.end(); ++cell)
-          for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
+          for (unsigned int i : GeometryInfo<3>::face_indices())
             if (cell->at_boundary(i))
               cell->face(i)->set_all_boundary_ids(2);
 
@@ -5560,7 +5560,7 @@ namespace GridGenerator
         // boundary. Then decide whether the center is nearer to the inner
         // or outer boundary to set the correct boundary id.
         for (cell = tria.begin(); cell != tria.end(); ++cell)
-          for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
+          for (unsigned int i : GeometryInfo<3>::face_indices())
             if (cell->at_boundary(i))
               {
                 const Triangulation<3>::face_iterator face = cell->face(i);
@@ -6256,7 +6256,7 @@ namespace GridGenerator
     std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
     for (; cell != endc; ++cell)
       {
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->face(f)->at_boundary())
             {
               for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
@@ -6295,7 +6295,7 @@ namespace GridGenerator
     cell       = triangulation.begin_active();
     for (; cell != endc; ++cell)
       {
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->face(f)->at_boundary())
             {
               double dx = cell->face(f)->center()(0) - center(0);
@@ -6474,7 +6474,7 @@ namespace GridGenerator
     std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
     for (; cell != endc; ++cell)
       {
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->face(f)->at_boundary())
             {
               for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
@@ -6517,7 +6517,7 @@ namespace GridGenerator
     cell       = triangulation.begin_active();
     for (; cell != endc; ++cell)
       {
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->face(f)->at_boundary())
             {
               double dx = cell->face(f)->center()(0);
@@ -6750,7 +6750,7 @@ namespace GridGenerator
            volume_mesh.begin(0);
          cell != volume_mesh.end(0);
          ++cell)
-      for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+      for (unsigned int i : GeometryInfo<dim>::face_indices())
         {
           const typename MeshType<dim, spacedim>::face_iterator face =
             cell->face(i);
diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc
index bf259751ff..547a407077 100644
--- a/source/grid/grid_in.cc
+++ b/source/grid/grid_in.cc
@@ -67,7 +67,7 @@ namespace
              triangulation.begin_active();
            cell != triangulation.end();
            ++cell)
-        for (unsigned int f = 0; f < GeometryInfo<1>::faces_per_cell; ++f)
+        for (unsigned int f : GeometryInfo<1>::face_indices())
           if (boundary_ids.find(cell->vertex_index(f)) != boundary_ids.end())
             {
               AssertThrow(
diff --git a/source/grid/grid_out.cc b/source/grid/grid_out.cc
index 25f4d54259..8757b4c723 100644
--- a/source/grid/grid_out.cc
+++ b/source/grid/grid_out.cc
@@ -910,7 +910,7 @@ GridOut::write_dx(const Triangulation<dim, spacedim> &tria,
 
       for (const auto &cell : tria.active_cell_iterators())
         {
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             {
               typename Triangulation<dim, spacedim>::face_iterator face =
                 cell->face(f);
@@ -940,7 +940,7 @@ GridOut::write_dx(const Triangulation<dim, spacedim> &tria,
       for (const auto &cell : tria.active_cell_iterators())
         {
           // Little trick to get -1 for the interior
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             {
               out << ' '
                   << static_cast<std::make_signed<types::boundary_id>::type>(
@@ -3357,7 +3357,7 @@ GridOut::write_mesh_per_processor_as_vtu(
           patch.data(3, vertex) = tria.locally_owned_subdomain();
         }
 
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (auto f : GeometryInfo<dim>::face_indices())
         patch.neighbors[f] = numbers::invalid_unsigned_int;
       patches.push_back(patch);
     }
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 8dca0dc61a..dbf5adf3f5 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -1250,7 +1250,7 @@ namespace GridTools
       endc = tria.end();
     for (; cell != endc; ++cell)
       {
-        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+        for (unsigned int i : GeometryInfo<dim>::face_indices())
           {
             const typename Triangulation<dim, spacedim>::face_iterator &face =
               cell->face(i);
@@ -2356,7 +2356,7 @@ namespace GridTools
     cell = triangulation.begin_active();
     for (; cell != endc; ++cell)
       {
-        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+        for (unsigned int i : GeometryInfo<dim>::face_indices())
           {
             if ((cell->at_boundary(i) == false) &&
                 (cell->neighbor(i)->is_active()))
@@ -2503,7 +2503,7 @@ namespace GridTools
         // received.
         if (cell->is_ghost())
           {
-            for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+            for (unsigned int i : GeometryInfo<dim>::face_indices())
               {
                 if (cell->at_boundary(i) == false)
                   {
@@ -2774,7 +2774,7 @@ namespace GridTools
       {
         const unsigned int index = cell->active_cell_index();
         cell_connectivity.add(index, index);
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if ((cell->at_boundary(f) == false) &&
               (cell->neighbor(f)->has_children() == false))
             {
@@ -3808,7 +3808,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 (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           {
             Assert(cell->face(f)->has_children(), ExcInternalError());
             Assert(cell->face(f)->refinement_case() ==
@@ -3922,7 +3922,7 @@ namespace GridTools
              tria.begin_active();
            cell != tria.end();
            ++cell)
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (auto f : GeometryInfo<dim>::face_indices())
           if (cell->face(f)->at_boundary())
             for (unsigned int e = 0; e < GeometryInfo<dim>::lines_per_face; ++e)
               {
@@ -3941,7 +3941,7 @@ namespace GridTools
            tria.begin_active();
          cell != tria.end();
          ++cell)
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (auto f : GeometryInfo<dim>::face_indices())
         if (cell->face(f)->at_boundary())
           {
             const auto         bid = cell->face(f)->boundary_id();
@@ -3988,7 +3988,7 @@ namespace GridTools
         cell->set_manifold_id(cell->material_id());
         if (compute_face_ids == true)
           {
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               {
                 if (cell->at_boundary(f) == false)
                   cell->face(f)->set_manifold_id(
@@ -4210,7 +4210,7 @@ namespace GridTools
         for (const auto &cell : tria.active_cell_iterators())
           {
             unsigned int boundary_face_counter = 0;
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               if (cell->face(f)->at_boundary())
                 boundary_face_counter++;
             if (boundary_face_counter > dim)
diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc
index 1b6334018a..47584b90e7 100644
--- a/source/grid/grid_tools_dof_handlers.cc
+++ b/source/grid/grid_tools_dof_handlers.cc
@@ -2235,7 +2235,7 @@ namespace GridTools
          cell != mesh.end(0);
          ++cell)
       {
-        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+        for (unsigned int i : GeometryInfo<dim>::face_indices())
           {
             const typename MeshType::face_iterator face = cell->face(i);
             if (face->at_boundary() && face->boundary_id() == b_id1)
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 94beff7582..2333a0b342 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -504,7 +504,7 @@ namespace
           cell = triangulation.begin(),
           endc = triangulation.end();
         for (; cell != endc; ++cell)
-          for (unsigned int q = 0; q < GeometryInfo<dim>::faces_per_cell; ++q)
+          for (unsigned int q : GeometryInfo<dim>::face_indices())
             ++quad_cell_count[cell->quad_index(q)];
         return quad_cell_count;
       }
@@ -854,7 +854,7 @@ namespace
                                                          endc =
                                                            triangulation.end();
     for (; cell != endc; ++cell)
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (auto f : GeometryInfo<dim>::face_indices())
         {
           const typename Triangulation<dim, spacedim>::face_iterator face =
             cell->face(f);
@@ -934,7 +934,7 @@ namespace
     // left_right_offset in this case as we want
     // the offset of the neighbor, not our own.
     for (cell = triangulation.begin(); cell != endc; ++cell)
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (auto f : GeometryInfo<dim>::face_indices())
         {
           const unsigned int offset =
             (cell->direction_flag() ?
@@ -1938,7 +1938,7 @@ namespace internal
                triangulation.begin_active();
              cell != triangulation.end();
              ++cell)
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             {
               (*triangulation.vertex_to_manifold_id_map_1d)
                 [cell->face(f)->vertex_index()] = numbers::flat_manifold_id;
@@ -3661,7 +3661,7 @@ namespace internal
               cell->child(c);
             for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
               --line_cell_count[child->line_index(l)];
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               --quad_cell_count[child->quad_index(f)];
           }
 
@@ -3752,7 +3752,7 @@ namespace internal
             cell->child(child)->clear_user_data();
             cell->child(child)->clear_user_flag();
 
-            for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+            for (auto f : GeometryInfo<dim>::face_indices())
               {
                 // set flags denoting deviations from
                 // standard orientation of faces back
@@ -10061,7 +10061,7 @@ namespace internal
           return true;
 
         const RefinementCase<dim> ref_case = cell->refinement_case();
-        for (unsigned int n = 0; n < GeometryInfo<dim>::faces_per_cell; ++n)
+        for (unsigned int n : GeometryInfo<dim>::face_indices())
           {
             // if the cell is not refined along that face, coarsening
             // will not change anything, so do nothing. the same
@@ -10345,7 +10345,7 @@ Triangulation<dim, spacedim>::set_all_manifold_ids_on_boundary(
     endc = this->end();
 
   for (; cell != endc; ++cell)
-    for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    for (auto f : GeometryInfo<dim>::face_indices())
       if (cell->face(f)->at_boundary())
         cell->face(f)->set_all_manifold_ids(m_number);
 }
@@ -10370,7 +10370,7 @@ Triangulation<dim, spacedim>::set_all_manifold_ids_on_boundary(
   for (; cell != endc; ++cell)
     {
       // loop on faces
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (auto f : GeometryInfo<dim>::face_indices())
         if (cell->face(f)->at_boundary() &&
             cell->face(f)->boundary_id() == b_id)
           {
@@ -13963,7 +13963,7 @@ namespace
     // count all neighbors that will be refined along the face of our
     // cell after the next step
     unsigned int count = 0;
-    for (unsigned int n = 0; n < GeometryInfo<dim>::faces_per_cell; ++n)
+    for (unsigned int n : GeometryInfo<dim>::face_indices())
       {
         const typename Triangulation<dim, spacedim>::cell_iterator neighbor =
           cell->neighbor(n);
diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc
index 78e8290a2c..be33ba35e2 100644
--- a/source/hp/dof_handler.cc
+++ b/source/hp/dof_handler.cc
@@ -1448,7 +1448,7 @@ namespace hp
     for (; cell != endc; ++cell)
       if (cell->is_locally_owned() && cell->at_boundary())
         {
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             if (cell->at_boundary(f))
               {
                 const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
@@ -1495,7 +1495,7 @@ namespace hp
     for (; cell != endc; ++cell)
       if (cell->is_locally_owned() && cell->at_boundary())
         {
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+          for (auto f : GeometryInfo<dim>::face_indices())
             if (cell->at_boundary(f) &&
                 (boundary_ids.find(cell->face(f)->boundary_id()) !=
                  boundary_ids.end()))
diff --git a/source/multigrid/mg_tools.cc b/source/multigrid/mg_tools.cc
index e9cd8a6e33..dbd96edf6c 100644
--- a/source/multigrid/mg_tools.cc
+++ b/source/multigrid/mg_tools.cc
@@ -785,7 +785,7 @@ namespace MGTools
         DoFTools::dof_couplings_from_component_couplings(fe, flux_mask);
 
     for (unsigned int i = 0; i < total_dofs; ++i)
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (auto f : GeometryInfo<dim>::face_indices())
         support_on_face(i, f) = fe.has_support_on_face(i, f);
 
     // Clear user flags because we will
@@ -974,7 +974,7 @@ namespace MGTools
       DoFTools::dof_couplings_from_component_couplings(fe, flux_mask);
 
     for (unsigned int i = 0; i < dofs_per_cell; ++i)
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (auto f : GeometryInfo<dim>::face_indices())
         support_on_face(i, f) = fe.has_support_on_face(i, f);
 
     for (; cell != endc; ++cell)
diff --git a/source/opencascade/utilities.cc b/source/opencascade/utilities.cc
index 42229eea32..4c5872724e 100644
--- a/source/opencascade/utilities.cc
+++ b/source/opencascade/utilities.cc
@@ -612,7 +612,7 @@ namespace OpenCASCADE
     unsigned int face_index;
 
     for (const auto &cell : triangulation.active_cell_iterators())
-      for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
+      for (unsigned int f : GeometryInfo<2>::face_indices())
         if (cell->face(f)->at_boundary())
           {
             // get global face and vertex indices
-- 
2.39.5