template <>
struct InnerPreconditioner<2>
{
- typedef SparseDirectUMFPACK type;
+ using type = SparseDirectUMFPACK;
};
template <>
struct InnerPreconditioner<3>
{
- typedef SparseILU<double> type;
+ using type = SparseILU<double>;
};
@endcode
<h4>Using imhomogeneous constraints for implementing Dirichlet boundary conditions</h4>
-In all the previous tutorial programs, we used the ConstraintMatrix merely
+In all the previous tutorial programs, we used the AffineConstraints object merely
for handling hanging node constraints (with exception of step-11). However,
the class can also be used to implement Dirichlet boundary conditions, as we
will show in this program, by fixing some node values $x_i = b_i$. Note that
these are inhomogeneous constraints, and we have to pay some special
attention to that. The way we are going to implement this is to first read
-in the boundary values into the ConstraintMatrix object by using the call
+in the boundary values into the AffineConstraints object by using the call
@code
VectorTools::interpolate_boundary_values (dof_handler,
very similar to how we were making the list of boundary nodes
before (note that we set Dirichlet conditions only on boundaries with
boundary flag 1). The actual application of the boundary values is then
-handled by the ConstraintMatrix object directly, without any additional
+handled by the AffineConstraints object directly, without any additional
interference.
We could then proceed as before, namely by filling the matrix, and then
simultaneously, since resolving inhomogeneous constraints requires knowledge
about both the matrix entries and the right hand side. For efficiency
reasons, though, we choose another strategy: all the constraints collected
-in the ConstraintMatrix can be resolved on the fly while writing local data
+in the AffineConstraints object can be resolved on the fly while writing local data
into the global matrix, by using the call
@code
constraints.distribute_local_to_global (local_matrix, local_rhs,
implements (inhomogeneous) Dirichlet boundary conditions. That's nice, isn't
it?
-We can conclude that the ConstraintMatrix provides an alternative to using
+We can conclude that the AffineConstraints class provides an alternative to using
MatrixTools::apply_boundary_values for implementing Dirichlet boundary
conditions.
<a name="constraint-matrix">
-<h4>Using ConstraintMatrix for increasing performance</h4>
+<h4>Using AffineConstraints for increasing performance</h4>
</a>
Frequently, a sparse matrix contains a substantial amount of elements that
SparseILU class was very inefficient and has been replaced by one that is
about 10 times faster. Even the vmult function of the SparseILU has been
improved to save about twenty percent of time. Small improvements were
-applied here and there. Moreover, the ConstraintMatrix object has been used
+applied here and there. Moreover, the AffineConstraints object has been used
to eliminate a lot of entries in the sparse matrix that are eventually going
to be zero, see <a href="#constraint-matrix">the section on using advanced
-features of the ConstraintMatrix class</a>.
+features of the AffineConstraints class</a>.
A profile of how many CPU instructions are spent at the various
different places in the program during refinement cycles
// distinguish between them by the use of the spatial dimension as a
// template parameter. See step-4 for details on templates. We are not going
// to create any preconditioner object here, all we do is to create class
- // that holds a local typedef determining the preconditioner class so we can
+ // that holds a local alias determining the preconditioner class so we can
// write our program in a dimension-independent way.
template <int dim>
struct InnerPreconditioner;
template <>
struct InnerPreconditioner<2>
{
- typedef SparseDirectUMFPACK type;
+ using type = SparseDirectUMFPACK;
};
// And the ILU preconditioning in 3D, called by SparseILU:
template <>
struct InnerPreconditioner<3>
{
- typedef SparseILU<double> type;
+ using type = SparseILU<double>;
};
// <code>preconditioner_sparsity_pattern</code>.
// In this example we also use adaptive grid refinement, which is handled
// in analogy to step-6. According to the discussion in the introduction,
- // we are also going to use the ConstraintMatrix for implementing Dirichlet
- // boundary conditions. Hence, we change the name
+ // we are also going to use the AffineConstraints object for implementing
+ // Dirichlet boundary conditions. Hence, we change the name
// <code>hanging_node_constraints</code> into <code>constraints</code>.
template <int dim>
class StokesProblem
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
- ConstraintMatrix constraints;
+ AffineConstraints<double> constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
DoFTools::count_dofs_per_block(dof_handler,
dofs_per_block,
block_component);
- const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
+ const unsigned int n_u = dofs_per_block[0];
+ const unsigned int n_p = dofs_per_block[1];
- std::cout << " Number of active cells: " << triangulation.n_active_cells()
- << std::endl
- << " Number of degrees of freedom: " << dof_handler.n_dofs()
- << " (" << n_u << '+' << n_p << ')' << std::endl;
+ std::cout << " Number of active cells: " //
+ << triangulation.n_active_cells() //
+ << std::endl //
+ << " Number of degrees of freedom: " //
+ << dof_handler.n_dofs() //
+ << " (" << n_u << '+' << n_p << ')' //
+ << std::endl;
// The next task is to allocate a sparsity pattern for the system matrix we
// will create and one for the preconditioner matrix. We could do this in
QGauss<dim> quadrature_formula(degree + 2);
- FEValues<dim> fe_values(fe,
- quadrature_formula,
- update_values | update_quadrature_points |
- update_JxW_values | update_gradients);
+ FEValues<dim> fe_values(fe, //
+ quadrature_formula, //
+ update_values | //
+ update_quadrature_points | //
+ update_JxW_values | //
+ update_gradients);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
std::vector<double> div_phi_u(dofs_per_cell);
std::vector<double> phi_p(dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell != endc; ++cell)
+ for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
local_matrix = 0;
for (unsigned int j = 0; j <= i; ++j)
{
local_matrix(i, j) +=
- (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) -
- div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
- fe_values.JxW(q);
+ (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) //
+ - div_phi_u[i] * phi_p[j] //
+ - phi_p[i] * div_phi_u[j]) //
+ * fe_values.JxW(q);
local_preconditioner_matrix(i, j) +=
(phi_p[i] * phi_p[j]) * fe_values.JxW(q);
const unsigned int component_i =
fe.system_to_component_index(i).first;
- local_rhs(i) += fe_values.shape_value(i, q) *
- rhs_values[q](component_i) * fe_values.JxW(q);
+ local_rhs(i) += fe_values.shape_value(i, q) * //
+ rhs_values[q](component_i) * //
+ fe_values.JxW(q);
}
}
// line of the local matrix contribution.
// Before we can write the local data into the global matrix (and
- // simultaneously use the ConstraintMatrix object to apply Dirichlet
- // boundary conditions and eliminate hanging node constraints, as we
- // discussed in the introduction), we have to be careful about one
+ // simultaneously use the AffineConstraints object to apply
+ // Dirichlet boundary conditions and eliminate hanging node constraints,
+ // as we discussed in the introduction), we have to be careful about one
// thing, though. We have only built half of the local matrices
// because of symmetry, but we're going to save the full matrices
// in order to use the standard functions for solving. This is done
// preconditioner for the velocity-velocity matrix, i.e.,
// <code>block(0,0)</code> in the system matrix. As mentioned above, this
// depends on the spatial dimension. Since the two classes described by
- // the <code>InnerPreconditioner::type</code> typedef have the same
+ // the <code>InnerPreconditioner::type</code> alias have the same
// interface, we do not have to do anything different whether we want to
// use a sparse direct solver or an ILU:
std::cout << " Computing preconditioner..." << std::endl << std::flush;
// pressure field.
constraints.distribute(solution);
- std::cout << " " << solver_control.last_step()
- << " outer CG Schur complement iterations for pressure"
+ std::cout << " " //
+ << solver_control.last_step() //
+ << " outer CG Schur complement iterations for pressure" //
<< std::endl;
}
std::vector<unsigned int> subdivisions(dim, 1);
subdivisions[0] = 4;
- const Point<dim> bottom_left =
- (dim == 2 ? Point<dim>(-2, -1) : Point<dim>(-2, 0, -1));
- const Point<dim> top_right =
- (dim == 2 ? Point<dim>(2, 0) : Point<dim>(2, 1, 0));
+ const Point<dim> bottom_left = (dim == 2 ? //
+ Point<dim>(-2, -1) : //
+ Point<dim>(-2, 0, -1));
+ const Point<dim> top_right = (dim == 2 ? //
+ Point<dim>(2, 0) : //
+ Point<dim>(2, 1, 0));
GridGenerator::subdivided_hyper_rectangle(triangulation,
subdivisions,
// Dirichlet boundary conditions, i.e. to faces that are located at 0 in
// the last coordinate direction. See the example description above for
// details.
- for (typename Triangulation<dim>::active_cell_iterator cell =
- triangulation.begin_active();
- cell != triangulation.end();
- ++cell)
+ 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);
// As first seen in step-6, we cycle over the different refinement levels
// and refine (except for the first cycle), setup the degrees of freedom
// and matrices, assemble, solve and create output:
- for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
+ for (unsigned int refinement_cycle = 0; //
+ refinement_cycle < 6; //
++refinement_cycle)
{
std::cout << "Refinement cycle " << refinement_cycle << std::endl;