From: Wolfgang Bangerth Date: Mon, 3 Feb 2020 01:47:22 +0000 (-0700) Subject: Use cell->face_iterators() or GeometryInfo::face_indices() in the tutorials. X-Git-Tag: v9.2.0-rc1~576^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=448b0a5964fa12feb64adae8a77d7815e33f1a78;p=dealii.git Use cell->face_iterators() or GeometryInfo::face_indices() in the tutorials. --- diff --git a/examples/step-14/step-14.cc b/examples/step-14/step-14.cc index 56068529fd..2813ba7893 100644 --- a/examples/step-14/step-14.cc +++ b/examples/step-14/step-14.cc @@ -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::faces_per_cell; - ++face_no) + for (unsigned int face_no : GeometryInfo::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 diff --git a/examples/step-20/doc/intro.dox b/examples/step-20/doc/intro.dox index 22b5ea333e..f89ce1546e 100644 --- a/examples/step-20/doc/intro.dox +++ b/examples/step-20/doc/intro.dox @@ -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::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; qface_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 diff --git a/examples/step-21/step-21.cc b/examples/step-21/step-21.cc index aaef851b44..4104101f0e 100644 --- a/examples/step-21/step-21.cc +++ b/examples/step-21/step-21.cc @@ -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::faces_per_cell; - ++face_no) + for (unsigned int face_no : GeometryInfo::face_indices()) { fe_face_values.reinit(cell, face_no); diff --git a/examples/step-29/step-29.cc b/examples/step-29/step-29.cc index 3146ed1307..0b196c4306 100644 --- a/examples/step-29/step-29.cc +++ b/examples/step-29/step-29.cc @@ -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::faces_per_cell; - ++face) - if (cell->face(face)->at_boundary() && - (cell->face(face)->boundary_id() == 0)) + for (unsigned int face_no : GeometryInfo::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 diff --git a/examples/step-30/step-30.cc b/examples/step-30/step-30.cc index d2ae115271..6d78192309 100644 --- a/examples/step-30/step-30.cc +++ b/examples/step-30/step-30.cc @@ -451,9 +451,7 @@ namespace Step30 cell->get_dof_indices(dofs); - for (unsigned int face_no = 0; - face_no < GeometryInfo::faces_per_cell; - ++face_no) + for (unsigned int face_no : GeometryInfo::face_indices()) { const auto face = cell->face(face_no); @@ -726,9 +724,7 @@ namespace Step30 Point jump; Point area; - for (unsigned int face_no = 0; - face_no < GeometryInfo::faces_per_cell; - ++face_no) + for (unsigned int face_no : GeometryInfo::face_indices()) { const auto face = cell->face(face_no); diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc index 35cb2cd473..1f5d3bf6c6 100644 --- a/examples/step-33/step-33.cc +++ b/examples/step-33/step-33.cc @@ -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::faces_per_cell; - ++face_no) + for (unsigned int face_no : GeometryInfo::face_indices()) if (cell->at_boundary(face_no)) { fe_v_face.reinit(cell, face_no); diff --git a/examples/step-46/doc/intro.dox b/examples/step-46/doc/intro.dox index a0fa0fdac7..b77ed76f19 100644 --- a/examples/step-46/doc/intro.dox +++ b/examples/step-46/doc/intro.dox @@ -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 local_face_dof_indices (stokes_fe.dofs_per_face); - for (typename hp::DoFHandler::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::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; sfface(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 (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::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; sfface(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; iconstraints.add_line(t) tells the diff --git a/examples/step-46/step-46.cc b/examples/step-46/step-46.cc index 72f1f0c497..e9d1ee270a 100644 --- a/examples/step-46/step-46.cc +++ b/examples/step-46/step-46.cc @@ -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::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::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::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 - // f 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::face_indices()) + { + // At this point we know that the current cell is a candidate + // for integration and that a neighbor behind face + // f 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 // assemble_system. for (const auto &cell : dof_handler.active_cell_iterators()) - for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + for (unsigned int f : GeometryInfo::face_indices()) if (cell_is_in_solid_domain(cell)) { if ((cell->at_boundary(f) == false) && diff --git a/examples/step-49/doc/results.dox b/examples/step-49/doc/results.dox index 396cbd7ee2..f550618ec3 100644 --- a/examples/step-49/doc/results.dox +++ b/examples/step-49/doc/results.dox @@ -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); diff --git a/examples/step-51/step-51.cc b/examples/step-51/step-51.cc index 1178fd47ad..fa0ac2fa59 100644 --- a/examples/step-51/step-51.cc +++ b/examples/step-51/step-51.cc @@ -547,20 +547,18 @@ namespace Step51 , fe_support_on_face(GeometryInfo::faces_per_cell) , exact_solution() { - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; - ++face) + for (unsigned int face_no : GeometryInfo::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::faces_per_cell; - ++face) + for (unsigned int face_no : GeometryInfo::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::faces_per_cell; - ++face) + for (unsigned int face_no : GeometryInfo::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)) * diff --git a/examples/step-60/step-60.cc b/examples/step-60/step-60.cc index a856da79f1..5bde582176 100644 --- a/examples/step-60/step-60.cc +++ b/examples/step-60/step-60.cc @@ -772,14 +772,9 @@ namespace Step60 for (auto &cell : cells) { cell->set_refine_flag(); - for (unsigned int face_no = 0; - face_no < GeometryInfo::faces_per_cell; - ++face_no) + for (unsigned int face_no : GeometryInfo::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(); } diff --git a/examples/step-61/step-61.cc b/examples/step-61/step-61.cc index 97c0b84a66..1e987d2dec 100644 --- a/examples/step-61/step-61.cc +++ b/examples/step-61/step-61.cc @@ -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::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); diff --git a/examples/step-9/step-9.cc b/examples/step-9/step-9.cc index 0bc4b4810c..37f4dbb5dd 100644 --- a/examples/step-9/step-9.cc +++ b/examples/step-9/step-9.cc @@ -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::faces_per_cell; - ++face_n) + for (unsigned int face_n : GeometryInfo::face_indices()) if (!cell->at_boundary(face_n)) { // First define an abbreviation for the iterator to the face and