Use GeometryInfo::face_indices().

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)

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>());

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);
 

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()] =

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 &&

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 =

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])

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);

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())
                 {

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));
       }

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;

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;

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)

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)

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);

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)

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()

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;