// 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)
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,
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;
}
{
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:
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);
// 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);
// 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
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)
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
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);
}
}
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;
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)
}
}
- 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);
}
}
- 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);
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);
}
}
- 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);
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);
// +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);
}
}
}
// 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);
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())
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()]++;
}
}
// 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}
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
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;
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
// @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
// 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)
{
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));
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;
}
}
- 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)
{
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));
// 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)
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);
// <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
// 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.
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
// 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: