From d40fe5775c01d0318c559b8160506dce64debb9f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 7 Jan 2020 14:40:59 -0700 Subject: [PATCH] Convert step-61 to use face loops with cell->face_iterators(). This functionality was previously implemented in e2d7a2ff171a111480307d97dce2aa4db063a953 by @masterleinad but had crossed paths with 3f96459cacd673c231357e87c15e95315857e749 (in #8999). --- examples/step-61/step-61.cc | 37 +++++++++++++++++-------------------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/examples/step-61/step-61.cc b/examples/step-61/step-61.cc index ff77d4ecd1..97c0b84a66 100644 --- a/examples/step-61/step-61.cc +++ b/examples/step-61/step-61.cc @@ -442,8 +442,8 @@ namespace Step61 // We need FEValuesExtractors to access the @p interior and // @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 pressure_interior(0); + const FEValuesExtractors::Scalar pressure_face(1); // This finally gets us in position to loop over all cells. On // each cell, we will first calculate the various cell matrices @@ -505,7 +505,8 @@ namespace Step61 fe_values_dgrt[velocities].divergence(i, q); for (unsigned int j = 0; j < dofs_per_cell; ++j) { - const double phi_j_interior = fe_values[interior].value(j, q); + const double phi_j_interior = + fe_values[pressure_interior].value(j, q); cell_matrix_G(i, j) -= (div_v_i * phi_j_interior * fe_values.JxW(q)); @@ -518,12 +519,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::faces_per_cell; - ++face_n) + for (const auto &face : cell->face_iterators()) { - fe_face_values.reinit(cell, face_n); - fe_face_values_dgrt.reinit(cell_dgrt, face_n); + fe_face_values.reinit(cell, face); + fe_face_values_dgrt.reinit(cell_dgrt, face); for (unsigned int q = 0; q < n_face_q_points; ++q) { @@ -536,7 +535,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[pressure_face].value(j, q); cell_matrix_G(i, j) += ((v_i * normal) * phi_j_face * fe_face_values.JxW(q)); @@ -582,7 +581,7 @@ namespace Step61 for (unsigned int q = 0; q < n_q_points; ++q) for (unsigned int i = 0; i < dofs_per_cell; ++i) { - cell_rhs(i) += (fe_values[interior].value(i, q) * + cell_rhs(i) += (fe_values[pressure_interior].value(i, q) * right_hand_side_values[q] * fe_values.JxW(q)); } @@ -694,9 +693,8 @@ namespace Step61 std::vector> coefficient_values(n_q_points_dgrt); const FEValuesExtractors::Vector velocities(0); - const FEValuesExtractors::Scalar pressure(dim); - const FEValuesExtractors::Scalar interior(0); - const FEValuesExtractors::Scalar face(1); + const FEValuesExtractors::Scalar pressure_interior(0); + const FEValuesExtractors::Scalar pressure_face(1); // In the introduction, we explained how to calculate the numerical velocity // on the cell. We need the pressure solution values on each cell, @@ -760,19 +758,18 @@ namespace Step61 fe_values_dgrt[velocities].divergence(i, q); for (unsigned int j = 0; j < dofs_per_cell; ++j) { - const double phi_j_interior = fe_values[interior].value(j, q); + const double phi_j_interior = + fe_values[pressure_interior].value(j, q); cell_matrix_G(i, j) -= (div_v_i * phi_j_interior * fe_values.JxW(q)); } } - for (unsigned int face_n = 0; - face_n < GeometryInfo::faces_per_cell; - ++face_n) + for (const auto &face : cell->face_iterators()) { - fe_face_values.reinit(cell, face_n); - fe_face_values_dgrt.reinit(cell_dgrt, face_n); + fe_face_values.reinit(cell, face); + fe_face_values_dgrt.reinit(cell_dgrt, face); for (unsigned int q = 0; q < n_face_q_points; ++q) { @@ -785,7 +782,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[pressure_face].value(j, q); cell_matrix_G(i, j) += ((v_i * normal) * phi_j_face * fe_face_values.JxW(q)); -- 2.39.5