]> https://gitweb.dealii.org/ - dealii.git/commitdiff
port step-16 to mesh_loop
authorTimo Heister <timo.heister@gmail.com>
Tue, 24 Jul 2018 15:42:57 +0000 (17:42 +0200)
committerTimo Heister <timo.heister@gmail.com>
Wed, 5 Sep 2018 15:00:55 +0000 (11:00 -0400)
examples/step-16/step-16.cc

index 544633329188db68fb4bee2cc56da57decd253df..f28d688d3c0d0d6eb0d3239cbc7fd2cfd8786f3a 100644 (file)
@@ -16,6 +16,7 @@
  * Authors: Guido Kanschat, University of Heidelberg, 2003
  *          Baerbel Janssen, University of Heidelberg, 2010
  *          Wolfgang Bangerth, Texas A&M University, 2010
+ *          Timo Heister, Clemson University, 2018
  */
 
 
 #include <deal.II/multigrid/mg_smoother.h>
 #include <deal.II/multigrid/mg_matrix.h>
 
-// Finally we include the MeshWorker framework. This framework through its
-// function loop() and integration_loop(), automates loops over cells and
-// assembling of data into vectors, matrices, etc. It obeys constraints
-// automatically. Since we have to build several matrices and have to be aware
-// of several sets of constraints, this will save us a lot of headache.
-#include <deal.II/meshworker/dof_info.h>
-#include <deal.II/meshworker/integration_info.h>
-#include <deal.II/meshworker/simple.h>
-#include <deal.II/meshworker/output.h>
-#include <deal.II/meshworker/loop.h>
-
-// In order to save effort, we use the pre-implemented Laplacian found in
-#include <deal.II/integrators/laplace.h>
-#include <deal.II/integrators/l2.h>
+// We will be using MeshWorker::mesh_loop to loop over the cells, so include it
+// here:
+#include <deal.II/meshworker/mesh_loop.h>
+
 
 // This is C++:
 #include <iostream>
@@ -91,85 +82,59 @@ using namespace dealii;
 
 namespace Step16
 {
-  // @sect3{The integrator on each cell}
-
-  // The MeshWorker::integration_loop() expects a class that provides functions
-  // for integration on cells and boundary and interior faces. This is done by
-  // the following class. In the constructor, we tell the loop that cell
-  // integrals should be computed (the 'true'), but integrals should not be
-  // computed on boundary and interior faces (the two 'false'). Accordingly, we
-  // only need a cell function, but none for the faces.
+  // @sect3{The Scratch and Copy objects}
+  //
+  // We use MeshWorker::mesh_loop() to assemble our matrices. For this, we need
+  // a ScratchData object to store temporary data on each cell (this is just the
+  // FEValues object) and a CopyData object that will contain the output of each
+  // cell assembly.
   template <int dim>
-  class LaplaceIntegrator : public MeshWorker::LocalIntegrator<dim>
+  struct ScratchData
   {
-  public:
-    LaplaceIntegrator();
-    virtual void cell(MeshWorker::DoFInfo<dim> &        dinfo,
-                      MeshWorker::IntegrationInfo<dim> &info) const override;
+    ScratchData(const Mapping<dim> &      mapping,
+                const FiniteElement<dim> &fe,
+                const unsigned int        quadrature_degree,
+                const UpdateFlags         update_flags)
+      : fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags)
+    {}
+
+    ScratchData(const ScratchData<dim> &scratch_data)
+      : fe_values(scratch_data.fe_values.get_mapping(),
+                  scratch_data.fe_values.get_fe(),
+                  scratch_data.fe_values.get_quadrature(),
+                  scratch_data.fe_values.get_update_flags())
+    {}
+
+    FEValues<dim> fe_values;
   };
 
-
-  template <int dim>
-  LaplaceIntegrator<dim>::LaplaceIntegrator()
-    : MeshWorker::LocalIntegrator<dim>(true, false, false)
-  {}
-
-
-  // Next the actual integrator on each cell. We solve a Poisson problem with a
-  // coefficient one in the right half plane and one tenth in the left
-  // half plane.
-
-  // The MeshWorker::LocalResults base class of MeshWorker::DoFInfo contains
-  // objects that can be filled in this local integrator. How many objects are
-  // created is determined inside the MeshWorker framework by the assembler
-  // class. Here, we test for instance that one matrix is required
-  // (MeshWorker::LocalResults::n_matrices()). The matrices are accessed through
-  // MeshWorker::LocalResults::matrix(), which takes the number of the matrix as
-  // its first argument. The second argument is only used for integrals over
-  // faces when there are two matrices for each test function used. Then, a
-  // second matrix with indicator 'true' would exist with the same index.
-
-  // MeshWorker::IntegrationInfo provides one or several FEValues objects, which
-  // below are used by LocalIntegrators::Laplace::cell_matrix() or
-  // LocalIntegrators::L2::L2(). Since we are assembling only a single PDE,
-  // there is also only one of these objects with index zero.
-
-  // In addition, we note that this integrator serves to compute the matrices
-  // for the multilevel preconditioner as well as the matrix and the right hand
-  // side for the global system. Since the assembler for a system requires an
-  // additional vector, MeshWorker::LocalResults::n_vectors() is returning a
-  // nonzero value. Accordingly, we fill a right hand side vector at the end of
-  // this function. Since LocalResults can deal with several BlockVector
-  // objects, but we are again in the simplest case here, we enter the
-  // information into block zero of vector zero.
-  template <int dim>
-  void
-  LaplaceIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &        dinfo,
-                               MeshWorker::IntegrationInfo<dim> &info) const
+  struct CopyData
   {
-    AssertDimension(dinfo.n_matrices(), 1);
-    const double coefficient = (dinfo.cell->center()(0) > 0.) ? .1 : 1.;
+    unsigned int                         level;
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
 
-    LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0, false).matrix,
-                                           info.fe_values(0),
-                                           coefficient);
-
-    if (dinfo.n_vectors() > 0)
-      {
-        std::vector<double> rhs(info.fe_values(0).n_quadrature_points, 1.);
-        LocalIntegrators::L2::L2(dinfo.vector(0).block(0),
-                                 info.fe_values(0),
-                                 rhs);
-      }
-  }
+    template <class Iterator>
+    void reinit(const Iterator &cell, unsigned int dofs_per_cell)
+    {
+      cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+      cell_rhs.reinit(dofs_per_cell);
 
+      local_dof_indices.resize(dofs_per_cell);
+      cell->get_active_or_mg_dof_indices(local_dof_indices);
+      level = cell->level();
+    }
+  };
 
   // @sect3{The <code>LaplaceProblem</code> class template}
 
-  // This main class is basically the same class as in step-6. As far as
-  // member functions is concerned, the only addition is the
-  // <code>assemble_multigrid</code> function that assembles the matrices that
-  // correspond to the discrete operators on intermediate levels:
+  // This main class is similar to the same class in step-6. As far as
+  // member functions is concerned, the only additions are:
+  // - The <code>assemble_multigrid</code> function that assembles the matrices
+  // that correspond to the discrete operators on intermediate levels.
+  // - The <code>cell_worker</code> function that assembles our PDE on a single
+  // cell.
   template <int dim>
   class LaplaceProblem
   {
@@ -178,6 +143,11 @@ namespace Step16
     void run();
 
   private:
+    template <class Iterator>
+    void cell_worker(const Iterator &  cell,
+                     ScratchData<dim> &scratch_data,
+                     CopyData &        copy_data);
+
     void setup_system();
     void assemble_system();
     void assemble_multigrid();
@@ -200,9 +170,9 @@ namespace Step16
     const unsigned int degree;
 
     // The following members are the essential data structures for the multigrid
-    // method. The first two represent the sparsity patterns and the matrices on
-    // individual levels of the multilevel hierarchy, very much like the objects
-    // for the global mesh above.
+    // method. The first four represent the sparsity patterns and the matrices
+    // on individual levels of the multilevel hierarchy, very much like the
+    // objects for the global mesh above.
     //
     // Then we have two new matrices only needed for multigrid methods with
     // local smoothing on adaptive meshes. They convey data between the interior
@@ -213,10 +183,11 @@ namespace Step16
     // level and information about indices lying on a refinement edge between
     // two different refinement levels. It thus serves a similar purpose as
     // AffineConstraints, but on each level.
-    MGLevelObject<SparsityPattern>      mg_sparsity_patterns;
+    MGLevelObject<SparsityPattern> mg_sparsity_patterns;
+    MGLevelObject<SparsityPattern> mg_interface_sparsity_patterns;
+
     MGLevelObject<SparseMatrix<double>> mg_matrices;
-    MGLevelObject<SparseMatrix<double>> mg_interface_in;
-    MGLevelObject<SparseMatrix<double>> mg_interface_out;
+    MGLevelObject<SparseMatrix<double>> mg_interface_matrices;
     MGConstrainedDoFs                   mg_constrained_dofs;
   };
 
@@ -257,15 +228,13 @@ namespace Step16
     dof_handler.distribute_dofs(fe);
     dof_handler.distribute_mg_dofs();
 
-    deallog << "   Number of degrees of freedom: " << dof_handler.n_dofs()
-            << " (by level: ";
+    std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << " (by level: ";
     for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
-      deallog << dof_handler.n_dofs(level)
-              << (level == triangulation.n_levels() - 1 ? ")" : ", ");
-    deallog << std::endl;
+      std::cout << dof_handler.n_dofs(level)
+                << (level == triangulation.n_levels() - 1 ? ")" : ", ");
+    std::cout << std::endl;
 
-    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
-    DoFTools::make_sparsity_pattern(dof_handler, dsp);
 
     solution.reinit(dof_handler.n_dofs());
     system_rhs.reinit(dof_handler.n_dofs());
@@ -278,18 +247,20 @@ namespace Step16
     const std::map<types::boundary_id, const Function<dim> *>
       dirichlet_boundary_functions = {
         {types::boundary_id(0), &homogeneous_dirichlet_bc}};
-    VectorTools::interpolate_boundary_values(
-      static_cast<const DoFHandler<dim> &>(dof_handler),
-      dirichlet_boundary_functions,
-      constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             dirichlet_boundary_functions,
+                                             constraints);
     constraints.close();
-    constraints.condense(dsp);
-    sparsity_pattern.copy_from(dsp);
+
+    {
+      DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+      DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
+      sparsity_pattern.copy_from(dsp);
+    }
     system_matrix.reinit(sparsity_pattern);
 
     // The multigrid constraints have to be initialized. They need to know
-    // about the boundary values as well, so we pass the
-    // <code>dirichlet_boundary</code> here as well.
+    // where Dirichlet boundary conditions are prescribed.
     mg_constrained_dofs.clear();
     mg_constrained_dofs.initialize(dof_handler);
     mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
@@ -306,143 +277,218 @@ namespace Step16
     // their SparsityPattern before the can be destroyed upon resizing.
     const unsigned int n_levels = triangulation.n_levels();
 
-    mg_interface_in.resize(0, n_levels - 1);
-    mg_interface_in.clear_elements();
-    mg_interface_out.resize(0, n_levels - 1);
-    mg_interface_out.clear_elements();
+    mg_interface_matrices.resize(0, n_levels - 1);
     mg_matrices.resize(0, n_levels - 1);
-    mg_matrices.clear_elements();
     mg_sparsity_patterns.resize(0, n_levels - 1);
+    mg_interface_sparsity_patterns.resize(0, n_levels - 1);
 
     // Now, we have to provide a matrix on each level. To this end, we first use
     // the MGTools::make_sparsity_pattern function to generate a preliminary
     // compressed sparsity pattern on each level (see the @ref Sparsity module
     // for more information on this topic) and then copy it over to the one we
-    // really want. The next step is to initialize both kinds of level matrices
-    // with these sparsity patterns.
+    // really want. The next step is to initialize the interface matrices with
+    // the fitting sparsity pattern.
     //
     // It may be worth pointing out that the interface matrices only have
     // entries for degrees of freedom that sit at or next to the interface
     // between coarser and finer levels of the mesh. They are therefore even
     // sparser than the matrices on the individual levels of our multigrid
-    // hierarchy. If we were more concerned about memory usage (and possibly the
-    // speed with which we can multiply with these matrices), we should use
-    // separate and different sparsity patterns for these two kinds of matrices.
+    // hierarchy. Therefore, we use a function specifically build for this
+    // purpose to generate it.
     for (unsigned int level = 0; level < n_levels; ++level)
       {
-        DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
-                                   dof_handler.n_dofs(level));
-        MGTools::make_sparsity_pattern(dof_handler, dsp, level);
+        {
+          DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
+                                     dof_handler.n_dofs(level));
+          MGTools::make_sparsity_pattern(dof_handler, dsp, level);
+
+          mg_sparsity_patterns[level].copy_from(dsp);
+          mg_matrices[level].reinit(mg_sparsity_patterns[level]);
+        }
+        {
+          DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
+                                     dof_handler.n_dofs(level));
+          MGTools::make_interface_sparsity_pattern(dof_handler,
+                                                   mg_constrained_dofs,
+                                                   dsp,
+                                                   level);
+          mg_interface_sparsity_patterns[level].copy_from(dsp);
+          mg_interface_matrices[level].reinit(
+            mg_interface_sparsity_patterns[level]);
+        }
+      }
+  }
+
+
+  // @sect4{LaplaceProblem::cell_worker}
 
-        mg_sparsity_patterns[level].copy_from(dsp);
+  // The cell_worker function is used to assemble the matrix and right-hand side
+  // on the given cell. This function is used for the active cells to generate
+  // the system_matrix and on each level to build the level matrices.
+  //
+  // Note that we also assemble a right-hand side when called from
+  // assemble_multigrid() even though it is not used.
+  template <int dim>
+  template <class Iterator>
+  void LaplaceProblem<dim>::cell_worker(const Iterator &  cell,
+                                        ScratchData<dim> &scratch_data,
+                                        CopyData &        copy_data)
+  {
+    FEValues<dim> &fe_values = scratch_data.fe_values;
+    fe_values.reinit(cell);
+
+    const unsigned int dofs_per_cell = fe_values.get_fe().dofs_per_cell;
+    const unsigned int n_q_points    = fe_values.get_quadrature().size();
+
+    copy_data.reinit(cell, dofs_per_cell);
+
+    const std::vector<double> &JxW = fe_values.get_JxW_values();
+
+    for (unsigned int q = 0; q < n_q_points; ++q)
+      {
+        const double coefficient =
+          (fe_values.get_quadrature_points()[q][0] < 0.0) ? 1.0 : 0.1;
+        //(cell->center().square() < 0.5 * 0.5) ? 10.0:1.0;
 
-        mg_matrices[level].reinit(mg_sparsity_patterns[level]);
-        mg_interface_in[level].reinit(mg_sparsity_patterns[level]);
-        mg_interface_out[level].reinit(mg_sparsity_patterns[level]);
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              {
+                copy_data.cell_matrix(i, j) +=
+                  coefficient *
+                  (fe_values.shape_grad(i, q) * fe_values.shape_grad(j, q)) *
+                  JxW[q];
+              }
+            copy_data.cell_rhs(i) += 1.0 * fe_values.shape_value(i, q) * JxW[q];
+          }
       }
   }
 
 
+
   // @sect4{LaplaceProblem::assemble_system}
 
-  // The following function assembles the linear system on the finest level of
-  // the mesh. Since we want to reuse the code here for the level assembling
-  // below, we use the local integrator class LaplaceIntegrator and leave the
-  // loops to the MeshWorker framework. Thus, this function first sets up the
-  // objects necessary for this framework, namely
-  //   - a MeshWorker::IntegrationInfoBox object, which will provide all the
-  //     required data in quadrature points on the cell. This object can be seen
-  //     as an extension of FEValues, providing a lot more useful information,
-  //   - a MeshWorker::DoFInfo object, which on the one hand side extends the
-  //     functionality of cell iterators, but also provides space for return
-  //     values in its base class LocalResults,
-  //   - an assembler, in this case for the whole system. The term 'simple' here
-  //     refers to the fact that the global system does not have a block
-  //     structure,
-  //   - the local integrator, which implements the actual forms.
-  //
-  // After the loop has combined all of these into a matrix and a right hand
-  // side, there is one thing left to do: the assemblers leave matrix rows and
-  // columns of constrained degrees of freedom untouched. Therefore, we put a
-  // one on the diagonal to make the whole system well posed. The value one, or
-  // any fixed value has the advantage, that its effect on the spectrum of the
-  // matrix is easily understood. Since the corresponding eigenvectors form an
-  // invariant subspace, the value chosen does not affect the convergence of
-  // Krylov space solvers.
+  // The following function assembles the linear system on the active cells of
+  // the mesh. For this, we pass two lambda functions to the mesh_loop()
+  // function. The cell_worker function redirects to the class member function
+  // of the same name, while the copier is specific to this function and copies
+  // local matrix and vector to the corresponding global ones using the
+  // constraints.
   template <int dim>
   void LaplaceProblem<dim>::assemble_system()
   {
-    MappingQ1<dim>                      mapping;
-    MeshWorker::IntegrationInfoBox<dim> info_box;
-    UpdateFlags                         update_flags =
-      update_values | update_gradients | update_hessians;
-    info_box.add_update_flags_all(update_flags);
-    info_box.initialize(fe, mapping);
-
-    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
-    MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double>>
-      assembler;
-    assembler.initialize(constraints);
-    assembler.initialize(system_matrix, system_rhs);
-
-    LaplaceIntegrator<dim> matrix_integrator;
-    MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
-                                           dof_handler.end(),
-                                           dof_info,
-                                           info_box,
-                                           matrix_integrator,
-                                           assembler);
-
-    for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
-      if (constraints.is_constrained(i))
-        system_matrix.set(i, i, 1.);
+    MappingQ1<dim> mapping;
+
+    typedef decltype(dof_handler.begin_active()) Iterator;
+
+    auto cell_worker = [&](const Iterator &  cell,
+                           ScratchData<dim> &scratch_data,
+                           CopyData &        copy_data) {
+      this->cell_worker(cell, scratch_data, copy_data);
+    };
+
+    auto copier = [&](const CopyData &cd) {
+      this->constraints.distribute_local_to_global(cd.cell_matrix,
+                                                   cd.cell_rhs,
+                                                   cd.local_dof_indices,
+                                                   system_matrix,
+                                                   system_rhs);
+    };
+
+    const unsigned int n_gauss_points = degree + 1;
+
+    ScratchData<dim> scratch_data(mapping,
+                                  fe,
+                                  n_gauss_points,
+                                  update_values | update_gradients |
+                                    update_JxW_values |
+                                    update_quadrature_points);
+
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          CopyData(),
+                          MeshWorker::assemble_own_cells);
   }
 
 
   // @sect4{LaplaceProblem::assemble_multigrid}
 
-  // The next function is the one that builds the linear operators (matrices)
+  // The next function is the one that builds the matrices
   // that define the multigrid method on each level of the mesh. The integration
   // core is the same as above, but the loop below will go over all existing
   // cells instead of just the active ones, and the results must be entered into
   // the correct level matrices. Fortunately, MeshWorker hides most of that from
   // us, and thus the difference between this function and the previous lies
   // only in the setup of the assembler and the different iterators in the loop.
-  // Also, fixing up the matrices in the end is a little more complicated.
+  //
+  // We generate an AffineConstraints<> object
   template <int dim>
   void LaplaceProblem<dim>::assemble_multigrid()
   {
-    MappingQ1<dim>                      mapping;
-    MeshWorker::IntegrationInfoBox<dim> info_box;
-    UpdateFlags                         update_flags =
-      update_values | update_gradients | update_hessians;
-    info_box.add_update_flags_all(update_flags);
-    info_box.initialize(fe, mapping);
-
-    MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
-    MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double>> assembler;
-    assembler.initialize(mg_constrained_dofs);
-    assembler.initialize(mg_matrices);
-    assembler.initialize_interfaces(mg_interface_in, mg_interface_out);
-
-    LaplaceIntegrator<dim> matrix_integrator;
-    MeshWorker::integration_loop<dim, dim>(dof_handler.begin_mg(),
-                                           dof_handler.end_mg(),
-                                           dof_info,
-                                           info_box,
-                                           matrix_integrator,
-                                           assembler);
-
-    const unsigned int nlevels = triangulation.n_levels();
-    for (unsigned int level = 0; level < nlevels; ++level)
+    MappingQ1<dim>     mapping;
+    const unsigned int n_levels = triangulation.n_levels();
+
+    std::vector<AffineConstraints<>> boundary_constraints(n_levels);
+    for (unsigned int level = 0; level < n_levels; ++level)
       {
-        for (unsigned int i = 0; i < dof_handler.n_dofs(level); ++i)
-          if (mg_constrained_dofs.is_boundary_index(level, i) ||
-              mg_constrained_dofs.at_refinement_edge(level, i))
-            mg_matrices[level].set(i, i, 1.);
+        IndexSet dofset;
+        DoFTools::extract_locally_relevant_level_dofs(dof_handler,
+                                                      level,
+                                                      dofset);
+        boundary_constraints[level].reinit(dofset);
+        boundary_constraints[level].add_lines(
+          mg_constrained_dofs.get_refinement_edge_indices(level));
+        boundary_constraints[level].add_lines(
+          mg_constrained_dofs.get_boundary_indices(level));
+        boundary_constraints[level].close();
       }
+
+    typedef decltype(dof_handler.begin_mg()) Iterator;
+
+    auto cell_worker = [&](const Iterator &  cell,
+                           ScratchData<dim> &scratch_data,
+                           CopyData &        copy_data) {
+      this->cell_worker(cell, scratch_data, copy_data);
+    };
+
+    auto copier = [&](const CopyData &cd) {
+      boundary_constraints[cd.level].distribute_local_to_global(
+        cd.cell_matrix, cd.local_dof_indices, mg_matrices[cd.level]);
+
+      const unsigned int dofs_per_cell = cd.local_dof_indices.size();
+
+      // TODO EXPLAIN:
+
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+          if (mg_constrained_dofs.is_interface_matrix_entry(
+                cd.level, cd.local_dof_indices[i], cd.local_dof_indices[j]))
+            {
+              mg_interface_matrices[cd.level].add(cd.local_dof_indices[i],
+                                                  cd.local_dof_indices[j],
+                                                  cd.cell_matrix(i, j));
+            }
+    };
+
+    const unsigned int n_gauss_points = degree + 1;
+
+    ScratchData<dim> scratch_data(mapping,
+                                  fe,
+                                  n_gauss_points,
+                                  update_values | update_gradients |
+                                    update_JxW_values |
+                                    update_quadrature_points);
+
+    MeshWorker::mesh_loop(dof_handler.begin_mg(),
+                          dof_handler.end_mg(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          CopyData(),
+                          MeshWorker::assemble_own_cells);
   }
 
 
@@ -521,8 +567,8 @@ namespace Step16
     // initialize both up and down versions of the operator with the matrices
     // we already built:
     mg::Matrix<Vector<double>> mg_matrix(mg_matrices);
-    mg::Matrix<Vector<double>> mg_interface_up(mg_interface_in);
-    mg::Matrix<Vector<double>> mg_interface_down(mg_interface_out);
+    mg::Matrix<Vector<double>> mg_interface_up(mg_interface_matrices);
+    mg::Matrix<Vector<double>> mg_interface_down(mg_interface_matrices);
 
     // Now, we are ready to set up the V-cycle operator and the multilevel
     // preconditioner.
@@ -541,6 +587,9 @@ namespace Step16
     solution = 0;
 
     solver.solve(system_matrix, solution, system_rhs, preconditioner);
+    std::cout << "  Number of CG iterations: " << solver_control.last_step()
+              << "\n"
+              << std::endl;
     constraints.distribute(solution);
   }
 
@@ -551,10 +600,7 @@ namespace Step16
   // The following two functions postprocess a solution once it is
   // computed. In particular, the first one refines the mesh at the beginning
   // of each cycle while the second one outputs results at the end of each
-  // such cycle. The functions are almost unchanged from those in step-6, with
-  // the exception of one minor difference: we generate output in VTK
-  // format, to use the more modern visualization programs available today
-  // compared to those that were available when step-6 was written.
+  // such cycle. The functions are almost unchanged from those in step-6.
   template <int dim>
   void LaplaceProblem<dim>::refine_grid()
   {
@@ -562,7 +608,7 @@ namespace Step16
 
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
-      QGauss<dim - 1>(3),
+      QGauss<dim - 1>(degree + 2),
       std::map<types::boundary_id, const Function<dim> *>(),
       solution,
       estimated_error_per_cell);
@@ -600,18 +646,18 @@ namespace Step16
   {
     for (unsigned int cycle = 0; cycle < 8; ++cycle)
       {
-        deallog << "Cycle " << cycle << std::endl;
+        std::cout << "Cycle " << cycle << std::endl;
 
         if (cycle == 0)
           {
             GridGenerator::hyper_ball(triangulation);
-            triangulation.refine_global(1);
+            triangulation.refine_global(2);
           }
         else
           refine_grid();
 
-        deallog << "   Number of active cells:       "
-                << triangulation.n_active_cells() << std::endl;
+        std::cout << "   Number of active cells:       "
+                  << triangulation.n_active_cells() << std::endl;
 
         setup_system();
 
@@ -634,8 +680,6 @@ int main()
     {
       using namespace Step16;
 
-      deallog.depth_console(2);
-
       LaplaceProblem<2> laplace_problem(1);
       laplace_problem.run();
     }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.