From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Wed, 4 Dec 2013 17:21:32 +0000 (+0000)
Subject: Two more test cases. Use common naming scheme for parallel assembly tests
X-Git-Tag: v8.1.0~97
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=016f23c3cd2e2416d70d2886ae3ad9017e7bdfcb;p=dealii.git

Two more test cases. Use common naming scheme for parallel assembly tests

git-svn-id: https://svn.dealii.org/trunk@31882 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/tests/deal.II/assemble_matrix_parallel.cc b/tests/deal.II/assemble_matrix_parallel_01.cc
similarity index 100%
rename from tests/deal.II/assemble_matrix_parallel.cc
rename to tests/deal.II/assemble_matrix_parallel_01.cc
diff --git a/tests/deal.II/assemble_block_matrix_parallel.output b/tests/deal.II/assemble_matrix_parallel_01.output
similarity index 100%
rename from tests/deal.II/assemble_block_matrix_parallel.output
rename to tests/deal.II/assemble_matrix_parallel_01.output
diff --git a/tests/deal.II/assemble_matrix_parallel_02.cc b/tests/deal.II/assemble_matrix_parallel_02.cc
new file mode 100644
index 0000000000..a2cd47d8dc
--- /dev/null
+++ b/tests/deal.II/assemble_matrix_parallel_02.cc
@@ -0,0 +1,520 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests that assembly in parallel works properly even when using two separate
+// loops on the same constraint matrix, otherwise similar to
+// assemble_matrix_parallel_01
+
+#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/work_stream.h>
+#include <deal.II/base/graph_coloring.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/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+  namespace Scratch
+  {
+    template <int dim>
+    struct Data
+    {
+      Data (const hp::FECollection<dim> &fe,
+            const hp::QCollection<dim>  &quadrature)
+        :
+        hp_fe_values(fe,
+                     quadrature,
+                     update_values    |  update_gradients |
+                     update_quadrature_points  |  update_JxW_values)
+      {}
+
+      Data (const Data &data)
+        :
+        hp_fe_values(data.hp_fe_values.get_mapping_collection(),
+                     data.hp_fe_values.get_fe_collection(),
+                     data.hp_fe_values.get_quadrature_collection(),
+                     data.hp_fe_values.get_update_flags())
+      {}
+
+      hp::FEValues<dim>               hp_fe_values;
+    };
+  }
+
+  namespace Copy
+  {
+    struct Data
+    {
+      Data (bool second_test = false)
+        :
+        second_test (second_test)
+      {}
+
+      Data (const Data &data)
+        :
+        second_test (data.second_test)
+      {}
+
+      std::vector<types::global_dof_index> local_dof_indices;
+      FullMatrix<double> local_matrix;
+      Vector<double> local_rhs;
+      const bool second_test;
+    };
+  }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem ();
+  ~LaplaceProblem ();
+
+  void run ();
+
+private:
+  void setup_system ();
+  void test_equality ();
+  void assemble_reference ();
+  void assemble_test ();
+  void assemble_test_1();
+  void assemble_test_2();
+  void solve ();
+  void create_coarse_grid ();
+  void postprocess ();
+
+  void local_assemble (const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                       Assembly::Scratch::Data<dim>  &scratch,
+                       Assembly::Copy::Data          &data);
+  void copy_local_to_global (const Assembly::Copy::Data &data);
+
+  std::vector<types::global_dof_index>
+  get_conflict_indices (typename hp::DoFHandler<dim>::active_cell_iterator const &cell) const;
+
+  Triangulation<dim>   triangulation;
+
+  hp::DoFHandler<dim>      dof_handler;
+  hp::FECollection<dim>    fe_collection;
+  hp::QCollection<dim>     quadrature_collection;
+  hp::QCollection<dim-1>   face_quadrature_collection;
+
+  ConstraintMatrix     constraints;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> reference_matrix;
+  SparseMatrix<double> test_matrix;
+  SparseMatrix<double> test_matrix_2;
+
+  Vector<double>       solution;
+  Vector<double>       reference_rhs;
+  Vector<double>       test_rhs;
+  Vector<double>       test_rhs_2;
+
+  std::vector<std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> > graph;
+
+  const unsigned int max_degree;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                            const unsigned int  /*component*/) const
+{
+  double sum = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    sum += std::sin(deal_II_numbers::PI*p[d]);
+  return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                           const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+  :
+  dof_handler (triangulation),
+  max_degree (5)
+{
+  if (dim == 2)
+    for (unsigned int degree=2; degree<=max_degree; ++degree)
+      {
+        fe_collection.push_back (FE_Q<dim>(degree));
+        quadrature_collection.push_back (QGauss<dim>(degree+1));
+        face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
+  else
+    for (unsigned int degree=1; degree<max_degree-1; ++degree)
+      {
+        fe_collection.push_back (FE_Q<dim>(degree));
+        quadrature_collection.push_back (QGauss<dim>(degree+1));
+        face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (typename hp::DoFHandler<dim>::active_cell_iterator const &cell) const
+{
+  std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+  cell->get_dof_indices(local_dof_indices);
+
+  constraints.resolve_indices(local_dof_indices);
+  return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  reference_matrix.clear();
+  test_matrix.clear();
+  dof_handler.distribute_dofs (fe_collection);
+
+  solution.reinit (dof_handler.n_dofs());
+  reference_rhs.reinit (dof_handler.n_dofs());
+  test_rhs.reinit (dof_handler.n_dofs());
+  test_rhs_2.reinit (dof_handler.n_dofs());
+
+  constraints.clear ();
+
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+  // add boundary conditions as inhomogeneous constraints here, do it after
+  // having added the hanging node constraints in order to be consistent and
+  // skip dofs that are already constrained (i.e., are hanging nodes on the
+  // boundary in 3D). In contrast to step-27, we choose a sine function.
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close ();
+
+  graph = GraphColoring::make_graph_coloring(dof_handler.begin_active(),dof_handler.end(),
+                                             static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+                                             (typename hp::DoFHandler<dim>::active_cell_iterator const &)> >
+                                             (std_cxx1x::bind(&LaplaceProblem<dim>::get_conflict_indices, this,std_cxx1x::_1)));
+
+
+  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+                                       dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                   constraints, false);
+  sparsity_pattern.copy_from (csp);
+
+  reference_matrix.reinit (sparsity_pattern);
+  test_matrix.reinit (sparsity_pattern);
+  test_matrix_2.reinit (sparsity_pattern);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                                     Assembly::Scratch::Data<dim>  &scratch,
+                                     Assembly::Copy::Data          &data)
+{
+  const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+  data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+  data.local_matrix = 0;
+
+  data.local_rhs.reinit (dofs_per_cell);
+  data.local_rhs = 0;
+
+  scratch.hp_fe_values.reinit (cell);
+
+  const FEValues<dim> &fe_values = scratch.hp_fe_values.get_present_fe_values ();
+
+  const RightHandSide<dim> rhs_function;
+
+  for (unsigned int q_point=0;
+       q_point<fe_values.n_quadrature_points;
+       ++q_point)
+    {
+      const double scale_mat = data.second_test ? numbers::PI : 1.;
+      const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            data.local_matrix(i,j) += (scale_mat *
+                                       fe_values.shape_grad(i,q_point) *
+                                       fe_values.shape_grad(j,q_point) *
+                                       fe_values.JxW(q_point));
+
+            data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+                                  rhs_value * scale_mat *
+                                  fe_values.JxW(q_point));
+        }
+    }
+
+  data.local_dof_indices.resize (dofs_per_cell);
+  cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::copy_local_to_global (const Assembly::Copy::Data &data)
+{
+  if (data.second_test)
+    constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+                                           data.local_dof_indices,
+                                           test_matrix_2, test_rhs_2);
+  else
+    constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+                                           data.local_dof_indices,
+                                           test_matrix, test_rhs);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  Assembly::Copy::Data copy_data;
+  Assembly::Scratch::Data<dim> assembly_data(fe_collection, quadrature_collection);
+
+  for (unsigned int color=0; color<graph.size(); ++color)
+    for (typename std::vector<typename hp::DoFHandler<dim>::active_cell_iterator>::const_iterator p = graph[color].begin();
+         p != graph[color].end(); ++p)
+      {
+        local_assemble(*p, assembly_data, copy_data);
+        copy_local_to_global(copy_data);
+      }
+
+  reference_matrix.add(1., test_matrix);
+  reference_rhs = test_rhs;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test_1 ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe_collection, quadrature_collection),
+         Assembly::Copy::Data (),
+         1);
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test_2 ()
+{
+  test_matrix_2 = 0;
+  test_rhs_2 = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe_collection, quadrature_collection),
+         Assembly::Copy::Data (true),
+         1);
+}
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test()
+{
+  // start two tasks that each run an assembly
+  Threads::TaskGroup<void> tasks;
+  tasks += Threads::new_task (&LaplaceProblem<dim>::assemble_test_1,
+                              *this);
+  tasks += Threads::new_task (&LaplaceProblem<dim>::assemble_test_2,
+                              *this);
+  tasks.join_all();
+
+  test_matrix.add(-1, reference_matrix);
+  const double frobenius_norm = test_matrix.frobenius_norm();
+
+  // the data should add up exactly (unfortunately, there is some roundoff due
+  // to the cell similarity detection, but there should not be any similarity
+  // for the hypershell geometry)
+  deallog << "log error in matrix norm: " << std::log(frobenius_norm) << std::endl;
+  test_rhs.add(-1., reference_rhs);
+  deallog << "log error in vector norm: " << std::log(test_rhs.l2_norm()) << std::endl;
+
+  test_matrix_2.add(-numbers::PI, reference_matrix);
+  deallog << "absolute error matrix assembly 2: " << test_matrix_2.frobenius_norm() << std::endl;
+  test_rhs_2.add(-numbers::PI, reference_rhs);
+  deallog << "absolute error vector assembly 2: " << test_rhs_2.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+    estimated_error_per_cell(i) = i;
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                   estimated_error_per_cell,
+                                                   0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+
+  for (typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
+       cell != dof_handler.end(); ++cell)
+    cell->set_active_fe_index (rand() % fe_collection.size());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_shell(triangulation,
+                                     Point<dim>(),
+                                     0.5, 1., (dim==3) ? 96 : 12, false);
+          triangulation.refine_global(2);
+        }
+
+      setup_system ();
+
+      assemble_reference ();
+      assemble_test ();
+
+      if (cycle < 2)
+        postprocess ();
+    }
+}
+
+
+
+int main ()
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-8);
+
+  {
+    deallog.push("2d");
+    LaplaceProblem<2> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/deal.II/assemble_matrix_parallel_02.output b/tests/deal.II/assemble_matrix_parallel_02.output
new file mode 100644
index 0000000000..627a39f770
--- /dev/null
+++ b/tests/deal.II/assemble_matrix_parallel_02.output
@@ -0,0 +1,13 @@
+
+DEAL:2d::log error in matrix norm: -inf
+DEAL:2d::log error in vector norm: -inf
+DEAL:2d::absolute error matrix assembly 2: 0
+DEAL:2d::absolute error vector assembly 2: 0
+DEAL:2d::log error in matrix norm: -inf
+DEAL:2d::log error in vector norm: -inf
+DEAL:2d::absolute error matrix assembly 2: 0
+DEAL:2d::absolute error vector assembly 2: 0
+DEAL:2d::log error in matrix norm: -inf
+DEAL:2d::log error in vector norm: -inf
+DEAL:2d::absolute error matrix assembly 2: 0
+DEAL:2d::absolute error vector assembly 2: 0
diff --git a/tests/deal.II/assemble_block_matrix_parallel.cc b/tests/deal.II/assemble_matrix_parallel_03.cc
similarity index 99%
rename from tests/deal.II/assemble_block_matrix_parallel.cc
rename to tests/deal.II/assemble_matrix_parallel_03.cc
index 2e20166c55..a0e1c59be0 100644
--- a/tests/deal.II/assemble_block_matrix_parallel.cc
+++ b/tests/deal.II/assemble_matrix_parallel_03.cc
@@ -16,7 +16,7 @@
 
 
 
-// same as assemble_matrix_parallel, but now for a BlockSparseMatrix
+// same as assemble_matrix_parallel_01, but now for a BlockSparseMatrix
 
 #include "../tests.h"
 
diff --git a/tests/deal.II/assemble_matrix_parallel.output b/tests/deal.II/assemble_matrix_parallel_03.output
similarity index 100%
rename from tests/deal.II/assemble_matrix_parallel.output
rename to tests/deal.II/assemble_matrix_parallel_03.output
diff --git a/tests/deal.II/assemble_matrix_parallel_04.cc b/tests/deal.II/assemble_matrix_parallel_04.cc
new file mode 100644
index 0000000000..148cdc137a
--- /dev/null
+++ b/tests/deal.II/assemble_matrix_parallel_04.cc
@@ -0,0 +1,525 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// same as assemble_matrix_parallel_02, but using a constraint that should
+// create conflicts on most cells (similar to mean value constraint)
+
+#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/work_stream.h>
+#include <deal.II/base/graph_coloring.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/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+  namespace Scratch
+  {
+    template <int dim>
+    struct Data
+    {
+      Data (const hp::FECollection<dim> &fe,
+            const hp::QCollection<dim>  &quadrature)
+        :
+        hp_fe_values(fe,
+                     quadrature,
+                     update_values    |  update_gradients |
+                     update_quadrature_points  |  update_JxW_values)
+      {}
+
+      Data (const Data &data)
+        :
+        hp_fe_values(data.hp_fe_values.get_mapping_collection(),
+                     data.hp_fe_values.get_fe_collection(),
+                     data.hp_fe_values.get_quadrature_collection(),
+                     data.hp_fe_values.get_update_flags())
+      {}
+
+      hp::FEValues<dim>               hp_fe_values;
+    };
+  }
+
+  namespace Copy
+  {
+    struct Data
+    {
+      Data (bool second_test = false)
+        :
+        second_test (second_test)
+      {}
+
+      Data (const Data &data)
+        :
+        second_test (data.second_test)
+      {}
+
+      std::vector<types::global_dof_index> local_dof_indices;
+      FullMatrix<double> local_matrix;
+      Vector<double> local_rhs;
+      const bool second_test;
+    };
+  }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+  LaplaceProblem ();
+  ~LaplaceProblem ();
+
+  void run ();
+
+private:
+  void setup_system ();
+  void test_equality ();
+  void assemble_reference ();
+  void assemble_test ();
+  void assemble_test_1();
+  void assemble_test_2();
+  void solve ();
+  void create_coarse_grid ();
+  void postprocess ();
+
+  void local_assemble (const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                       Assembly::Scratch::Data<dim>  &scratch,
+                       Assembly::Copy::Data          &data);
+  void copy_local_to_global (const Assembly::Copy::Data &data);
+
+  std::vector<types::global_dof_index>
+  get_conflict_indices (typename hp::DoFHandler<dim>::active_cell_iterator const &cell) const;
+
+  Triangulation<dim>   triangulation;
+
+  hp::DoFHandler<dim>      dof_handler;
+  hp::FECollection<dim>    fe_collection;
+  hp::QCollection<dim>     quadrature_collection;
+  hp::QCollection<dim-1>   face_quadrature_collection;
+
+  ConstraintMatrix     constraints;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> reference_matrix;
+  SparseMatrix<double> test_matrix;
+  SparseMatrix<double> test_matrix_2;
+
+  Vector<double>       solution;
+  Vector<double>       reference_rhs;
+  Vector<double>       test_rhs;
+  Vector<double>       test_rhs_2;
+
+  std::vector<std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> > graph;
+
+  const unsigned int max_degree;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                            const unsigned int  /*component*/) const
+{
+  double sum = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    sum += std::sin(deal_II_numbers::PI*p[d]);
+  return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+                           const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+  :
+  dof_handler (triangulation),
+  max_degree (5)
+{
+  if (dim == 2)
+    for (unsigned int degree=2; degree<=max_degree; ++degree)
+      {
+        fe_collection.push_back (FE_Q<dim>(degree));
+        quadrature_collection.push_back (QGauss<dim>(degree+1));
+        face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
+  else
+    for (unsigned int degree=1; degree<max_degree-1; ++degree)
+      {
+        fe_collection.push_back (FE_Q<dim>(degree));
+        quadrature_collection.push_back (QGauss<dim>(degree+1));
+        face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+  dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (typename hp::DoFHandler<dim>::active_cell_iterator const &cell) const
+{
+  std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+  cell->get_dof_indices(local_dof_indices);
+
+  constraints.resolve_indices(local_dof_indices);
+  return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  reference_matrix.clear();
+  test_matrix.clear();
+  dof_handler.distribute_dofs (fe_collection);
+
+  solution.reinit (dof_handler.n_dofs());
+  reference_rhs.reinit (dof_handler.n_dofs());
+  test_rhs.reinit (dof_handler.n_dofs());
+  test_rhs_2.reinit (dof_handler.n_dofs());
+
+  constraints.clear ();
+
+  // add a constraint that expands over most of the domain, which should
+  // create many conflicts
+  constraints.add_line(3);
+  for (unsigned int i=5; i<dof_handler.n_dofs(); i+=11)
+    constraints.add_entry(3, i, -0.5);
+
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+  // add boundary conditions as inhomogeneous constraints here, do it after
+  // having added the hanging node constraints in order to be consistent and
+  // skip dofs that are already constrained (i.e., are hanging nodes on the
+  // boundary in 3D). In contrast to step-27, we choose a sine function.
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+  constraints.close ();
+
+  graph = GraphColoring::make_graph_coloring(dof_handler.begin_active(),dof_handler.end(),
+                                             static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+                                             (typename hp::DoFHandler<dim>::active_cell_iterator const &)> >
+                                             (std_cxx1x::bind(&LaplaceProblem<dim>::get_conflict_indices, this,std_cxx1x::_1)));
+
+
+  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+                                       dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                   constraints, false);
+  sparsity_pattern.copy_from (csp);
+
+  reference_matrix.reinit (sparsity_pattern);
+  test_matrix.reinit (sparsity_pattern);
+  test_matrix_2.reinit (sparsity_pattern);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                                     Assembly::Scratch::Data<dim>  &scratch,
+                                     Assembly::Copy::Data          &data)
+{
+  const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+  data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+  data.local_matrix = 0;
+
+  data.local_rhs.reinit (dofs_per_cell);
+  data.local_rhs = 0;
+
+  scratch.hp_fe_values.reinit (cell);
+
+  const FEValues<dim> &fe_values = scratch.hp_fe_values.get_present_fe_values ();
+
+  const RightHandSide<dim> rhs_function;
+
+  for (unsigned int q_point=0;
+       q_point<fe_values.n_quadrature_points;
+       ++q_point)
+    {
+      const double scale_mat = data.second_test ? numbers::PI : 1.;
+      const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            data.local_matrix(i,j) += (scale_mat *
+                                       fe_values.shape_grad(i,q_point) *
+                                       fe_values.shape_grad(j,q_point) *
+                                       fe_values.JxW(q_point));
+
+            data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+                                  rhs_value * scale_mat *
+                                  fe_values.JxW(q_point));
+        }
+    }
+
+  data.local_dof_indices.resize (dofs_per_cell);
+  cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::copy_local_to_global (const Assembly::Copy::Data &data)
+{
+  if (data.second_test)
+    constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+                                           data.local_dof_indices,
+                                           test_matrix_2, test_rhs_2);
+  else
+    constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+                                           data.local_dof_indices,
+                                           test_matrix, test_rhs);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  Assembly::Copy::Data copy_data;
+  Assembly::Scratch::Data<dim> assembly_data(fe_collection, quadrature_collection);
+
+  for (unsigned int color=0; color<graph.size(); ++color)
+    for (typename std::vector<typename hp::DoFHandler<dim>::active_cell_iterator>::const_iterator p = graph[color].begin();
+         p != graph[color].end(); ++p)
+      {
+        local_assemble(*p, assembly_data, copy_data);
+        copy_local_to_global(copy_data);
+      }
+
+  reference_matrix.add(1., test_matrix);
+  reference_rhs = test_rhs;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test_1 ()
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe_collection, quadrature_collection),
+         Assembly::Copy::Data (),
+         1);
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test_2 ()
+{
+  test_matrix_2 = 0;
+  test_rhs_2 = 0;
+
+  WorkStream::
+    run (graph,
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          local_assemble,
+                          this,
+                          std_cxx1x::_1,
+                          std_cxx1x::_2,
+                          std_cxx1x::_3),
+         std_cxx1x::bind (&LaplaceProblem<dim>::
+                          copy_local_to_global,
+                          this,
+                          std_cxx1x::_1),
+         Assembly::Scratch::Data<dim>(fe_collection, quadrature_collection),
+         Assembly::Copy::Data (true),
+         1);
+}
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test()
+{
+  // start two tasks that each run an assembly
+  Threads::TaskGroup<void> tasks;
+  tasks += Threads::new_task (&LaplaceProblem<dim>::assemble_test_1,
+                              *this);
+  tasks += Threads::new_task (&LaplaceProblem<dim>::assemble_test_2,
+                              *this);
+  tasks.join_all();
+
+  test_matrix.add(-1, reference_matrix);
+  const double frobenius_norm = test_matrix.frobenius_norm();
+
+  // the data should add up exactly (unfortunately, there is some roundoff due
+  // to the cell similarity detection, but there should not be any similarity
+  // for the hypershell geometry)
+  deallog << "log error in matrix norm: " << std::log(frobenius_norm) << std::endl;
+  test_rhs.add(-1., reference_rhs);
+  deallog << "log error in vector norm: " << std::log(test_rhs.l2_norm()) << std::endl;
+
+  test_matrix_2.add(-numbers::PI, reference_matrix);
+  deallog << "absolute error matrix assembly 2: " << test_matrix_2.frobenius_norm() << std::endl;
+  test_rhs_2.add(-numbers::PI, reference_rhs);
+  deallog << "absolute error vector assembly 2: " << test_rhs_2.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+    estimated_error_per_cell(i) = i;
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                   estimated_error_per_cell,
+                                                   0.3, 0.03);
+  triangulation.execute_coarsening_and_refinement ();
+
+  for (typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
+       cell != dof_handler.end(); ++cell)
+    cell->set_active_fe_index (rand() % fe_collection.size());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_shell(triangulation,
+                                     Point<dim>(),
+                                     0.5, 1., (dim==3) ? 96 : 12, false);
+          triangulation.refine_global(2);
+        }
+
+      setup_system ();
+
+      assemble_reference ();
+      assemble_test ();
+
+      if (cycle < 2)
+        postprocess ();
+    }
+}
+
+
+
+int main ()
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-8);
+
+  {
+    deallog.push("2d");
+    LaplaceProblem<2> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/deal.II/assemble_matrix_parallel_04.output b/tests/deal.II/assemble_matrix_parallel_04.output
new file mode 100644
index 0000000000..627a39f770
--- /dev/null
+++ b/tests/deal.II/assemble_matrix_parallel_04.output
@@ -0,0 +1,13 @@
+
+DEAL:2d::log error in matrix norm: -inf
+DEAL:2d::log error in vector norm: -inf
+DEAL:2d::absolute error matrix assembly 2: 0
+DEAL:2d::absolute error vector assembly 2: 0
+DEAL:2d::log error in matrix norm: -inf
+DEAL:2d::log error in vector norm: -inf
+DEAL:2d::absolute error matrix assembly 2: 0
+DEAL:2d::absolute error vector assembly 2: 0
+DEAL:2d::log error in matrix norm: -inf
+DEAL:2d::log error in vector norm: -inf
+DEAL:2d::absolute error matrix assembly 2: 0
+DEAL:2d::absolute error vector assembly 2: 0