// 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
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
//
// 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);
// 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
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);
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);
// 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);
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
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
// 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);
+ }
+ }
}
}
// 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) &&
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);
, 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);
}
}
// 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.
// 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] =
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) -=
}
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) *
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;
}
// 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;
}
// 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)) *
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();
}
// 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);
// 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
projects / dealii.git / commitdiff