From: tcclevenger <tcleven@clemson.edu>
Date: Mon, 9 Oct 2017 21:57:10 +0000 (-0400)
Subject: added is_interface_matrix_entry and step-50-mesh_loop
X-Git-Tag: v9.0.0-rc1~807^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5406%2Fhead;p=dealii.git

added is_interface_matrix_entry and step-50-mesh_loop
---

diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h
index a135930089..d319a859f1 100644
--- a/include/deal.II/multigrid/mg_constrained_dofs.h
+++ b/include/deal.II/multigrid/mg_constrained_dofs.h
@@ -110,6 +110,17 @@ public:
   bool at_refinement_edge (const unsigned int level,
                            const types::global_dof_index index) const;
 
+
+  /**
+   * Determine whether the (i,j) entry of the interface matrix
+   * on a given level should be set. This is taken in terms of
+   * dof i, that is, return true if i is at a refinement edge,
+   * j is not, and both are not on the external boundary.
+   */
+  bool is_interface_matrix_entry (const unsigned int level,
+                                  const types::global_dof_index i,
+                                  const types::global_dof_index j) const;
+
   /**
    * Return the indices of level dofs on the given level that are subject to
    * Dirichlet boundary conditions (as set by the @p function_map parameter in
@@ -312,6 +323,24 @@ MGConstrainedDoFs::at_refinement_edge (const unsigned int level,
   return refinement_edge_indices[level].is_element(index);
 }
 
+inline
+bool
+MGConstrainedDoFs::is_interface_matrix_entry (const unsigned int level,
+                                              const types::global_dof_index i,
+                                              const types::global_dof_index j) const
+{
+  const IndexSet &interface_dofs_on_level
+    = this->get_refinement_edge_indices(level);
+
+  return interface_dofs_on_level.is_element(i)   // at_refinement_edge(i)
+         &&
+         !interface_dofs_on_level.is_element(j)  // !at_refinement_edge(j)
+         &&
+         !this->is_boundary_index(level, i)      // !on_boundary(i)
+         &&
+         !this->is_boundary_index(level, j);     // !on_boundary(j)
+}
+
 
 
 
diff --git a/include/deal.II/multigrid/mg_tools.h b/include/deal.II/multigrid/mg_tools.h
index 0488040e85..fdcf167e2d 100644
--- a/include/deal.II/multigrid/mg_tools.h
+++ b/include/deal.II/multigrid/mg_tools.h
@@ -28,6 +28,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 template <int dim, int spacedim> class DoFHandler;
+class MGConstrainedDoFs;
 
 /* !@addtogroup mg */
 /* @{ */
@@ -137,6 +138,20 @@ namespace MGTools
                                    const unsigned int                 level,
                                    const Table<2,DoFTools::Coupling> &flux_mask);
 
+
+  /**
+   * Create sparsity pattern for interface_in/out matrices used in a multigrid computation.
+   * These matrices contain an entry representing the coupling of degrees of
+   * freedom on a refinement edge to those not on the refinement edge of a certain level.
+   */
+  template <typename DoFHandlerType, typename SparsityPatternType>
+  void
+  make_interface_sparsity_pattern (const DoFHandlerType    &dof_handler,
+                                   const MGConstrainedDoFs &mg_constrained_dofs,
+                                   SparsityPatternType     &sparsity,
+                                   const unsigned int      level);
+
+
   /**
    * Count the dofs block-wise on each level.
    *
diff --git a/source/multigrid/mg_tools.cc b/source/multigrid/mg_tools.cc
index e69cd338d3..b0ec9651fe 100644
--- a/source/multigrid/mg_tools.cc
+++ b/source/multigrid/mg_tools.cc
@@ -28,6 +28,7 @@
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/multigrid/mg_tools.h>
 #include <deal.II/multigrid/mg_base.h>
+#include <deal.II/multigrid/mg_constrained_dofs.h>
 #include <deal.II/base/mg_level_object.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.h>
@@ -982,6 +983,39 @@ namespace MGTools
 
 
 
+  template <typename DoFHandlerType, typename SparsityPatternType>
+  void
+  make_interface_sparsity_pattern (const DoFHandlerType    &dof,
+                                   const MGConstrainedDoFs &mg_constrained_dofs,
+                                   SparsityPatternType     &sparsity,
+                                   const unsigned int      level)
+  {
+    const types::global_dof_index n_dofs = dof.n_dofs(level);
+    (void)n_dofs;
+
+    Assert (sparsity.n_rows() == n_dofs,
+            ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
+    Assert (sparsity.n_cols() == n_dofs,
+            ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
+
+    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
+    std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
+    typename DoFHandlerType::cell_iterator cell = dof.begin(level),
+                                           endc = dof.end(level);
+    for (; cell!=endc; ++cell)
+      if (cell->is_locally_owned_on_level())
+        {
+          cell->get_mg_dof_indices (dofs_on_this_cell);
+          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(level,dofs_on_this_cell[i],dofs_on_this_cell[j]))
+                sparsity.add (dofs_on_this_cell[i],
+                              dofs_on_this_cell[j]);
+        }
+  }
+
+
+
   template <int dim, int spacedim>
   void
   count_dofs_per_component (const DoFHandler<dim,spacedim> &dof_handler,
diff --git a/source/multigrid/mg_tools.inst.in b/source/multigrid/mg_tools.inst.in
index e5afa6b37c..45e5c10cad 100644
--- a/source/multigrid/mg_tools.inst.in
+++ b/source/multigrid/mg_tools.inst.in
@@ -27,6 +27,13 @@ for (PATTERN : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS;
         const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
         PATTERN &,
         const unsigned int);
+
+    template void
+    make_interface_sparsity_pattern<DoFHandler<deal_II_dimension, deal_II_space_dimension>, PATTERN> (
+        const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+        const MGConstrainedDoFs &,
+        PATTERN &,
+        const unsigned int);
 #endif
 
 #if deal_II_dimension == deal_II_space_dimension
diff --git a/tests/meshworker/step-50-mesh_loop.cc b/tests/meshworker/step-50-mesh_loop.cc
new file mode 100644
index 0000000000..14ed3165c8
--- /dev/null
+++ b/tests/meshworker/step-50-mesh_loop.cc
@@ -0,0 +1,558 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2016 - 2017 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ *
+ */
+
+// Step-50 using mesh_loop()
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/conditional_ostream.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+
+#include <deal.II/lac/generic_linear_algebra.h>
+
+
+namespace LA
+{
+  using namespace dealii::LinearAlgebraTrilinos;
+}
+
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/base/work_stream.h>
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+namespace Step50
+{
+  using namespace dealii;
+
+
+  template <int dim>
+  struct ScratchData
+  {
+    ScratchData (const FiniteElement<dim> &fe,
+                 const unsigned int quadrature_degree)
+      :
+      fe_values (fe,
+                 QGauss<dim>(quadrature_degree),
+                 update_values   | update_gradients |
+                 update_quadrature_points | update_JxW_values)
+    {}
+
+    ScratchData (const ScratchData<dim> &scratch_data)
+      :
+      fe_values (scratch_data.fe_values.get_fe(),
+                 scratch_data.fe_values.get_quadrature(),
+                 update_values   | update_gradients |
+                 update_quadrature_points | update_JxW_values)
+    {}
+
+    FEValues<dim>     fe_values;
+
+  };
+  struct CopyData
+  {
+    unsigned int level;
+    unsigned int dofs_per_cell;
+
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+  };
+
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem (const unsigned int deg);
+    void run ();
+
+  private:
+    void setup_system ();
+
+    template <class IteratorType>
+    void assemble_cell (const IteratorType &cell,
+                        ScratchData<dim> &scratch_data,
+                        CopyData &copy_data);
+
+    void assemble_system_and_multigrid ();
+    void solve ();
+    void refine_grid ();
+
+    parallel::distributed::Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>    mg_dof_handler;
+
+    typedef LA::MPI::SparseMatrix matrix_t;
+    typedef LA::MPI::Vector vector_t;
+
+    matrix_t system_matrix;
+
+    IndexSet locally_relevant_set;
+
+    ConstraintMatrix     constraints;
+
+    vector_t       solution;
+    vector_t       system_rhs;
+
+    const unsigned int degree;
+
+    MGLevelObject<matrix_t> mg_matrices;
+    MGLevelObject<matrix_t> mg_interface_matrices;
+    MGConstrainedDoFs       mg_constrained_dofs;
+
+
+  };
+
+  template <int dim>
+  class Coefficient : public Function<dim>
+  {
+  public:
+    Coefficient () : Function<dim>() {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void value_list (const std::vector<Point<dim> > &points,
+                             std::vector<double>            &values,
+                             const unsigned int              component = 0) const;
+  };
+
+
+
+  template <int dim>
+  double Coefficient<dim>::value (const Point<dim> &p,
+                                  const unsigned int) const
+  {
+    if (p.square() < 0.5*0.5)
+      return 5;
+    else
+      return 1;
+  }
+
+
+
+  template <int dim>
+  void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+                                     std::vector<double>            &values,
+                                     const unsigned int              component) const
+  {
+    (void)component;
+    const unsigned int n_points = points.size();
+
+    Assert (values.size() == n_points,
+            ExcDimensionMismatch (values.size(), n_points));
+
+    Assert (component == 0,
+            ExcIndexRange (component, 0, 1));
+
+    for (unsigned int i=0; i<n_points; ++i)
+      values[i] = Coefficient<dim>::value (points[i]);
+  }
+
+
+
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
+    :
+    triangulation (MPI_COMM_WORLD,Triangulation<dim>::
+                   limit_level_difference_at_vertices,
+                   parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
+    fe (degree),
+    mg_dof_handler (triangulation),
+    degree(degree)
+  {}
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::setup_system ()
+  {
+    mg_dof_handler.distribute_dofs (fe);
+    mg_dof_handler.distribute_mg_dofs ();
+
+    DoFTools::extract_locally_relevant_dofs (mg_dof_handler,
+                                             locally_relevant_set);
+
+    solution.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
+    system_rhs.reinit(mg_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
+
+    constraints.reinit (locally_relevant_set);
+    DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
+
+    std::set<types::boundary_id>         dirichlet_boundary_ids;
+    typename FunctionMap<dim>::type      dirichlet_boundary;
+    Functions::ConstantFunction<dim>                    homogeneous_dirichlet_bc (1.0);
+    dirichlet_boundary_ids.insert(0);
+    dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
+    VectorTools::interpolate_boundary_values (mg_dof_handler,
+                                              dirichlet_boundary,
+                                              constraints);
+    constraints.close ();
+
+    DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(), mg_dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern (mg_dof_handler, dsp, constraints);
+    system_matrix.reinit (mg_dof_handler.locally_owned_dofs(), dsp, MPI_COMM_WORLD, true);
+
+
+    mg_constrained_dofs.clear();
+    mg_constrained_dofs.initialize(mg_dof_handler);
+    mg_constrained_dofs.make_zero_boundary_constraints(mg_dof_handler, dirichlet_boundary_ids);
+
+
+    const unsigned int n_levels = triangulation.n_global_levels();
+
+    mg_interface_matrices.resize(0, n_levels-1);
+    mg_interface_matrices.clear_elements ();
+    mg_matrices.resize(0, n_levels-1);
+    mg_matrices.clear_elements ();
+
+    for (unsigned int level=0; level<n_levels; ++level)
+      {
+        {
+          DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(level),
+                                     mg_dof_handler.n_dofs(level));
+          MGTools::make_sparsity_pattern(mg_dof_handler, dsp, level);
+
+          mg_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
+                                    mg_dof_handler.locally_owned_mg_dofs(level),
+                                    dsp,
+                                    MPI_COMM_WORLD, true);
+        }
+
+        {
+          DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(level),
+                                     mg_dof_handler.n_dofs(level));
+          MGTools::make_interface_sparsity_pattern(mg_dof_handler, mg_constrained_dofs, dsp, level);
+
+          mg_interface_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
+                                              mg_dof_handler.locally_owned_mg_dofs(level),
+                                              dsp,
+                                              MPI_COMM_WORLD, true);
+        }
+      }
+  }
+
+
+
+  template <int dim>
+  template <class IteratorType>
+  void LaplaceProblem<dim>::assemble_cell(const IteratorType &cell,
+                                          ScratchData<dim> &scratch_data,
+                                          CopyData &copy_data)
+  {
+    const unsigned int level = cell->level();
+    copy_data.level = level;
+
+    const unsigned int dofs_per_cell = scratch_data.fe_values.get_fe().dofs_per_cell;
+    copy_data.dofs_per_cell = dofs_per_cell;
+    const unsigned int n_q_points      = scratch_data.fe_values.get_quadrature().size();
+
+    copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+    if (!cell->is_level_cell())
+      copy_data.cell_rhs.reinit (dofs_per_cell);
+
+    copy_data.local_dof_indices.resize(dofs_per_cell);
+    cell->get_active_or_mg_dof_indices (copy_data.local_dof_indices);
+
+    scratch_data.fe_values.reinit (cell);
+
+    const Coefficient<dim> coefficient;
+    std::vector<double>    coefficient_values (n_q_points);
+    coefficient.value_list (scratch_data.fe_values.get_quadrature_points(),
+                            coefficient_values);
+
+    for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+      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_values[q_point] *
+                                           scratch_data.fe_values.shape_grad (i, q_point) *
+                                           scratch_data.fe_values.shape_grad (j, q_point) *
+                                           scratch_data.fe_values.JxW (q_point));
+          if (!cell->is_level_cell())
+            copy_data.cell_rhs(i) += (scratch_data.fe_values.shape_value(i,q_point) *
+                                      10.0 *
+                                      scratch_data.fe_values.JxW (q_point));
+        }
+  }
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::assemble_system_and_multigrid ()
+  {
+    std::vector<ConstraintMatrix> boundary_constraints (triangulation.n_global_levels());
+    for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
+      {
+        IndexSet dofset;
+        DoFTools::extract_locally_relevant_level_dofs (mg_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 ();
+      }
+
+    auto cell_worker_active = [&] (const decltype(mg_dof_handler.begin_active()) & cell,
+                                   ScratchData<dim> &scratch_data,
+                                   CopyData &copy_data)
+    {
+      this->assemble_cell(cell, scratch_data, copy_data);
+    };
+
+    auto cell_worker_mg = [&] (const decltype(mg_dof_handler.begin_mg()) & cell,
+                               ScratchData<dim> &scratch_data,
+                               CopyData &copy_data)
+    {
+      this->assemble_cell(cell, scratch_data, copy_data);
+    };
+
+    auto copier_active = [&](const CopyData &c)
+    {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             system_matrix,
+                                             system_rhs);
+    };
+
+    auto copier_mg = [&](const CopyData &c)
+    {
+      boundary_constraints[c.level]
+      .distribute_local_to_global (c.cell_matrix,
+                                   c.local_dof_indices,
+                                   mg_matrices[c.level]);
+
+      for (unsigned int i=0; i<c.dofs_per_cell; ++i)
+        for (unsigned int j=0; j<c.dofs_per_cell; ++j)
+          if (mg_constrained_dofs.is_interface_matrix_entry(c.level, c.local_dof_indices[i], c.local_dof_indices[j]))
+            mg_interface_matrices[c.level].add(c.local_dof_indices[i],c.local_dof_indices[j],c.cell_matrix(i,j));
+    };
+
+    MeshWorker::mesh_loop(mg_dof_handler.begin_active(),
+                          mg_dof_handler.end(),
+                          cell_worker_active,
+                          copier_active,
+                          ScratchData<dim>(fe, degree+1),
+                          CopyData(),
+                          MeshWorker::assemble_own_cells);
+
+    MeshWorker::mesh_loop(mg_dof_handler.begin_mg(),
+                          mg_dof_handler.end_mg(),
+                          cell_worker_mg,
+                          copier_mg,
+                          ScratchData<dim>(fe, degree+1),
+                          CopyData(),
+                          MeshWorker::assemble_own_cells);
+
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+
+    for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
+      {
+        mg_matrices[level].compress(VectorOperation::add);
+        deallog << "mg_matrices[" << level << "]: " << mg_matrices[level].frobenius_norm() << std::endl;
+        mg_interface_matrices[level].compress(VectorOperation::add);
+        deallog << "mg_interface_matrices[" << level << "]: " << mg_interface_matrices[level].frobenius_norm() << std::endl;
+      }
+  }
+
+
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::solve ()
+  {
+    MGTransferPrebuilt<vector_t> mg_transfer(mg_constrained_dofs);
+    mg_transfer.build_matrices(mg_dof_handler);
+
+    matrix_t &coarse_matrix = mg_matrices[0];
+
+    SolverControl coarse_solver_control (1000, 1e-10, false, false);
+    SolverCG<vector_t> coarse_solver(coarse_solver_control);
+    PreconditionIdentity id;
+    MGCoarseGridIterativeSolver<vector_t, SolverCG<vector_t>, matrix_t, PreconditionIdentity>
+    coarse_grid_solver(coarse_solver, coarse_matrix, id);
+
+    typedef LA::MPI::PreconditionJacobi Smoother;
+    MGSmootherPrecondition<matrix_t, Smoother, vector_t> mg_smoother;
+    mg_smoother.initialize(mg_matrices, Smoother::AdditionalData(0.5));
+    mg_smoother.set_steps(2);
+
+    mg::Matrix<vector_t> mg_matrix(mg_matrices);
+    mg::Matrix<vector_t> mg_interface_up(mg_interface_matrices);
+    mg::Matrix<vector_t> mg_interface_down(mg_interface_matrices);
+
+    Multigrid<vector_t > mg(mg_matrix,
+                            coarse_grid_solver,
+                            mg_transfer,
+                            mg_smoother,
+                            mg_smoother);
+    mg.set_edge_matrices(mg_interface_down, mg_interface_up);
+
+    PreconditionMG<dim, vector_t, MGTransferPrebuilt<vector_t> >
+    preconditioner(mg_dof_handler, mg, mg_transfer);
+
+
+    SolverControl solver_control (500, 1e-8*system_rhs.l2_norm(), false);
+    SolverCG<vector_t> solver (solver_control);
+
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+
+    deallog << "   CG converged in " << solver_control.last_step() << " iterations." << std::endl;
+
+    constraints.distribute (solution);
+  }
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::refine_grid ()
+  {
+    for (typename Triangulation<dim>::active_cell_iterator cell=triangulation.begin_active(); cell != triangulation.end(); ++cell)
+      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+        if (cell->vertex(v)[0] <= 0.5 && cell->vertex(v)[1] <= 0.5)
+          cell->set_refine_flag();
+    triangulation.execute_coarsening_and_refinement ();
+  }
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::run ()
+  {
+    for (unsigned int cycle=0; cycle<3; ++cycle)
+      {
+        deallog << "Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube (triangulation);
+            triangulation.refine_global (4);
+          }
+        else
+          refine_grid ();
+
+        deallog << "   Number of active cells:       "
+                << triangulation.n_global_active_cells()
+                << std::endl;
+
+        setup_system ();
+
+        deallog << "   Number of degrees of freedom: "
+                << mg_dof_handler.n_dofs()
+                << " (by level: ";
+        for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
+          deallog << mg_dof_handler.n_dofs(level)
+                  << (level == triangulation.n_global_levels()-1
+                      ? ")" : ", ");
+        deallog << std::endl;
+
+        assemble_system_and_multigrid ();
+
+        solve ();
+      }
+  }
+}
+
+
+int main (int argc, char *argv[])
+{
+  dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  mpi_initlog(true);
+
+  try
+    {
+      using namespace dealii;
+      using namespace Step50;
+
+      LaplaceProblem<2> laplace_problem(1/*degree*/);
+      laplace_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      throw;
+    }
+
+  return 0;
+}
diff --git a/tests/meshworker/step-50-mesh_loop.with_trilinos=true.with_p4est=true.mpirun=3.output b/tests/meshworker/step-50-mesh_loop.with_trilinos=true.with_p4est=true.mpirun=3.output
new file mode 100644
index 0000000000..7be1e95d04
--- /dev/null
+++ b/tests/meshworker/step-50-mesh_loop.with_trilinos=true.with_p4est=true.mpirun=3.output
@@ -0,0 +1,55 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       256
+DEAL::   Number of degrees of freedom: 289 (by level: 4, 9, 25, 81, 289)
+DEAL::mg_matrices[0]: 2.78887
+DEAL::mg_interface_matrices[0]: 0.00000
+DEAL::mg_matrices[1]: 7.37274
+DEAL::mg_interface_matrices[1]: 0.00000
+DEAL::mg_matrices[2]: 20.4478
+DEAL::mg_interface_matrices[2]: 0.00000
+DEAL::mg_matrices[3]: 47.8346
+DEAL::mg_interface_matrices[3]: 0.00000
+DEAL::mg_matrices[4]: 101.551
+DEAL::mg_interface_matrices[4]: 0.00000
+DEAL:cg::Starting value 21.0544
+DEAL:cg::Convergence step 7 value 3.37709e-08
+DEAL::   CG converged in 7 iterations.
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       499
+DEAL::   Number of degrees of freedom: 550 (by level: 4, 9, 25, 81, 289, 361)
+DEAL::mg_matrices[0]: 2.78887
+DEAL::mg_interface_matrices[0]: 0.00000
+DEAL::mg_matrices[1]: 7.37274
+DEAL::mg_interface_matrices[1]: 0.00000
+DEAL::mg_matrices[2]: 20.4478
+DEAL::mg_interface_matrices[2]: 0.00000
+DEAL::mg_matrices[3]: 47.8346
+DEAL::mg_interface_matrices[3]: 0.00000
+DEAL::mg_matrices[4]: 101.551
+DEAL::mg_interface_matrices[4]: 0.00000
+DEAL::mg_matrices[5]: 195.872
+DEAL::mg_interface_matrices[5]: 3.31662
+DEAL:cg::Starting value 29.0111
+DEAL:cg::Convergence step 7 value 9.42221e-08
+DEAL::   CG converged in 7 iterations.
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       1366
+DEAL::   Number of degrees of freedom: 1451 (by level: 4, 9, 25, 81, 289, 361, 1225)
+DEAL::mg_matrices[0]: 2.78887
+DEAL::mg_interface_matrices[0]: 0.00000
+DEAL::mg_matrices[1]: 7.37274
+DEAL::mg_interface_matrices[1]: 0.00000
+DEAL::mg_matrices[2]: 20.4478
+DEAL::mg_interface_matrices[2]: 0.00000
+DEAL::mg_matrices[3]: 47.8346
+DEAL::mg_interface_matrices[3]: 0.00000
+DEAL::mg_matrices[4]: 101.551
+DEAL::mg_interface_matrices[4]: 0.00000
+DEAL::mg_matrices[5]: 195.872
+DEAL::mg_interface_matrices[5]: 3.31662
+DEAL::mg_matrices[6]: 397.385
+DEAL::mg_interface_matrices[6]: 4.65475
+DEAL:cg::Starting value 40.5102
+DEAL:cg::Convergence step 7 value 2.64596e-07
+DEAL::   CG converged in 7 iterations.
diff --git a/tests/multigrid/interface_matrix_entry_01.cc b/tests/multigrid/interface_matrix_entry_01.cc
new file mode 100644
index 0000000000..89101a0c47
--- /dev/null
+++ b/tests/multigrid/interface_matrix_entry_01.cc
@@ -0,0 +1,130 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// check functions MGTools::make_interface_sparsity_pattern()
+// and MGConstrainedDoFs::is_interface_matrix_entry()
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+#include <deal.II/multigrid/multigrid.h>
+
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/generic_linear_algebra.h>
+
+using namespace dealii;
+
+template <int dim>
+void test ()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD,Triangulation<dim>::
+                                                 limit_level_difference_at_vertices,
+                                                 parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::hyper_cube(tria,0,1);
+  tria.refine_global(1);
+  for (typename parallel::distributed::Triangulation<dim>::active_cell_iterator cell = tria.begin_active();
+       cell != tria.end(); ++cell)
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+      {
+        if (dim == 2)
+          if (cell->vertex(v)[0] < 0.5 && cell->vertex(v)[1] < 0.5)
+            {
+              cell->set_refine_flag();
+              break;
+            }
+        if (dim == 3)
+          if (cell->vertex(v)[0] < 0.5 && cell->vertex(v)[1] < 0.5 && cell->vertex(v)[2] < 0.5)
+            {
+              cell->set_refine_flag();
+              break;
+            }
+      }
+  tria.execute_coarsening_and_refinement();
+
+  FE_Q<dim>            fe(2);
+  DoFHandler<dim>    mg_dof_handler(tria);
+  IndexSet locally_relevant_set;
+  ConstraintMatrix     constraints;
+  MGConstrainedDoFs    mg_constrained_dofs;
+
+  mg_dof_handler.distribute_dofs (fe);
+  mg_dof_handler.distribute_mg_dofs ();
+
+  DoFTools::extract_locally_relevant_dofs (mg_dof_handler,
+                                           locally_relevant_set);
+
+  constraints.reinit (locally_relevant_set);
+  DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints);
+
+  std::set<types::boundary_id>         dirichlet_boundary_ids;
+  typename FunctionMap<dim>::type      dirichlet_boundary;
+  Functions::ConstantFunction<dim>                    homogeneous_dirichlet_bc (0.0);
+  dirichlet_boundary_ids.insert(0);
+  dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
+  VectorTools::interpolate_boundary_values (mg_dof_handler,
+                                            dirichlet_boundary,
+                                            constraints);
+  constraints.close ();
+
+  mg_constrained_dofs.clear();
+  mg_constrained_dofs.initialize(mg_dof_handler);
+  mg_constrained_dofs.make_zero_boundary_constraints(mg_dof_handler, dirichlet_boundary_ids);
+
+  for (unsigned int level=0; level<tria.n_levels(); ++level)
+    {
+      DynamicSparsityPattern dsp(mg_dof_handler.n_dofs(level),
+                                 mg_dof_handler.n_dofs(level));
+      MGTools::make_interface_sparsity_pattern(mg_dof_handler, mg_constrained_dofs, dsp, level);
+      dsp.print(deallog.get_file_stream());
+    }
+}
+
+
+
+int main (int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll all;
+
+  deallog.push("2d");
+  test<2>();
+  deallog.pop();
+  deallog.push("3d");
+  test<3>();
+  deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/multigrid/interface_matrix_entry_01.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/multigrid/interface_matrix_entry_01.with_mpi=true.with_p4est=true.mpirun=1.output
new file mode 100644
index 0000000000..7be9b7be7b
--- /dev/null
+++ b/tests/multigrid/interface_matrix_entry_01.with_mpi=true.with_p4est=true.mpirun=1.output
@@ -0,0 +1,338 @@
+
+[0]
+[1]
+[2]
+[3]
+[4]
+[5]
+[6]
+[7]
+[8]
+[0]
+[1]
+[2]
+[3]
+[4]
+[5]
+[6]
+[7]
+[8]
+[9]
+[10]
+[11]
+[12]
+[13]
+[14]
+[15]
+[16]
+[17]
+[18]
+[19]
+[20]
+[21]
+[22]
+[23]
+[24]
+[0]
+[1]
+[2]
+[3]
+[4]
+[5]
+[6]
+[7]
+[8]
+[9]
+[10,3,5,13,14,18,24]
+[11,3,5,13,14]
+[12]
+[13]
+[14]
+[15]
+[16,3,7,13,18,20,24]
+[17]
+[18]
+[19,3,7,18,20]
+[20]
+[21,3,13,18,24]
+[22,3,13,18,24]
+[23,3,13,18,24]
+[24]
+[0]
+[1]
+[2]
+[3]
+[4]
+[5]
+[6]
+[7]
+[8]
+[9]
+[10]
+[11]
+[12]
+[13]
+[14]
+[15]
+[16]
+[17]
+[18]
+[19]
+[20]
+[21]
+[22]
+[23]
+[24]
+[25]
+[26]
+[0]
+[1]
+[2]
+[3]
+[4]
+[5]
+[6]
+[7]
+[8]
+[9]
+[10]
+[11]
+[12]
+[13]
+[14]
+[15]
+[16]
+[17]
+[18]
+[19]
+[20]
+[21]
+[22]
+[23]
+[24]
+[25]
+[26]
+[27]
+[28]
+[29]
+[30]
+[31]
+[32]
+[33]
+[34]
+[35]
+[36]
+[37]
+[38]
+[39]
+[40]
+[41]
+[42]
+[43]
+[44]
+[45]
+[46]
+[47]
+[48]
+[49]
+[50]
+[51]
+[52]
+[53]
+[54]
+[55]
+[56]
+[57]
+[58]
+[59]
+[60]
+[61]
+[62]
+[63]
+[64]
+[65]
+[66]
+[67]
+[68]
+[69]
+[70]
+[71]
+[72]
+[73]
+[74]
+[75]
+[76]
+[77]
+[78]
+[79]
+[80]
+[81]
+[82]
+[83]
+[84]
+[85]
+[86]
+[87]
+[88]
+[89]
+[90]
+[91]
+[92]
+[93]
+[94]
+[95]
+[96]
+[97]
+[98]
+[99]
+[100]
+[101]
+[102]
+[103]
+[104]
+[105]
+[106]
+[107]
+[108]
+[109]
+[110]
+[111]
+[112]
+[113]
+[114]
+[115]
+[116]
+[117]
+[118]
+[119]
+[120]
+[121]
+[122]
+[123]
+[124]
+[0]
+[1]
+[2]
+[3]
+[4]
+[5]
+[6]
+[7]
+[8]
+[9]
+[10]
+[11]
+[12]
+[13]
+[14]
+[15]
+[16]
+[17]
+[18]
+[19]
+[20]
+[21]
+[22]
+[23]
+[24]
+[25]
+[26]
+[27]
+[28]
+[29]
+[30,7,13,19,21,36,41,43,44,53,58,73,74,86,88,102,104,113,124]
+[31]
+[32]
+[33]
+[34,7,13,19,21,36,41,43,44,86,88,102,104]
+[35]
+[36]
+[37]
+[38,7,13,19,21,36,41,43,44,53,58,73,74]
+[39,7,13,19,21,36,41,43,44]
+[40]
+[41]
+[42]
+[43]
+[44]
+[45]
+[46]
+[47]
+[48,7,15,19,23,36,41,53,58,61,62,73,74,86,90,102,113,116,124]
+[49]
+[50]
+[51]
+[52]
+[53]
+[54,7,15,19,23,53,58,61,62,86,90,113,116]
+[55]
+[56,7,15,19,23,36,41,53,58,61,62,73,74]
+[57]
+[58]
+[59,7,15,19,23,53,58,61,62]
+[60]
+[61]
+[62]
+[63]
+[64,7,19,36,41,53,58,73,74,86,102,113,124]
+[65]
+[66]
+[67,7,19,36,41,53,58,73,74,86,102,113,124]
+[68,7,19,36,41,53,58,73,74,86,102,113,124]
+[69,7,19,36,41,53,58,73,74]
+[70,7,19,36,41,53,58,73,74]
+[71,7,19,36,41,53,58,73,74]
+[72]
+[73]
+[74]
+[75]
+[76]
+[77]
+[78,7,13,15,25,36,43,53,61,73,86,88,90,92,102,104,113,116,124]
+[79]
+[80,7,13,15,25,36,43,86,88,90,92,102,104]
+[81]
+[82,7,13,15,25,53,61,86,88,90,92,113,116]
+[83]
+[84]
+[85]
+[86]
+[87]
+[88]
+[89]
+[90]
+[91,7,13,15,25,86,88,90,92]
+[92]
+[93]
+[94,7,13,36,43,53,73,86,88,102,104,113,124]
+[95,7,13,36,43,86,88,102,104]
+[96]
+[97,7,13,36,43,53,73,86,88,102,104,113,124]
+[98]
+[99,7,13,36,43,53,73,86,88,102,104,113,124]
+[100,7,13,36,43,86,88,102,104]
+[101]
+[102]
+[103,7,13,36,43,86,88,102,104]
+[104]
+[105]
+[106,7,15,36,53,61,73,86,90,102,113,116,124]
+[107]
+[108,7,15,36,53,61,73,86,90,102,113,116,124]
+[109,7,15,53,61,86,90,113,116]
+[110]
+[111,7,15,36,53,61,73,86,90,102,113,116,124]
+[112]
+[113]
+[114,7,15,53,61,86,90,113,116]
+[115,7,15,53,61,86,90,113,116]
+[116]
+[117,7,36,53,73,86,102,113,124]
+[118,7,36,53,73,86,102,113,124]
+[119,7,36,53,73,86,102,113,124]
+[120,7,36,53,73,86,102,113,124]
+[121,7,36,53,73,86,102,113,124]
+[122,7,36,53,73,86,102,113,124]
+[123,7,36,53,73,86,102,113,124]
+[124]
+DEAL:0::OK