From 40072e92c38338791163a9df4d527764063e1ff0 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 21 Jun 2018 14:54:49 +0200 Subject: [PATCH] examples/step-22: Update indenting and modernize --- examples/step-22/doc/intro.dox | 20 ++++---- examples/step-22/step-22.cc | 84 ++++++++++++++++++---------------- 2 files changed, 55 insertions(+), 49 deletions(-) diff --git a/examples/step-22/doc/intro.dox b/examples/step-22/doc/intro.dox index ead7896cc7..74a11a1dde 100644 --- a/examples/step-22/doc/intro.dox +++ b/examples/step-22/doc/intro.dox @@ -554,13 +554,13 @@ struct InnerPreconditioner; template <> struct InnerPreconditioner<2> { - typedef SparseDirectUMFPACK type; + using type = SparseDirectUMFPACK; }; template <> struct InnerPreconditioner<3> { - typedef SparseILU type; + using type = SparseILU; }; @endcode @@ -621,13 +621,13 @@ left by the outward motion of material at this location.

Using imhomogeneous constraints for implementing Dirichlet boundary conditions

-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, @@ -639,7 +639,7 @@ in the boundary values into the ConstraintMatrix object by using the call 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 @@ -652,7 +652,7 @@ Note that we call this on the system matrix and system right hand side 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, @@ -667,13 +667,13 @@ side, it distributes the hanging node constraints and additionally 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. -

Using ConstraintMatrix for increasing performance

+

Using AffineConstraints for increasing performance

 Frequently, a sparse matrix contains a substantial amount of elements that @@ -752,10 +752,10 @@ consumption. Likewise, the implementation of the decomposition step in the 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 the section on using advanced -features of the ConstraintMatrix class. +features of the AffineConstraints class. A profile of how many CPU instructions are spent at the various different places in the program during refinement cycles diff --git a/examples/step-22/step-22.cc b/examples/step-22/step-22.cc index dd023a8c73..d3d4926b01 100644 --- a/examples/step-22/step-22.cc +++ b/examples/step-22/step-22.cc @@ -79,7 +79,7 @@ namespace Step22 // 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 struct InnerPreconditioner; @@ -88,14 +88,14 @@ namespace Step22 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 type; + using type = SparseILU; }; @@ -108,8 +108,8 @@ namespace Step22 // preconditioner_sparsity_pattern. // 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 // hanging_node_constraints into constraints. template class StokesProblem @@ -131,7 +131,7 @@ namespace Step22 FESystem fe; DoFHandler dof_handler; - ConstraintMatrix constraints; + AffineConstraints constraints; BlockSparsityPattern sparsity_pattern; BlockSparseMatrix system_matrix; @@ -498,12 +498,16 @@ namespace Step22 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 @@ -615,10 +619,12 @@ namespace Step22 QGauss quadrature_formula(degree + 2); - FEValues fe_values(fe, - quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values | update_gradients); + FEValues 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; @@ -666,10 +672,7 @@ namespace Step22 std::vector div_phi_u(dofs_per_cell); std::vector phi_p(dofs_per_cell); - typename DoFHandler::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; @@ -694,9 +697,10 @@ namespace Step22 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); @@ -716,8 +720,9 @@ namespace Step22 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); } } @@ -726,9 +731,9 @@ namespace Step22 // 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 @@ -757,7 +762,7 @@ namespace Step22 // preconditioner for the velocity-velocity matrix, i.e., // block(0,0) in the system matrix. As mentioned above, this // depends on the spatial dimension. Since the two classes described by - // the InnerPreconditioner::type typedef have the same + // the InnerPreconditioner::type 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; @@ -847,8 +852,9 @@ namespace Step22 // 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; } @@ -976,10 +982,12 @@ namespace Step22 std::vector subdivisions(dim, 1); subdivisions[0] = 4; - const Point bottom_left = - (dim == 2 ? Point(-2, -1) : Point(-2, 0, -1)); - const Point top_right = - (dim == 2 ? Point(2, 0) : Point(2, 1, 0)); + const Point bottom_left = (dim == 2 ? // + Point(-2, -1) : // + Point(-2, 0, -1)); + const Point top_right = (dim == 2 ? // + Point(2, 0) : // + Point(2, 1, 0)); GridGenerator::subdivided_hyper_rectangle(triangulation, subdivisions, @@ -991,10 +999,7 @@ namespace Step22 // 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::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::faces_per_cell; ++f) if (cell->face(f)->center()[dim - 1] == 0) cell->face(f)->set_all_boundary_ids(1); @@ -1008,7 +1013,8 @@ namespace Step22 // 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; -- 2.39.5