]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use cell->face_iterators() or GeometryInfo::face_indices() in the tutorials. 9467/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 3 Feb 2020 01:47:22 +0000 (18:47 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 3 Feb 2020 01:47:22 +0000 (18:47 -0700)
13 files changed:
examples/step-14/step-14.cc
examples/step-20/doc/intro.dox
examples/step-21/step-21.cc
examples/step-29/step-29.cc
examples/step-30/step-30.cc
examples/step-33/step-33.cc
examples/step-46/doc/intro.dox
examples/step-46/step-46.cc
examples/step-49/doc/results.dox
examples/step-51/step-51.cc
examples/step-60/step-60.cc
examples/step-61/step-61.cc
examples/step-9/step-9.cc

index 56068529fd95828741f163ddb15d117aba0fe5a1..2813ba789398de05a7b13b32f9c87407b788060a 100644 (file)
@@ -2273,9 +2273,7 @@ namespace Step14
       // After computing the cell terms, turn to the face terms. For this,
       // loop over all faces of the present cell, and see whether
       // something needs to be computed on it:
-      for (unsigned int face_no = 0;
-           face_no < GeometryInfo<dim>::faces_per_cell;
-           ++face_no)
+      for (unsigned int face_no : GeometryInfo<dim>::face_indices())
         {
           // First, if this face is part of the boundary, then there is
           // nothing to do. However, to make things easier when summing up
index 22b5ea333e57b4a7f31a7868c07c7054d87dfe85..f89ce1546e99494a67dec75a9326b3082c2f8689 100644 (file)
@@ -331,24 +331,20 @@ to loop over all boundary faces and integrate there. The mechanism works in
 the same way as above, i.e. the extractor classes also work on FEFaceValues objects:
 
 @code
-      for (unsigned int face_no=0;
-           face_no<GeometryInfo<dim>::faces_per_cell;
-           ++face_no)
-        if (cell->at_boundary(face_no))
-          {
-            fe_face_values.reinit (cell, face_no);
-
-            pressure_boundary_values
-              .value_list (fe_face_values.get_quadrature_points(),
-                           boundary_values);
-
-            for (unsigned int q=0; q<n_face_q_points; ++q)
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                local_rhs(i) += -(fe_face_values[velocities].value (i, q) *
-                                  fe_face_values.normal_vector(q) *
-                                  boundary_values[q] *
-                                  fe_face_values.JxW(q));
-          }
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary())
+            {
+              fe_face_values.reinit(cell, face);
+
+              pressure_boundary_values.value_list(
+                fe_face_values.get_quadrature_points(), boundary_values);
+
+              for (unsigned int q = 0; q < n_face_q_points; ++q)
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  local_rhs(i) += -(fe_face_values[velocities].value(i, q) *
+                                    fe_face_values.normal_vector(q) *
+                                    boundary_values[q] *
+                                    fe_face_values.JxW(q));
 @endcode
 
 You will find the exact same code as above in the sources for the present
index aaef851b44c5900b7bc3ec0a55d061932036c2ee..4104101f0ecb3fcd6fe805578630faf7bf291cd6 100644 (file)
@@ -842,9 +842,7 @@ namespace Step21
         //
         // All this is a bit tricky, but has been explained in some detail
         // already in step-9. Take a look there how this is supposed to work!
-        for (unsigned int face_no = 0;
-             face_no < GeometryInfo<dim>::faces_per_cell;
-             ++face_no)
+        for (unsigned int face_no : GeometryInfo<dim>::face_indices())
           {
             fe_face_values.reinit(cell, face_no);
 
index 3146ed1307702d12d3d9bad9e942d2c544cadf17..0b196c4306d4b626dea87e243422e660f0d0fd4c 100644 (file)
@@ -651,27 +651,26 @@ namespace Step29
         // is at the boundary, and second has the correct boundary indicator
         // associated with $\Gamma_2$, the part of the boundary where we have
         // absorbing boundary conditions:
-        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-             ++face)
-          if (cell->face(face)->at_boundary() &&
-              (cell->face(face)->boundary_id() == 0))
+        for (unsigned int face_no : GeometryInfo<dim>::face_indices())
+          if (cell->face(face_no)->at_boundary() &&
+              (cell->face(face_no)->boundary_id() == 0))
             {
               // These faces will certainly contribute to the off-diagonal
               // blocks of the system matrix, so we ask the FEFaceValues
               // object to provide us with the shape function values on this
               // face:
-              fe_face_values.reinit(cell, face);
+              fe_face_values.reinit(cell, face_no);
 
 
               // Next, we loop through all DoFs of the current cell to find
               // pairs that belong to different components and both have
-              // support on the current face:
+              // support on the current face_no:
               for (unsigned int i = 0; i < dofs_per_cell; ++i)
                 for (unsigned int j = 0; j < dofs_per_cell; ++j)
                   if ((fe.system_to_component_index(i).first !=
                        fe.system_to_component_index(j).first) &&
-                      fe.has_support_on_face(i, face) &&
-                      fe.has_support_on_face(j, face))
+                      fe.has_support_on_face(i, face_no) &&
+                      fe.has_support_on_face(j, face_no))
                     // The check whether shape functions have support on a
                     // face is not strictly necessary: if we don't check for
                     // it we would simply add up terms to the local cell
index d2ae11527118073346748e4bef81866d5297a90d..6d78192309f7c0aeb192d2f0c56690caa8d56046 100644 (file)
@@ -451,9 +451,7 @@ namespace Step30
 
         cell->get_dof_indices(dofs);
 
-        for (unsigned int face_no = 0;
-             face_no < GeometryInfo<dim>::faces_per_cell;
-             ++face_no)
+        for (unsigned int face_no : GeometryInfo<dim>::face_indices())
           {
             const auto face = cell->face(face_no);
 
@@ -726,9 +724,7 @@ namespace Step30
           Point<dim> jump;
           Point<dim> area;
 
-          for (unsigned int face_no = 0;
-               face_no < GeometryInfo<dim>::faces_per_cell;
-               ++face_no)
+          for (unsigned int face_no : GeometryInfo<dim>::face_indices())
             {
               const auto face = cell->face(face_no);
 
index 35cb2cd473937b42980fdf64f8c19bc30c78a903..1f5d3bf6c653e3fbdbad96e314652cea8c53eb14 100644 (file)
@@ -1486,9 +1486,7 @@ namespace Step33
         // whether we are working on an external or internal face; if it is an
         // external face, the fourth argument denoting the degrees of freedom
         // indices of the neighbor is ignored, so we pass an empty vector):
-        for (unsigned int face_no = 0;
-             face_no < GeometryInfo<dim>::faces_per_cell;
-             ++face_no)
+        for (unsigned int face_no : GeometryInfo<dim>::face_indices())
           if (cell->at_boundary(face_no))
             {
               fe_v_face.reinit(cell, face_no);
index a0fa0fdac7e63b9c4119e4362bc5156b35bd5543..b77ed76f19cf81336e9c4f2ef88c08f62becedc9 100644 (file)
@@ -465,44 +465,38 @@ ensure that the velocity degrees of freedom on this face are
 zero. Some care is necessary to deal with the case that the adjacent
 solid cell is refined, yielding the following code:
 @code
-    std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
-    for (typename hp::DoFHandler<dim>::active_cell_iterator
-          cell = dof_handler.begin_active();
-        cell != dof_handler.end(); ++cell)
-      if (cell_is_in_fluid_domain (cell))
-       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-         if (!cell->at_boundary(f))
-           {
-             bool face_is_on_interface = false;
-
-             if ((cell->neighbor(f)->has_children() == false)
-                 &&
-                 (cell_is_in_solid_domain (cell->neighbor(f))))
-               face_is_on_interface = true;
-             else if (cell->neighbor(f)->has_children() == true)
-               {
-                                                  // The neighbor does
-                                                  // have
-                                                  // children. See if
-                                                  // any of the cells
-                                                  // on the other
-                                                  // side are elastic
-                 for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
-                   if (cell_is_in_solid_domain (cell->neighbor_child_on_subface(f, sf)))
-                     {
-                       face_is_on_interface = true;
-                       break;
-                     }
-               }
-
-             if (face_is_on_interface)
-               {
-                 cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
-                 for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
-                   if (stokes_fe.face_system_to_component_index(i).first < dim)
-                     constraints.add_line (local_face_dof_indices[i]);
-               }
-           }
+std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
+for (const auto &cell: dof_handler.active_cell_iterators())
+  if (cell_is_in_fluid_domain (cell))
+    for (unsigned int f : GeometryInfo<dim>::face_indices())
+      if (!cell->at_boundary(f))
+        {
+          bool face_is_on_interface = false;
+
+          if ((cell->neighbor(f)->has_children() == false)
+                 &&
+                 (cell_is_in_solid_domain (cell->neighbor(f))))
+               face_is_on_interface = true;
+          else if (cell->neighbor(f)->has_children() == true)
+               {
+              // The neighbor does have children. See if any of the cells
+              // on the other side are elastic
+                 for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
+                   if (cell_is_in_solid_domain (cell->neighbor_child_on_subface(f, sf)))
+                     {
+                   face_is_on_interface = true;
+                           break;
+                     }
+               }
+
+          if (face_is_on_interface)
+           {
+             cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
+             for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
+             if (stokes_fe.face_system_to_component_index(i).first < dim)
+               constraints.add_line (local_face_dof_indices[i]);
+           }
+        }
 @endcode
 
 The call <code>constraints.add_line(t)</code> tells the
index 72f1f0c497e0b71228d1e910712eff637829046b..e9d1ee270a8f16ce8ca4ae9cf5d4791da8edad60 100644 (file)
@@ -367,36 +367,36 @@ namespace Step46
         stokes_fe.dofs_per_face);
       for (const auto &cell : dof_handler.active_cell_iterators())
         if (cell_is_in_fluid_domain(cell))
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
-            if (!cell->at_boundary(f))
-              {
-                bool face_is_on_interface = false;
-
-                if ((cell->neighbor(f)->has_children() == false) &&
-                    (cell_is_in_solid_domain(cell->neighbor(f))))
-                  face_is_on_interface = true;
-                else if (cell->neighbor(f)->has_children() == true)
-                  {
-                    for (unsigned int sf = 0; sf < cell->face(f)->n_children();
-                         ++sf)
-                      if (cell_is_in_solid_domain(
-                            cell->neighbor_child_on_subface(f, sf)))
-                        {
-                          face_is_on_interface = true;
-                          break;
-                        }
-                  }
-
-                if (face_is_on_interface)
-                  {
-                    cell->face(f)->get_dof_indices(local_face_dof_indices, 0);
-                    for (unsigned int i = 0; i < local_face_dof_indices.size();
-                         ++i)
-                      if (stokes_fe.face_system_to_component_index(i).first <
-                          dim)
-                        constraints.add_line(local_face_dof_indices[i]);
-                  }
-              }
+          for (unsigned int face_no : GeometryInfo<dim>::face_indices())
+            {
+              bool face_is_on_interface = false;
+
+              if ((cell->neighbor(face_no)->has_children() == false) &&
+                  (cell_is_in_solid_domain(cell->neighbor(face_no))))
+                face_is_on_interface = true;
+              else if (cell->neighbor(face_no)->has_children() == true)
+                {
+                  for (unsigned int sf = 0;
+                       sf < cell->face(face_no)->n_children();
+                       ++sf)
+                    if (cell_is_in_solid_domain(
+                          cell->neighbor_child_on_subface(face_no, sf)))
+                      {
+                        face_is_on_interface = true;
+                        break;
+                      }
+                }
+
+              if (face_is_on_interface)
+                {
+                  cell->face(face_no)->get_dof_indices(local_face_dof_indices,
+                                                       0);
+                  for (unsigned int i = 0; i < local_face_dof_indices.size();
+                       ++i)
+                    if (stokes_fe.face_system_to_component_index(i).first < dim)
+                      constraints.add_line(local_face_dof_indices[i]);
+                }
+            }
     }
 
     // At the end of all this, we can declare to the constraints object that
@@ -640,122 +640,119 @@ namespace Step46
         // boundary and the potential neighbor behind it is part of the fluid
         // domain. Let's start with these conditions:
         if (cell_is_in_solid_domain(cell))
-          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
-            if (cell->at_boundary(f) == false)
-              {
-                // At this point we know that the current cell is a candidate
-                // for integration and that a neighbor behind face
-                // <code>f</code> exists. There are now three possibilities:
-                //
-                // - The neighbor is at the same refinement level and has no
-                //   children.
-                // - The neighbor has children.
-                // - The neighbor is coarser.
-                //
-                // In all three cases, we are only interested in it if it is
-                // part of the fluid subdomain. So let us start with the first
-                // and simplest case: if the neighbor is at the same level,
-                // has no children, and is a fluid cell, then the two cells
-                // share a boundary that is part of the interface along which
-                // we want to integrate interface terms. All we have to do is
-                // initialize two FEFaceValues object with the current face
-                // and the face of the neighboring cell (note how we find out
-                // which face of the neighboring cell borders on the current
-                // cell) and pass things off to the function that evaluates
-                // the interface terms (the third through fifth arguments to
-                // this function provide it with scratch arrays). The result
-                // is then again copied into the global matrix, using a
-                // function that knows that the DoF indices of rows and
-                // columns of the local matrix result from different cells:
-                if ((cell->neighbor(f)->level() == cell->level()) &&
-                    (cell->neighbor(f)->has_children() == false) &&
-                    cell_is_in_fluid_domain(cell->neighbor(f)))
-                  {
-                    elasticity_fe_face_values.reinit(cell, f);
-                    stokes_fe_face_values.reinit(cell->neighbor(f),
-                                                 cell->neighbor_of_neighbor(f));
-
-                    assemble_interface_term(elasticity_fe_face_values,
-                                            stokes_fe_face_values,
-                                            elasticity_phi,
-                                            stokes_symgrad_phi_u,
-                                            stokes_phi_p,
-                                            local_interface_matrix);
-
-                    cell->neighbor(f)->get_dof_indices(neighbor_dof_indices);
-                    constraints.distribute_local_to_global(
-                      local_interface_matrix,
-                      local_dof_indices,
-                      neighbor_dof_indices,
-                      system_matrix);
-                  }
-
-                // The second case is if the neighbor has further children. In
-                // that case, we have to loop over all the children of the
-                // neighbor to see if they are part of the fluid subdomain. If
-                // they are, then we integrate over the common interface,
-                // which is a face for the neighbor and a subface of the
-                // current cell, requiring us to use an FEFaceValues for the
-                // neighbor and an FESubfaceValues for the current cell:
-                else if ((cell->neighbor(f)->level() == cell->level()) &&
-                         (cell->neighbor(f)->has_children() == true))
-                  {
-                    for (unsigned int subface = 0;
-                         subface < cell->face(f)->n_children();
-                         ++subface)
-                      if (cell_is_in_fluid_domain(
-                            cell->neighbor_child_on_subface(f, subface)))
-                        {
-                          elasticity_fe_subface_values.reinit(cell, f, subface);
-                          stokes_fe_face_values.reinit(
-                            cell->neighbor_child_on_subface(f, subface),
-                            cell->neighbor_of_neighbor(f));
-
-                          assemble_interface_term(elasticity_fe_subface_values,
-                                                  stokes_fe_face_values,
-                                                  elasticity_phi,
-                                                  stokes_symgrad_phi_u,
-                                                  stokes_phi_p,
-                                                  local_interface_matrix);
-
-                          cell->neighbor_child_on_subface(f, subface)
-                            ->get_dof_indices(neighbor_dof_indices);
-                          constraints.distribute_local_to_global(
-                            local_interface_matrix,
-                            local_dof_indices,
-                            neighbor_dof_indices,
-                            system_matrix);
-                        }
-                  }
-
-                // The last option is that the neighbor is coarser. In that
-                // case we have to use an FESubfaceValues object for the
-                // neighbor and a FEFaceValues for the current cell; the rest
-                // is the same as before:
-                else if (cell->neighbor_is_coarser(f) &&
-                         cell_is_in_fluid_domain(cell->neighbor(f)))
-                  {
-                    elasticity_fe_face_values.reinit(cell, f);
-                    stokes_fe_subface_values.reinit(
-                      cell->neighbor(f),
-                      cell->neighbor_of_coarser_neighbor(f).first,
-                      cell->neighbor_of_coarser_neighbor(f).second);
-
-                    assemble_interface_term(elasticity_fe_face_values,
-                                            stokes_fe_subface_values,
-                                            elasticity_phi,
-                                            stokes_symgrad_phi_u,
-                                            stokes_phi_p,
-                                            local_interface_matrix);
-
-                    cell->neighbor(f)->get_dof_indices(neighbor_dof_indices);
-                    constraints.distribute_local_to_global(
-                      local_interface_matrix,
-                      local_dof_indices,
-                      neighbor_dof_indices,
-                      system_matrix);
-                  }
-              }
+          for (unsigned int f : GeometryInfo<dim>::face_indices())
+            {
+              // At this point we know that the current cell is a candidate
+              // for integration and that a neighbor behind face
+              // <code>f</code> exists. There are now three possibilities:
+              //
+              // - The neighbor is at the same refinement level and has no
+              //   children.
+              // - The neighbor has children.
+              // - The neighbor is coarser.
+              //
+              // In all three cases, we are only interested in it if it is
+              // part of the fluid subdomain. So let us start with the first
+              // and simplest case: if the neighbor is at the same level,
+              // has no children, and is a fluid cell, then the two cells
+              // share a boundary that is part of the interface along which
+              // we want to integrate interface terms. All we have to do is
+              // initialize two FEFaceValues object with the current face
+              // and the face of the neighboring cell (note how we find out
+              // which face of the neighboring cell borders on the current
+              // cell) and pass things off to the function that evaluates
+              // the interface terms (the third through fifth arguments to
+              // this function provide it with scratch arrays). The result
+              // is then again copied into the global matrix, using a
+              // function that knows that the DoF indices of rows and
+              // columns of the local matrix result from different cells:
+              if ((cell->neighbor(f)->level() == cell->level()) &&
+                  (cell->neighbor(f)->has_children() == false) &&
+                  cell_is_in_fluid_domain(cell->neighbor(f)))
+                {
+                  elasticity_fe_face_values.reinit(cell, f);
+                  stokes_fe_face_values.reinit(cell->neighbor(f),
+                                               cell->neighbor_of_neighbor(f));
+
+                  assemble_interface_term(elasticity_fe_face_values,
+                                          stokes_fe_face_values,
+                                          elasticity_phi,
+                                          stokes_symgrad_phi_u,
+                                          stokes_phi_p,
+                                          local_interface_matrix);
+
+                  cell->neighbor(f)->get_dof_indices(neighbor_dof_indices);
+                  constraints.distribute_local_to_global(local_interface_matrix,
+                                                         local_dof_indices,
+                                                         neighbor_dof_indices,
+                                                         system_matrix);
+                }
+
+              // The second case is if the neighbor has further children. In
+              // that case, we have to loop over all the children of the
+              // neighbor to see if they are part of the fluid subdomain. If
+              // they are, then we integrate over the common interface,
+              // which is a face for the neighbor and a subface of the
+              // current cell, requiring us to use an FEFaceValues for the
+              // neighbor and an FESubfaceValues for the current cell:
+              else if ((cell->neighbor(f)->level() == cell->level()) &&
+                       (cell->neighbor(f)->has_children() == true))
+                {
+                  for (unsigned int subface = 0;
+                       subface < cell->face(f)->n_children();
+                       ++subface)
+                    if (cell_is_in_fluid_domain(
+                          cell->neighbor_child_on_subface(f, subface)))
+                      {
+                        elasticity_fe_subface_values.reinit(cell, f, subface);
+                        stokes_fe_face_values.reinit(
+                          cell->neighbor_child_on_subface(f, subface),
+                          cell->neighbor_of_neighbor(f));
+
+                        assemble_interface_term(elasticity_fe_subface_values,
+                                                stokes_fe_face_values,
+                                                elasticity_phi,
+                                                stokes_symgrad_phi_u,
+                                                stokes_phi_p,
+                                                local_interface_matrix);
+
+                        cell->neighbor_child_on_subface(f, subface)
+                          ->get_dof_indices(neighbor_dof_indices);
+                        constraints.distribute_local_to_global(
+                          local_interface_matrix,
+                          local_dof_indices,
+                          neighbor_dof_indices,
+                          system_matrix);
+                      }
+                }
+
+              // The last option is that the neighbor is coarser. In that
+              // case we have to use an FESubfaceValues object for the
+              // neighbor and a FEFaceValues for the current cell; the rest
+              // is the same as before:
+              else if (cell->neighbor_is_coarser(f) &&
+                       cell_is_in_fluid_domain(cell->neighbor(f)))
+                {
+                  elasticity_fe_face_values.reinit(cell, f);
+                  stokes_fe_subface_values.reinit(
+                    cell->neighbor(f),
+                    cell->neighbor_of_coarser_neighbor(f).first,
+                    cell->neighbor_of_coarser_neighbor(f).second);
+
+                  assemble_interface_term(elasticity_fe_face_values,
+                                          stokes_fe_subface_values,
+                                          elasticity_phi,
+                                          stokes_symgrad_phi_u,
+                                          stokes_phi_p,
+                                          local_interface_matrix);
+
+                  cell->neighbor(f)->get_dof_indices(neighbor_dof_indices);
+                  constraints.distribute_local_to_global(local_interface_matrix,
+                                                         local_dof_indices,
+                                                         neighbor_dof_indices,
+                                                         system_matrix);
+                }
+            }
       }
   }
 
@@ -952,7 +949,7 @@ namespace Step46
     // encountered when assembling interface terms in
     // <code>assemble_system</code>.
     for (const auto &cell : dof_handler.active_cell_iterators())
-      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      for (unsigned int f : GeometryInfo<dim>::face_indices())
         if (cell_is_in_solid_domain(cell))
           {
             if ((cell->at_boundary(f) == false) &&
index 396cbd7ee2700af5a9095260a5f596d99c10a1d4..f550618ec3762eb6c8814c5390c853a9245bcf89 100644 (file)
@@ -304,11 +304,10 @@ for (auto &cell : triangulation.active_cell_iterators())
   if (cell->center()[1] >= 3.0)
     cell->set_all_manifold_ids(cylinder_id);
 
-for (auto &cell : triangulation.active_cell_iterators())
-  for (unsigned int face_n = 0; face_n < GeometryInfo<3>::faces_per_cell;
-       ++face_n)
+for (const auto &cell : triangulation.active_cell_iterators())
+  for (const auto &face : cell->face_iterators())
     {
-      const Point<3> face_center = cell->face(face_n)->center();
+      const Point<3> face_center = face->center();
       if (std::abs(face_center[0]) < 1.0e-5 &&
           std::abs(face_center[1] - 3.0) < 1.0e-5)
         cell->set_all_manifold_ids(numbers::flat_manifold_id);
index 1178fd47ad44cb365abf3fc6bbbb9306e9b5f887..fa0ac2fa5972e0e73e0f1a52ee9cdea1453d183a 100644 (file)
@@ -547,20 +547,18 @@ namespace Step51
       , fe_support_on_face(GeometryInfo<dim>::faces_per_cell)
       , exact_solution()
     {
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (unsigned int face_no : GeometryInfo<dim>::face_indices())
         for (unsigned int i = 0; i < fe_local.dofs_per_cell; ++i)
           {
-            if (fe_local.has_support_on_face(i, face))
-              fe_local_support_on_face[face].push_back(i);
+            if (fe_local.has_support_on_face(i, face_no))
+              fe_local_support_on_face[face_no].push_back(i);
           }
 
-      for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-           ++face)
+      for (unsigned int face_no : GeometryInfo<dim>::face_indices())
         for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
           {
-            if (fe.has_support_on_face(i, face))
-              fe_support_on_face[face].push_back(i);
+            if (fe.has_support_on_face(i, face_no))
+              fe_support_on_face[face_no].push_back(i);
           }
     }
 
@@ -759,11 +757,10 @@ namespace Step51
     // Face terms are assembled on all faces of all elements. This is in
     // contrast to more traditional DG methods, where each face is only visited
     // once in the assembly procedure.
-    for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-         ++face)
+    for (unsigned int face_no : GeometryInfo<dim>::face_indices())
       {
-        scratch.fe_face_values_local.reinit(loc_cell, face);
-        scratch.fe_face_values.reinit(cell, face);
+        scratch.fe_face_values_local.reinit(loc_cell, face_no);
+        scratch.fe_face_values.reinit(cell, face_no);
 
         // The already obtained $\hat{u}$ values are needed when solving for the
         // local variables.
@@ -792,11 +789,11 @@ namespace Step51
             // We store the non-zero flux and scalar values, making use of the
             // support_on_face information we created in @p ScratchData.
             for (unsigned int k = 0;
-                 k < scratch.fe_local_support_on_face[face].size();
+                 k < scratch.fe_local_support_on_face[face_no].size();
                  ++k)
               {
                 const unsigned int kk =
-                  scratch.fe_local_support_on_face[face][k];
+                  scratch.fe_local_support_on_face[face_no][k];
                 scratch.q_phi[k] =
                   scratch.fe_face_values_local[fluxes].value(kk, q);
                 scratch.u_phi[k] =
@@ -813,29 +810,29 @@ namespace Step51
             if (!task_data.trace_reconstruct)
               {
                 for (unsigned int k = 0;
-                     k < scratch.fe_support_on_face[face].size();
+                     k < scratch.fe_support_on_face[face_no].size();
                      ++k)
                   scratch.tr_phi[k] = scratch.fe_face_values.shape_value(
-                    scratch.fe_support_on_face[face][k], q);
+                    scratch.fe_support_on_face[face_no][k], q);
                 for (unsigned int i = 0;
-                     i < scratch.fe_local_support_on_face[face].size();
+                     i < scratch.fe_local_support_on_face[face_no].size();
                      ++i)
                   for (unsigned int j = 0;
-                       j < scratch.fe_support_on_face[face].size();
+                       j < scratch.fe_support_on_face[face_no].size();
                        ++j)
                     {
                       const unsigned int ii =
-                        scratch.fe_local_support_on_face[face][i];
+                        scratch.fe_local_support_on_face[face_no][i];
                       const unsigned int jj =
-                        scratch.fe_support_on_face[face][j];
+                        scratch.fe_support_on_face[face_no][j];
                       scratch.lf_matrix(ii, jj) +=
                         ((scratch.q_phi[i] * normal +
                           (convection * normal - tau_stab) * scratch.u_phi[i]) *
                          scratch.tr_phi[j]) *
                         JxW;
 
-                      // Note the sign of the face-local matrix.  We negate the
-                      // sign during assembly here so that we can use the
+                      // Note the sign of the face_no-local matrix.  We negate
+                      // the sign during assembly here so that we can use the
                       // FullMatrix::mmult with addition when computing the
                       // Schur complement.
                       scratch.fl_matrix(jj, ii) -=
@@ -846,24 +843,24 @@ namespace Step51
                     }
 
                 for (unsigned int i = 0;
-                     i < scratch.fe_support_on_face[face].size();
+                     i < scratch.fe_support_on_face[face_no].size();
                      ++i)
                   for (unsigned int j = 0;
-                       j < scratch.fe_support_on_face[face].size();
+                       j < scratch.fe_support_on_face[face_no].size();
                        ++j)
                     {
                       const unsigned int ii =
-                        scratch.fe_support_on_face[face][i];
+                        scratch.fe_support_on_face[face_no][i];
                       const unsigned int jj =
-                        scratch.fe_support_on_face[face][j];
+                        scratch.fe_support_on_face[face_no][j];
                       task_data.cell_matrix(ii, jj) +=
                         ((convection * normal - tau_stab) * scratch.tr_phi[i] *
                          scratch.tr_phi[j]) *
                         JxW;
                     }
 
-                if (cell->face(face)->at_boundary() &&
-                    (cell->face(face)->boundary_id() == 1))
+                if (cell->face(face_no)->at_boundary() &&
+                    (cell->face(face_no)->boundary_id() == 1))
                   {
                     const double neumann_value =
                       -scratch.exact_solution.gradient(quadrature_point) *
@@ -871,11 +868,11 @@ namespace Step51
                       convection * normal *
                         scratch.exact_solution.value(quadrature_point);
                     for (unsigned int i = 0;
-                         i < scratch.fe_support_on_face[face].size();
+                         i < scratch.fe_support_on_face[face_no].size();
                          ++i)
                       {
                         const unsigned int ii =
-                          scratch.fe_support_on_face[face][i];
+                          scratch.fe_support_on_face[face_no][i];
                         task_data.cell_vector(ii) +=
                           scratch.tr_phi[i] * neumann_value * JxW;
                       }
@@ -886,16 +883,16 @@ namespace Step51
             // u_h\right>_{\partial \mathcal T}$ to the local matrix. As opposed
             // to the face matrices above, we need it in both assembly stages.
             for (unsigned int i = 0;
-                 i < scratch.fe_local_support_on_face[face].size();
+                 i < scratch.fe_local_support_on_face[face_no].size();
                  ++i)
               for (unsigned int j = 0;
-                   j < scratch.fe_local_support_on_face[face].size();
+                   j < scratch.fe_local_support_on_face[face_no].size();
                    ++j)
                 {
                   const unsigned int ii =
-                    scratch.fe_local_support_on_face[face][i];
+                    scratch.fe_local_support_on_face[face_no][i];
                   const unsigned int jj =
-                    scratch.fe_local_support_on_face[face][j];
+                    scratch.fe_local_support_on_face[face_no][j];
                   scratch.ll_matrix(ii, jj) +=
                     tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
                 }
@@ -908,11 +905,11 @@ namespace Step51
             // since we have moved everything to the other side of the equation.
             if (task_data.trace_reconstruct)
               for (unsigned int i = 0;
-                   i < scratch.fe_local_support_on_face[face].size();
+                   i < scratch.fe_local_support_on_face[face_no].size();
                    ++i)
                 {
                   const unsigned int ii =
-                    scratch.fe_local_support_on_face[face][i];
+                    scratch.fe_local_support_on_face[face_no][i];
                   scratch.l_rhs(ii) -=
                     (scratch.q_phi[i] * normal +
                      scratch.u_phi[i] * (convection * normal - tau_stab)) *
index a856da79f1af756e02c7bdad4dde7f0ef035fd66..5bde5821763ff48d89a8d4883bc08b88b7e289ab 100644 (file)
@@ -772,14 +772,9 @@ namespace Step60
         for (auto &cell : cells)
           {
             cell->set_refine_flag();
-            for (unsigned int face_no = 0;
-                 face_no < GeometryInfo<spacedim>::faces_per_cell;
-                 ++face_no)
+            for (unsigned int face_no : GeometryInfo<spacedim>::face_indices())
               if (!cell->at_boundary(face_no))
-                {
-                  auto neighbor = cell->neighbor(face_no);
-                  neighbor->set_refine_flag();
-                }
+                cell->neighbor(face_no)->set_refine_flag();
           }
         space_grid->execute_coarsening_and_refinement();
       }
index 97c0b84a6605dc40d020a64104a0f8b1f0f53785..1e987d2decbae2fd7ec90b498122280f7eea758d 100644 (file)
@@ -929,12 +929,10 @@ namespace Step61
         // calculate the $L_2$ flux error on the cell by appropriately scaling
         // with face and cell areas and add it to the global error.
         const double cell_area = cell_dgrt->measure();
-        for (unsigned int face_n = 0;
-             face_n < GeometryInfo<dim>::faces_per_cell;
-             ++face_n)
+        for (const auto &face_dgrt : cell_dgrt->face_iterators())
           {
-            const double face_length = cell_dgrt->face(face_n)->measure();
-            fe_face_values_dgrt.reinit(cell_dgrt, face_n);
+            const double face_length = face_dgrt->measure();
+            fe_face_values_dgrt.reinit(cell_dgrt, face_dgrt);
             fe_face_values_dgrt[velocities].get_function_values(
               darcy_velocity, velocity_face_values);
 
index 0bc4b4810c65cc4f40889f265f95a3a87f20a72b..37f4dbb5ddae3d110986639883ef3f5098d5f9ae 100644 (file)
@@ -1043,8 +1043,7 @@ namespace Step9
     // have to clear the array storing the iterators to the active
     // neighbors, of course.
     scratch_data.active_neighbors.clear();
-    for (unsigned int face_n = 0; face_n < GeometryInfo<dim>::faces_per_cell;
-         ++face_n)
+    for (unsigned int face_n : GeometryInfo<dim>::face_indices())
       if (!cell->at_boundary(face_n))
         {
           // First define an abbreviation for the iterator to the face and

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.