From e2d7a2ff171a111480307d97dce2aa4db063a953 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <arndtd@ornl.gov>
Date: Mon, 4 Nov 2019 09:29:00 -0500
Subject: [PATCH] Use face_iterators in tutorials

---
 examples/step-10/step-10.cc |  8 ++---
 examples/step-14/step-14.cc | 16 ++++------
 examples/step-18/step-18.cc | 16 +++++-----
 examples/step-20/step-20.cc |  8 ++---
 examples/step-21/step-21.cc |  8 ++---
 examples/step-22/step-22.cc |  6 ++--
 examples/step-24/step-24.cc |  6 ++--
 examples/step-29/step-29.cc | 11 ++++---
 examples/step-42/step-42.cc | 59 ++++++++++++++-----------------------
 examples/step-43/step-43.cc | 18 +++++------
 examples/step-44/step-44.cc | 23 +++++++--------
 examples/step-46/step-46.cc |  7 ++---
 examples/step-49/step-49.cc |  7 ++---
 examples/step-51/step-51.cc | 11 ++++---
 examples/step-52/step-52.cc | 11 ++++---
 examples/step-53/step-53.cc |  4 +--
 examples/step-54/step-54.cc |  4 +--
 examples/step-61/step-61.cc | 34 +++++++++------------
 examples/step-65/step-65.cc |  7 ++---
 examples/step-7/step-7.cc   | 17 ++++-------
 examples/step-9/step-9.cc   |  7 ++---
 21 files changed, 118 insertions(+), 170 deletions(-)

diff --git a/examples/step-10/step-10.cc b/examples/step-10/step-10.cc
index 8b6215f52b..36662569b3 100644
--- a/examples/step-10/step-10.cc
+++ b/examples/step-10/step-10.cc
@@ -355,14 +355,12 @@ namespace Step10
             // added to the long double variable `perimeter'.
             long double perimeter = 0;
             for (const auto &cell : dof_handler.active_cell_iterators())
-              for (unsigned int face_no = 0;
-                   face_no < GeometryInfo<dim>::faces_per_cell;
-                   ++face_no)
-                if (cell->face(face_no)->at_boundary())
+              for (const auto &face : cell->face_iterators())
+                if (face->at_boundary())
                   {
                     // We reinit the FEFaceValues object with the cell
                     // iterator and the number of the face.
-                    fe_face_values.reinit(cell, face_no);
+                    fe_face_values.reinit(cell, face);
                     for (unsigned int i = 0;
                          i < fe_face_values.n_quadrature_points;
                          ++i)
diff --git a/examples/step-14/step-14.cc b/examples/step-14/step-14.cc
index 4f73f70132..56068529fd 100644
--- a/examples/step-14/step-14.cc
+++ b/examples/step-14/step-14.cc
@@ -2187,10 +2187,8 @@ namespace Step14
       FaceIntegrals face_integrals;
       for (const auto &cell :
            DualSolver<dim>::dof_handler.active_cell_iterators())
-        for (unsigned int face_no = 0;
-             face_no < GeometryInfo<dim>::faces_per_cell;
-             ++face_no)
-          face_integrals[cell->face(face_no)] = -1e20;
+        for (const auto &face : cell->face_iterators())
+          face_integrals[face] = -1e20;
 
       auto worker = [this,
                      &error_indicators,
@@ -2229,15 +2227,11 @@ namespace Step14
       for (const auto &cell :
            DualSolver<dim>::dof_handler.active_cell_iterators())
         {
-          for (unsigned int face_no = 0;
-               face_no < GeometryInfo<dim>::faces_per_cell;
-               ++face_no)
+          for (const auto &face : cell->face_iterators())
             {
-              Assert(face_integrals.find(cell->face(face_no)) !=
-                       face_integrals.end(),
+              Assert(face_integrals.find(face) != face_integrals.end(),
                      ExcInternalError());
-              error_indicators(present_cell) -=
-                0.5 * face_integrals[cell->face(face_no)];
+              error_indicators(present_cell) -= 0.5 * face_integrals[face];
             }
           ++present_cell;
         }
diff --git a/examples/step-18/step-18.cc b/examples/step-18/step-18.cc
index f8a8129214..436d000685 100644
--- a/examples/step-18/step-18.cc
+++ b/examples/step-18/step-18.cc
@@ -774,22 +774,22 @@ namespace Step18
   {
     const double inner_radius = 0.8, outer_radius = 1;
     GridGenerator::cylinder_shell(triangulation, 3, inner_radius, outer_radius);
-    for (auto &cell : triangulation.active_cell_iterators())
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->at_boundary())
+    for (const auto &cell : triangulation.active_cell_iterators())
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary())
           {
-            const Point<dim> face_center = cell->face(f)->center();
+            const Point<dim> face_center = face->center();
 
             if (face_center[2] == 0)
-              cell->face(f)->set_boundary_id(0);
+              face->set_boundary_id(0);
             else if (face_center[2] == 3)
-              cell->face(f)->set_boundary_id(1);
+              face->set_boundary_id(1);
             else if (std::sqrt(face_center[0] * face_center[0] +
                                face_center[1] * face_center[1]) <
                      (inner_radius + outer_radius) / 2)
-              cell->face(f)->set_boundary_id(2);
+              face->set_boundary_id(2);
             else
-              cell->face(f)->set_boundary_id(3);
+              face->set_boundary_id(3);
           }
 
     // Once all this is done, we can refine the mesh once globally:
diff --git a/examples/step-20/step-20.cc b/examples/step-20/step-20.cc
index c898d57861..c1713c06e7 100644
--- a/examples/step-20/step-20.cc
+++ b/examples/step-20/step-20.cc
@@ -523,12 +523,10 @@ namespace Step20
               local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q);
             }
 
-        for (unsigned int face_n = 0;
-             face_n < GeometryInfo<dim>::faces_per_cell;
-             ++face_n)
-          if (cell->at_boundary(face_n))
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary())
             {
-              fe_face_values.reinit(cell, face_n);
+              fe_face_values.reinit(cell, face);
 
               pressure_boundary_values.value_list(
                 fe_face_values.get_quadrature_points(), boundary_values);
diff --git a/examples/step-21/step-21.cc b/examples/step-21/step-21.cc
index a739a74432..aaef851b44 100644
--- a/examples/step-21/step-21.cc
+++ b/examples/step-21/step-21.cc
@@ -714,12 +714,10 @@ namespace Step21
 
         // Next, we also have to deal with the pressure boundary values. This,
         // again is as in step-20:
-        for (unsigned int face_no = 0;
-             face_no < GeometryInfo<dim>::faces_per_cell;
-             ++face_no)
-          if (cell->at_boundary(face_no))
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary())
             {
-              fe_face_values.reinit(cell, face_no);
+              fe_face_values.reinit(cell, face);
 
               pressure_boundary_values.value_list(
                 fe_face_values.get_quadrature_points(), boundary_values);
diff --git a/examples/step-22/step-22.cc b/examples/step-22/step-22.cc
index dcaff91ce8..e49db80492 100644
--- a/examples/step-22/step-22.cc
+++ b/examples/step-22/step-22.cc
@@ -1023,9 +1023,9 @@ namespace Step22
     // the last coordinate direction. See the example description above for
     // details.
     for (const auto &cell : triangulation.active_cell_iterators())
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->center()[dim - 1] == 0)
-          cell->face(f)->set_all_boundary_ids(1);
+      for (const auto &face : cell->face_iterators())
+        if (face->center()[dim - 1] == 0)
+          face->set_all_boundary_ids(1);
 
 
     // We then apply an initial refinement before solving for the first
diff --git a/examples/step-24/step-24.cc b/examples/step-24/step-24.cc
index 17fa4fb5ef..611e41c205 100644
--- a/examples/step-24/step-24.cc
+++ b/examples/step-24/step-24.cc
@@ -341,12 +341,12 @@ namespace Step24
       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
       for (const auto &cell : dof_handler.active_cell_iterators())
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
-          if (cell->at_boundary(f))
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary())
             {
               cell_matrix = 0;
 
-              fe_values.reinit(cell, f);
+              fe_values.reinit(cell, face);
 
               for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
diff --git a/examples/step-29/step-29.cc b/examples/step-29/step-29.cc
index 2cbc7efc06..3146ed1307 100644
--- a/examples/step-29/step-29.cc
+++ b/examples/step-29/step-29.cc
@@ -454,13 +454,12 @@ namespace Step29
     GridGenerator::subdivided_hyper_cube(triangulation, 5, 0, 1);
 
     for (auto &cell : triangulation.cell_iterators())
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
-        if (cell->face(face)->at_boundary() &&
-            ((cell->face(face)->center() - transducer).norm_square() < 0.01))
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary() &&
+            ((face->center() - transducer).norm_square() < 0.01))
           {
-            cell->face(face)->set_boundary_id(1);
-            cell->face(face)->set_manifold_id(1);
+            face->set_boundary_id(1);
+            face->set_manifold_id(1);
           }
     // For the circle part of the transducer lens, a SphericalManifold object
     // is used (which, of course, in 2D just represents a circle), with center
diff --git a/examples/step-42/step-42.cc b/examples/step-42/step-42.cc
index 799c59e8be..3b83bfeda2 100644
--- a/examples/step-42/step-42.cc
+++ b/examples/step-42/step-42.cc
@@ -934,23 +934,18 @@ namespace Step42
         GridGenerator::hyper_rectangle(triangulation, p1, p2);
 
         for (const auto &cell : triangulation.active_cell_iterators())
-          for (unsigned int face_no = 0;
-               face_no < GeometryInfo<dim>::faces_per_cell;
-               ++face_no)
-            if (cell->face(face_no)->at_boundary())
+          for (const auto &face : cell->face_iterators())
+            if (face->at_boundary())
               {
-                if (std::fabs(cell->face(face_no)->center()[2] - p2[2]) < 1e-12)
-                  cell->face(face_no)->set_boundary_id(1);
-                if (std::fabs(cell->face(face_no)->center()[0] - p1[0]) <
-                      1e-12 ||
-                    std::fabs(cell->face(face_no)->center()[0] - p2[0]) <
-                      1e-12 ||
-                    std::fabs(cell->face(face_no)->center()[1] - p1[1]) <
-                      1e-12 ||
-                    std::fabs(cell->face(face_no)->center()[1] - p2[1]) < 1e-12)
-                  cell->face(face_no)->set_boundary_id(8);
-                if (std::fabs(cell->face(face_no)->center()[2] - p1[2]) < 1e-12)
-                  cell->face(face_no)->set_boundary_id(6);
+                if (std::fabs(face->center()[2] - p2[2]) < 1e-12)
+                  face->set_boundary_id(1);
+                if (std::fabs(face->center()[0] - p1[0]) < 1e-12 ||
+                    std::fabs(face->center()[0] - p2[0]) < 1e-12 ||
+                    std::fabs(face->center()[1] - p1[1]) < 1e-12 ||
+                    std::fabs(face->center()[1] - p2[1]) < 1e-12)
+                  face->set_boundary_id(8);
+                if (std::fabs(face->center()[2] - p1[2]) < 1e-12)
+                  face->set_boundary_id(6);
               }
       }
 
@@ -1145,10 +1140,8 @@ namespace Step42
 
     for (const auto &cell : dof_handler.active_cell_iterators())
       if (cell->is_locally_owned())
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
-          if (cell->face(face)->at_boundary() &&
-              cell->face(face)->boundary_id() == 1)
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary() && face->boundary_id() == 1)
             {
               fe_values_face.reinit(cell, face);
               cell_matrix = 0;
@@ -1236,13 +1229,11 @@ namespace Step42
 
     for (const auto &cell : dof_handler.active_cell_iterators())
       if (!cell->is_artificial())
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
-          if (cell->face(face)->at_boundary() &&
-              cell->face(face)->boundary_id() == 1)
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary() && face->boundary_id() == 1)
             {
               fe_values_face.reinit(cell, face);
-              cell->face(face)->get_dof_indices(dof_indices);
+              face->get_dof_indices(dof_indices);
 
               for (unsigned int q_point = 0; q_point < n_face_q_points;
                    ++q_point)
@@ -1425,10 +1416,8 @@ namespace Step42
                 }
             }
 
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
-            if (cell->face(face)->at_boundary() &&
-                cell->face(face)->boundary_id() == 1)
+          for (const auto &face : cell->face_iterators())
+            if (face->at_boundary() && face->boundary_id() == 1)
               {
                 fe_values_face.reinit(cell, face);
 
@@ -1557,10 +1546,8 @@ namespace Step42
                 }
             }
 
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
-            if (cell->face(face)->at_boundary() &&
-                cell->face(face)->boundary_id() == 1)
+          for (const auto &face : cell->face_iterators())
+            if (face->at_boundary() && face->boundary_id() == 1)
               {
                 fe_values_face.reinit(cell, face);
 
@@ -2125,10 +2112,8 @@ namespace Step42
 
     for (const auto &cell : dof_handler.active_cell_iterators())
       if (cell->is_locally_owned())
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
-          if (cell->face(face)->at_boundary() &&
-              cell->face(face)->boundary_id() == 1)
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary() && face->boundary_id() == 1)
             {
               fe_values_face.reinit(cell, face);
 
diff --git a/examples/step-43/step-43.cc b/examples/step-43/step-43.cc
index a6feeedbca..f40725c55a 100644
--- a/examples/step-43/step-43.cc
+++ b/examples/step-43/step-43.cc
@@ -1125,12 +1125,10 @@ namespace Step43
               }
           }
 
-        for (unsigned int face_no = 0;
-             face_no < GeometryInfo<dim>::faces_per_cell;
-             ++face_no)
-          if (cell->at_boundary(face_no))
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary())
             {
-              darcy_fe_face_values.reinit(cell, face_no);
+              darcy_fe_face_values.reinit(cell, face);
 
               pressure_boundary_values.value_list(
                 darcy_fe_face_values.get_quadrature_points(), boundary_values);
@@ -1310,13 +1308,11 @@ namespace Step43
                                           global_S_variation,
                                           local_dof_indices);
 
-        for (unsigned int face_no = 0;
-             face_no < GeometryInfo<dim>::faces_per_cell;
-             ++face_no)
-          if (cell->at_boundary(face_no))
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary())
             {
-              darcy_fe_face_values.reinit(darcy_cell, face_no);
-              saturation_fe_face_values.reinit(cell, face_no);
+              darcy_fe_face_values.reinit(darcy_cell, face);
+              saturation_fe_face_values.reinit(cell, face);
               assemble_saturation_rhs_boundary_term(saturation_fe_face_values,
                                                     darcy_fe_face_values,
                                                     local_dof_indices);
diff --git a/examples/step-44/step-44.cc b/examples/step-44/step-44.cc
index aebee6b5c3..0030c41abb 100644
--- a/examples/step-44/step-44.cc
+++ b/examples/step-44/step-44.cc
@@ -1467,22 +1467,21 @@ namespace Step44
     // +y surface and will get boundary ID 6 (zero through five are already
     // used when creating the six faces of the cube domain):
     for (const auto &cell : triangulation.active_cell_iterators())
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (const auto &face : cell->face_iterators())
         {
-          if (cell->face(face)->at_boundary() == true &&
-              cell->face(face)->center()[1] == 1.0 * parameters.scale)
+          if (face->at_boundary() == true &&
+              face->center()[1] == 1.0 * parameters.scale)
             {
               if (dim == 3)
                 {
-                  if (cell->face(face)->center()[0] < 0.5 * parameters.scale &&
-                      cell->face(face)->center()[2] < 0.5 * parameters.scale)
-                    cell->face(face)->set_boundary_id(6);
+                  if (face->center()[0] < 0.5 * parameters.scale &&
+                      face->center()[2] < 0.5 * parameters.scale)
+                    face->set_boundary_id(6);
                 }
               else
                 {
-                  if (cell->face(face)->center()[0] < 0.5 * parameters.scale)
-                    cell->face(face)->set_boundary_id(6);
+                  if (face->center()[0] < 0.5 * parameters.scale)
+                    face->set_boundary_id(6);
                 }
             }
         }
@@ -2333,10 +2332,8 @@ namespace Step44
     // Next we assemble the Neumann contribution. We first check to see it the
     // cell face exists on a boundary on which a traction is applied and add
     // the contribution if this is the case.
-    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-         ++face)
-      if (cell->face(face)->at_boundary() == true &&
-          cell->face(face)->boundary_id() == 6)
+    for (const auto &face : cell->face_iterators())
+      if (face->at_boundary() == true && face->boundary_id() == 6)
         {
           scratch.fe_face_values.reinit(cell, face);
 
diff --git a/examples/step-46/step-46.cc b/examples/step-46/step-46.cc
index 58e5d0bf92..262fbcf0e5 100644
--- a/examples/step-46/step-46.cc
+++ b/examples/step-46/step-46.cc
@@ -281,10 +281,9 @@ namespace Step46
     GridGenerator::subdivided_hyper_cube(triangulation, 8, -1, 1);
 
     for (const auto &cell : triangulation.active_cell_iterators())
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
-        if (cell->face(f)->at_boundary() &&
-            (cell->face(f)->center()[dim - 1] == 1))
-          cell->face(f)->set_all_boundary_ids(1);
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary() && (face->center()[dim - 1] == 1))
+          face->set_all_boundary_ids(1);
 
 
     for (const auto &cell : dof_handler.active_cell_iterators())
diff --git a/examples/step-49/step-49.cc b/examples/step-49/step-49.cc
index 60eed23c7f..056311f212 100644
--- a/examples/step-49/step-49.cc
+++ b/examples/step-49/step-49.cc
@@ -74,11 +74,10 @@ void print_mesh_info(const Triangulation<dim> &triangulation,
     std::map<types::boundary_id, unsigned int> boundary_count;
     for (const auto &cell : triangulation.active_cell_iterators())
       {
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
+        for (const auto &face : cell->face_iterators())
           {
-            if (cell->face(face)->at_boundary())
-              boundary_count[cell->face(face)->boundary_id()]++;
+            if (face->at_boundary())
+              boundary_count[face->boundary_id()]++;
           }
       }
 
diff --git a/examples/step-51/step-51.cc b/examples/step-51/step-51.cc
index 433632338d..946063f9da 100644
--- a/examples/step-51/step-51.cc
+++ b/examples/step-51/step-51.cc
@@ -1333,12 +1333,11 @@ namespace Step51
     // refinement, the flags are set in every refinement step, not just at the
     // beginning.
     for (const auto &cell : triangulation.cell_iterators())
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
-        if (cell->face(face)->at_boundary())
-          if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12) ||
-              (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12))
-            cell->face(face)->set_boundary_id(1);
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary())
+          if ((std::fabs(face->center()(0) - (-1)) < 1e-12) ||
+              (std::fabs(face->center()(1) - (-1)) < 1e-12))
+            face->set_boundary_id(1);
   }
 
   // @sect4{HDG::run}
diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc
index 1dac0b7110..4f7548bb28 100644
--- a/examples/step-52/step-52.cc
+++ b/examples/step-52/step-52.cc
@@ -642,14 +642,13 @@ namespace Step52
     triangulation.refine_global(4);
 
     for (const auto &cell : triangulation.active_cell_iterators())
-      for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
-        if (cell->face(f)->at_boundary())
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary())
           {
-            if ((cell->face(f)->center()[0] == 0.) ||
-                (cell->face(f)->center()[0] == 5.))
-              cell->face(f)->set_boundary_id(1);
+            if ((face->center()[0] == 0.) || (face->center()[0] == 5.))
+              face->set_boundary_id(1);
             else
-              cell->face(f)->set_boundary_id(0);
+              face->set_boundary_id(0);
           }
 
     // Next, we set up the linear systems and fill them with content so that
diff --git a/examples/step-53/step-53.cc b/examples/step-53/step-53.cc
index 92f7b81286..03cc1fdc1c 100644
--- a/examples/step-53/step-53.cc
+++ b/examples/step-53/step-53.cc
@@ -385,8 +385,8 @@ namespace Step53
     for (unsigned int i = 0; i < 6; ++i)
       {
         for (const auto &cell : triangulation.active_cell_iterators())
-          for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f)
-            if (cell->face(f)->boundary_id() == 5)
+          for (const auto &face : cell->face_iterators())
+            if (face->boundary_id() == 5)
               {
                 cell->set_refine_flag();
                 break;
diff --git a/examples/step-54/step-54.cc b/examples/step-54/step-54.cc
index 9e4ee625d2..36074c8d79 100644
--- a/examples/step-54/step-54.cc
+++ b/examples/step-54/step-54.cc
@@ -231,8 +231,8 @@ namespace Step54
     Triangulation<2, 3>::active_cell_iterator cell = tria.begin_active();
     cell->set_manifold_id(1);
 
-    for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
-      cell->face(f)->set_manifold_id(2);
+    for (const auto &face : cell->face_iterators())
+      face->set_manifold_id(2);
 
     // Once both the CAD geometry and the initial mesh have been
     // imported and digested, we use the CAD surfaces and curves to
diff --git a/examples/step-61/step-61.cc b/examples/step-61/step-61.cc
index 32bf8dcdbf..ce3efe9d4e 100644
--- a/examples/step-61/step-61.cc
+++ b/examples/step-61/step-61.cc
@@ -441,7 +441,7 @@ namespace Step61
     // @p face component of the shape functions.
     const FEValuesExtractors::Vector velocities(0);
     const FEValuesExtractors::Scalar interior(0);
-    const FEValuesExtractors::Scalar face(1);
+    const FEValuesExtractors::Scalar face_extractor(1);
 
     // This finally gets us in position to loop over all cells. On
     // each cell, we will first calculate the various cell matrices
@@ -514,12 +514,10 @@ namespace Step61
         // of the polynomial space and the dot product of a basis function of
         // the Raviart-Thomas space and the normal vector. So we loop over all
         // the faces of the element and obtain the normal vector.
-        for (unsigned int face_n = 0;
-             face_n < GeometryInfo<dim>::faces_per_cell;
-             ++face_n)
+        for (const auto &face : cell->face_iterators())
           {
-            fe_face_values.reinit(cell, face_n);
-            fe_face_values_rt.reinit(cell_rt, face_n);
+            fe_face_values.reinit(cell, face);
+            fe_face_values_rt.reinit(cell_rt, face);
 
             for (unsigned int q = 0; q < n_face_q_points; ++q)
               {
@@ -532,7 +530,7 @@ namespace Step61
                     for (unsigned int j = 0; j < dofs_per_cell; ++j)
                       {
                         const double phi_j_face =
-                          fe_face_values[face].value(j, q);
+                          fe_face_values[face_extractor].value(j, q);
 
                         cell_matrix_G(i, j) +=
                           ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
@@ -747,7 +745,7 @@ namespace Step61
     const FEValuesExtractors::Vector velocities(0);
     const FEValuesExtractors::Scalar pressure(dim);
     const FEValuesExtractors::Scalar interior(0);
-    const FEValuesExtractors::Scalar face(1);
+    const FEValuesExtractors::Scalar face_extractor(1);
     const FEValuesExtractors::Vector velocities_dgrt(0);
 
     const ExactVelocity<dim> exact_velocity;
@@ -827,12 +825,10 @@ namespace Step61
                 }
             }
 
-        for (unsigned int face_n = 0;
-             face_n < GeometryInfo<dim>::faces_per_cell;
-             ++face_n)
+        for (const auto &face : cell->face_iterators())
           {
-            fe_face_values.reinit(cell, face_n);
-            fe_face_values_rt.reinit(cell_rt, face_n);
+            fe_face_values.reinit(cell, face);
+            fe_face_values_rt.reinit(cell_rt, face);
 
             for (unsigned int q = 0; q < n_face_q_points; ++q)
               {
@@ -845,7 +841,7 @@ namespace Step61
                     for (unsigned int j = 0; j < dofs_per_cell; ++j)
                       {
                         const double phi_j_face =
-                          fe_face_values[face].value(j, q);
+                          fe_face_values[face_extractor].value(j, q);
 
                         cell_matrix_G(i, j) +=
                           ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
@@ -909,13 +905,11 @@ namespace Step61
         // the $L_2$ flux error on the cell and add it to the global
         // error.
         const double cell_area = cell->measure();
-        for (unsigned int face_n = 0;
-             face_n < GeometryInfo<dim>::faces_per_cell;
-             ++face_n)
+        for (const auto &face : cell->face_iterators())
           {
-            const double face_length = cell->face(face_n)->measure();
-            fe_face_values.reinit(cell, face_n);
-            fe_face_values_rt.reinit(cell_rt, face_n);
+            const double face_length = face->measure();
+            fe_face_values.reinit(cell, face);
+            fe_face_values_rt.reinit(cell_rt, face);
 
             double L2_err_flux_face_sqr_local = 0;
             for (unsigned int q = 0; q < n_face_q_points_rt; ++q)
diff --git a/examples/step-65/step-65.cc b/examples/step-65/step-65.cc
index d9672c21a3..2bba62ff20 100644
--- a/examples/step-65/step-65.cc
+++ b/examples/step-65/step-65.cc
@@ -227,19 +227,18 @@ namespace Step65
 
     for (const auto &cell : triangulation.cell_iterators())
       {
-        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        for (const auto &face : cell->face_iterators())
           {
             bool face_at_sphere_boundary = true;
             for (unsigned int v = 0;
                  v < GeometryInfo<dim - 1>::vertices_per_cell;
                  ++v)
               {
-                if (std::abs(cell->face(f)->vertex(v).norm_square() - 0.25) >
-                    1e-12)
+                if (std::abs(face->vertex(v).norm_square() - 0.25) > 1e-12)
                   face_at_sphere_boundary = false;
               }
             if (face_at_sphere_boundary)
-              cell->face(f)->set_all_manifold_ids(1);
+              face->set_all_manifold_ids(1);
           }
         if (cell->center().norm_square() < 0.25)
           cell->set_material_id(1);
diff --git a/examples/step-7/step-7.cc b/examples/step-7/step-7.cc
index 10d8833015..b0fbccf438 100644
--- a/examples/step-7/step-7.cc
+++ b/examples/step-7/step-7.cc
@@ -641,11 +641,8 @@ namespace Step7
         // <code>run()</code> function further below. (The default value of
         // boundary indicators is <code>0</code>, so faces can only have an
         // indicator equal to <code>1</code> if we have explicitly set it.)
-        for (unsigned int face_number = 0;
-             face_number < GeometryInfo<dim>::faces_per_cell;
-             ++face_number)
-          if (cell->face(face_number)->at_boundary() &&
-              (cell->face(face_number)->boundary_id() == 1))
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary() && (face->boundary_id() == 1))
             {
               // If we came into here, then we have found an external face
               // belonging to Gamma2. Next, we have to compute the values of
@@ -653,7 +650,7 @@ namespace Step7
               // need for the computation of the contour integral. This is
               // done using the <code>reinit</code> function which we already
               // know from the FEValue class:
-              fe_face_values.reinit(cell, face_number);
+              fe_face_values.reinit(cell, face);
 
               // And we can then perform the integration by using a loop over
               // all quadrature points.
@@ -971,14 +968,12 @@ namespace Step7
             triangulation.refine_global(3);
 
             for (const auto &cell : triangulation.cell_iterators())
-              for (unsigned int face_number = 0;
-                   face_number < GeometryInfo<dim>::faces_per_cell;
-                   ++face_number)
+              for (const auto &face : cell->face_iterators())
                 {
-                  const auto center = cell->face(face_number)->center();
+                  const auto center = face->center();
                   if ((std::fabs(center(0) - (-1)) < 1e-12) ||
                       (std::fabs(center(1) - (-1)) < 1e-12))
-                    cell->face(face_number)->set_boundary_id(1);
+                    face->set_boundary_id(1);
                 }
           }
         else
diff --git a/examples/step-9/step-9.cc b/examples/step-9/step-9.cc
index 836afd2be5..8832e46f0f 100644
--- a/examples/step-9/step-9.cc
+++ b/examples/step-9/step-9.cc
@@ -675,15 +675,14 @@ namespace Step9
     // the direction of flow at this point; we obtain this information
     // using the FEFaceValues object and only decide within the main loop
     // whether a quadrature point is on the inflow boundary.
-    for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell;
-         ++face_n)
-      if (cell->face(face_n)->at_boundary())
+    for (const auto &face : cell->face_iterators())
+      if (face->at_boundary())
         {
           // Ok, this face of the present cell is on the boundary of the
           // domain. Just as for the usual FEValues object which we have
           // used in previous examples and also above, we have to
           // reinitialize the FEFaceValues object for the present face:
-          scratch_data.fe_face_values.reinit(cell, face_n);
+          scratch_data.fe_face_values.reinit(cell, face);
 
           // For the quadrature points at hand, we ask for the values of
           // the inflow function and for the direction of flow:
-- 
2.39.5