From: Denis Davydov <davydden@gmail.com>
Date: Thu, 23 Mar 2017 16:48:47 +0000 (+0100)
Subject: add a function to get DoFs that have support only within the given domain
X-Git-Tag: v9.0.0-rc1~1688^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c0c9940bc2916998fd2cbbd47bd9d42a09df1f06;p=dealii.git

add a function to get DoFs that have support only within the given domain
---

diff --git a/doc/doxygen/images/extract_dofs_with_support_on.png b/doc/doxygen/images/extract_dofs_with_support_on.png
new file mode 100644
index 0000000000..9d6c3328e9
Binary files /dev/null and b/doc/doxygen/images/extract_dofs_with_support_on.png differ
diff --git a/doc/news/changes/minor/20170324DenisDavydov b/doc/news/changes/minor/20170324DenisDavydov
new file mode 100644
index 0000000000..462f1a659c
--- /dev/null
+++ b/doc/news/changes/minor/20170324DenisDavydov
@@ -0,0 +1,4 @@
+New: DoFTools::extract_dofs_with_support_on() returns a set of degrees of
+freedom that have support on a subset of locally owned cells.
+<br>
+(Denis Davydov, 2017/03/24)
diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h
index aa4d0ababb..d0eead0241 100644
--- a/include/deal.II/dofs/dof_tools.h
+++ b/include/deal.II/dofs/dof_tools.h
@@ -1436,6 +1436,39 @@ namespace DoFTools
                                          std::vector<bool>      &selected_dofs,
                                          const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());
 
+  /**
+   * Extract all degrees of freedom (DoFs) that have support only within cells,
+   * for which @p predicate is <code>true</code>. The result is returned as
+   * an IndexSet.
+   *
+   * Consider the following FE space where predicate returns <code>true</code>
+   * for all cells on the left side:
+   *
+   * @image html extract_dofs_with_support_on.png
+   *
+   * This functions will return the union of all DoFs on those cells minus
+   * DoF 11, 13, 2 and 0; the result will be <code>[9,10], 12, [14,38]</code>.
+   * On the image above the returned DoFs are separated from the rest by the
+   * red line
+   *
+   * Essentially, the question this functions answers is the following:
+   * Given a subdomain with associated DoFs, what is the largest subset of
+   * these DoFs that are allowed to be non-zero such that after calling
+   * ConstraintMatrix::distribute() the resulting solution vector will have
+   * support only within the given domain. Here @p cm is the ConstraintMatrix
+   * containing hanging nodes constraints.
+   *
+   * In case of parallel::distributed::Triangulation @p predicate will be called
+   * only for locally owned and ghost cells. The resulting index set may contain
+   * DoFs that are associated to the locally owned or ghost cells, but are not
+   * owned by the current MPI core.
+   */
+  template <typename DoFHandlerType>
+  IndexSet
+  extract_dofs_with_support_on(const DoFHandlerType &dof_handler,
+                               const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> &predicate,
+                               const ConstraintMatrix &cm = ConstraintMatrix());
+
   /**
    * Extract a vector that represents the constant modes of the DoFHandler for
    * the components chosen by <tt>component_mask</tt> (see
diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc
index 19314c16c1..796335ef7f 100644
--- a/source/dofs/dof_tools.cc
+++ b/source/dofs/dof_tools.cc
@@ -765,6 +765,110 @@ namespace DoFTools
 
 
 
+  template <typename DoFHandlerType>
+  IndexSet
+  extract_dofs_with_support_on(const DoFHandlerType &dof_handler,
+                               const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> &predicate,
+                               const ConstraintMatrix &cm)
+  {
+    const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> predicate_local
+    = [=] (const typename DoFHandlerType::active_cell_iterator & cell) -> bool
+    {
+      return cell->is_locally_owned() && predicate(cell);
+    };
+
+    std::vector<types::global_dof_index> local_dof_indices;
+    local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
+
+    // Get all the dofs that live on the subdomain:
+    std::set<types::global_dof_index> predicate_dofs;
+
+    typename DoFHandlerType::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (!cell->is_artificial() && predicate(cell))
+        {
+          local_dof_indices.resize (cell->get_fe().dofs_per_cell);
+          cell->get_dof_indices (local_dof_indices);
+          predicate_dofs.insert(local_dof_indices.begin(),
+                                local_dof_indices.end());
+        }
+
+    // Get halo layer and accumulate its DoFs
+    std::set<types::global_dof_index> halo_dofs;
+
+    const std::vector<typename DoFHandlerType::active_cell_iterator> halo_cells =
+      GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local);
+    for (typename std::vector<typename DoFHandlerType::active_cell_iterator>::const_iterator it = halo_cells.begin();
+         it != halo_cells.end(); ++it)
+      {
+        // skip halo cells that satisfy the given predicate and thereby
+        // shall be a part of the index set on another MPI core.
+        // Those could only be ghost cells with p::d::Tria.
+        if (predicate(*it))
+          {
+            Assert ((*it)->is_ghost(),
+                    ExcInternalError());
+            continue;
+          }
+
+        const unsigned int dofs_per_cell = (*it)->get_fe().dofs_per_cell;
+        local_dof_indices.resize (dofs_per_cell);
+        (*it)->get_dof_indices (local_dof_indices);
+        halo_dofs.insert(local_dof_indices.begin(),
+                         local_dof_indices.end());
+      }
+
+    // A complication coming from the fact that DoFs living on locally
+    // owned cells which satisfy predicate may participate in constrains for
+    // DoFs outside of this region.
+    if (cm.n_constraints() > 0)
+      {
+        const std::vector<std::pair<types::global_dof_index,double> > *line_ptr;
+
+        // Remove DoFs that are part of constrains for DoFs outside
+        // of predicate. Since we will subtract halo_dofs from predicate_dofs,
+        // achieve this by extending halo_dofs with DoFs to which
+        // halo_dofs are constrained.
+        std::set<types::global_dof_index> extra_halo;
+        for (std::set<types::global_dof_index>::iterator it = halo_dofs.begin();
+             it != halo_dofs.end(); ++it)
+          {
+            line_ptr = cm.get_constraint_entries(*it);
+            // if halo DoF is constrained, add all DoFs to which it's constrained.
+            if (line_ptr!=NULL)
+              {
+                const unsigned int line_size = line_ptr->size();
+                for (unsigned int j=0; j<line_size; ++j)
+                  extra_halo.insert((*line_ptr)[j].first);
+              }
+          }
+
+        halo_dofs.insert(extra_halo.begin(),
+                         extra_halo.end());
+      }
+
+    // Rework our std::set's into IndexSet and subtract halo DoFs from the
+    // predicate_dofs:
+    IndexSet support_set(dof_handler.n_dofs());
+    support_set.add_indices(predicate_dofs.begin(), predicate_dofs.end());
+    support_set.compress();
+
+    IndexSet halo_set(dof_handler.n_dofs());
+    halo_set.add_indices(halo_dofs.begin(), halo_dofs.end());
+    halo_set.compress();
+
+    support_set.subtract_set(halo_set);
+
+    // we intentionally do not want to limit the output to locally owned DoFs.
+    // In other words, the output may contain DoFs that are on the locally
+    // owned cells, but are not owned by this MPI core.
+
+    return support_set;
+  }
+
+
 
   namespace internal
   {
diff --git a/source/dofs/dof_tools.inst.in b/source/dofs/dof_tools.inst.in
index a11b5b3077..08d536ae16 100644
--- a/source/dofs/dof_tools.inst.in
+++ b/source/dofs/dof_tools.inst.in
@@ -160,6 +160,13 @@ for (deal_II_dimension : DIMENSIONS)
      std::vector<bool>                        &,
      const std::set<types::boundary_id> &);
 
+    template
+    IndexSet
+    DoFTools::extract_dofs_with_support_on<DoFHandler<deal_II_dimension> >
+    (const DoFHandler<deal_II_dimension> &,
+     const std::function< bool(const typename DoFHandler<deal_II_dimension>::active_cell_iterator &)> &,
+     const ConstraintMatrix&);
+
     template
     void
     DoFTools::extract_hanging_node_dofs
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc
new file mode 100644
index 0000000000..741020cd8c
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc
@@ -0,0 +1,137 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+//
+
+// Test DoFTools::extract_dofs_with_support_on()
+
+#include "../tests.h"
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+template<int dim>
+bool
+pred_d(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) < 0.49 &&
+          cell->center()(1) < 0.49);
+}
+
+template<int dim>
+bool
+pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) < 0.49 &&
+          cell->center()(0) > 0.25 &&
+          cell->center()(1) < 0.49);
+}
+
+
+
+template <int dim>
+void test (const unsigned int flag)
+{
+
+  // Setup system
+  Triangulation<dim> triangulation;
+
+  GridGenerator::hyper_rectangle (triangulation,
+                                  Point<dim>(0,0),
+                                  Point<dim>(1,1));
+
+  if (flag == 0)
+    triangulation.refine_global(2);
+  else
+    triangulation.refine_global(1);
+
+  // Extra refinement to generate hanging nodes
+  for (typename Triangulation<dim>::active_cell_iterator
+       cell = triangulation.begin_active();
+       cell != triangulation.end(); ++cell)
+    if ( (flag==1 && pred_d<dim>(cell)) ||
+         (flag==2 && !pred_d<dim>(cell))  )
+      cell->set_refine_flag ();
+
+  triangulation.prepare_coarsening_and_refinement();
+  triangulation.execute_coarsening_and_refinement ();
+
+  DoFHandler<dim> dh (triangulation);
+
+  FE_Q<dim> fe(2);
+  dh.distribute_dofs (fe);
+
+  ConstraintMatrix cm;
+  DoFTools::make_hanging_node_constraints(dh,cm);
+
+  IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d<dim>, cm);
+  support.print(deallog);
+
+  // print grid and DoFs for visual inspection
+  if (false)
+    {
+      std::cout << "-------------------- " << Utilities::int_to_string(dim) << Utilities::int_to_string(flag) << std::endl;
+      cm.print(std::cout);
+
+      std::map<types::global_dof_index, Point<dim> > support_points;
+      MappingQ1<dim> mapping;
+      DoFTools::map_dofs_to_support_points(mapping, dh, support_points);
+
+      const std::string filename =
+        "grid" + Utilities::int_to_string(dim) + Utilities::int_to_string(flag) + ".gp";
+      std::ofstream f(filename.c_str());
+
+      f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl
+        << "set output \"grid" << Utilities::int_to_string(dim) << Utilities::int_to_string(flag) << ".png\"" << std::endl
+        << "set size square" << std::endl
+        << "set view equal xy" << std::endl
+        << "unset xtics" << std::endl
+        << "unset ytics" << std::endl
+        << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" << std::endl;
+      GridOut().write_gnuplot (triangulation, f);
+      f << "e" << std::endl;
+
+      DoFTools::write_gnuplot_dof_support_point_info(f,
+                                                     support_points);
+      f << "e" << std::endl;
+    }
+
+  dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+  initlog();
+
+  test<2>(0);
+  test<2>(1);
+  test<2>(2);
+
+  return 0;
+}
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output
new file mode 100644
index 0000000000..d326d0283a
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output
@@ -0,0 +1,4 @@
+
+DEAL::{[0,8], [12,14], [17,18], 20, 24}
+DEAL::{[21,42]}
+DEAL::{0, 4, 6, 8}
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc
new file mode 100644
index 0000000000..d74ba742fb
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc
@@ -0,0 +1,185 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+//
+
+// similar to dof_tools_23, but for the p::d::Tria. As an MPI test, make sure
+// that the total number of shape functions that are non-zero within the domain
+// is the same.
+
+#include "../tests.h"
+#include <deal.II/distributed/tria.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+template<int dim>
+bool
+pred_d(const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) < 0.5 &&
+          cell->center()(1) < 0.5);
+}
+
+
+std::string output_name(const unsigned int flag, const unsigned int subdomain)
+{
+  return "output" + Utilities::int_to_string(flag) + "_" +
+         Utilities::int_to_string(subdomain) +
+         ".vtu";
+}
+
+
+template <int dim>
+void test (const unsigned int flag)
+{
+  // Setup system
+  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+
+  GridGenerator::hyper_rectangle (triangulation,
+                                  Point<dim>(0,0),
+                                  Point<dim>(1,1));
+
+  if (flag == 0)
+    triangulation.refine_global(2);
+  else
+    triangulation.refine_global(1);
+
+  // Extra refinement to generate hanging nodes
+  for (typename parallel::distributed::Triangulation<dim>::active_cell_iterator
+       cell = triangulation.begin_active();
+       cell != triangulation.end(); ++cell)
+    if (cell->is_locally_owned() &&
+        ((flag==1 &&  pred_d<dim>(cell)) ||
+         (flag==2 && !pred_d<dim>(cell)))   )
+      cell->set_refine_flag ();
+
+  triangulation.prepare_coarsening_and_refinement();
+  triangulation.execute_coarsening_and_refinement ();
+
+  if (flag > 0)
+    triangulation.refine_global(1);
+
+  DoFHandler<dim> dh (triangulation);
+
+  FE_Q<dim> fe(2);
+  dh.distribute_dofs (fe);
+
+  IndexSet locally_relevant_set;
+  DoFTools::extract_locally_relevant_dofs (dh,
+                                           locally_relevant_set);
+
+  ConstraintMatrix cm;
+  cm.reinit (locally_relevant_set);
+  DoFTools::make_hanging_node_constraints(dh,cm);
+  cm.close ();
+
+  const IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d<dim>, cm);
+  const IndexSet support_local = support & dh.locally_owned_dofs();
+
+  deallog << support.n_elements() << std::endl;
+
+  // now accumulate the number of indices and make sure it's the same for
+  // various runs with different number of MPI cores
+  const unsigned int dofs_support = Utilities::MPI::sum(support_local.n_elements(), MPI_COMM_WORLD);
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      deallog << "accumulated: " << std::endl;
+      deallog << dofs_support << std::endl;
+    }
+
+  // print grid and DoFs for visual inspection
+  if (false)
+    {
+      DataOut<dim> data_out;
+      data_out.attach_dof_handler (dh);
+
+      std::vector<LinearAlgebra::distributed::Vector<double> > shape_functions(dh.n_dofs());
+      for (unsigned int i = 0; i < dh.n_dofs(); ++i)
+        {
+          LinearAlgebra::distributed::Vector<double> &s = shape_functions[i];
+          s.reinit(dh.locally_owned_dofs(),
+                   locally_relevant_set,
+                   MPI_COMM_WORLD);
+          s = 0.;
+          if (dh.locally_owned_dofs().is_element(i))
+            s[i] = 1.0;
+          s.compress(VectorOperation::insert);
+          cm.distribute(s);
+          s.update_ghost_values();
+
+          data_out.add_data_vector (s,
+                                    std::string("N_") +
+                                    Utilities::int_to_string(i));
+        }
+
+      Vector<float> subdomain (triangulation.n_active_cells());
+      for (unsigned int i=0; i<subdomain.size(); ++i)
+        subdomain(i) = triangulation.locally_owned_subdomain();
+      data_out.add_data_vector (subdomain, "subdomain");
+      data_out.build_patches ();
+
+      const std::string filename = output_name(flag, triangulation.locally_owned_subdomain());
+
+      std::ofstream output (filename.c_str());
+      data_out.write_vtu (output);
+
+      if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+        {
+          std::vector<std::string> filenames;
+          for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
+            filenames.push_back (output_name(flag, i));
+
+          const std::string master_name = "output" + Utilities::int_to_string(flag) + ".pvtu";
+          std::ofstream pvtu_master (master_name.c_str());
+          data_out.write_pvtu_record (pvtu_master, filenames);
+        }
+
+    }
+
+  dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());
+  const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+  MPILogInitAll log;
+
+  deallog.push(Utilities::int_to_string(myid));
+
+  test<2>(0);
+  test<2>(1);
+  test<2>(2);
+
+  deallog.pop();
+}
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output
new file mode 100644
index 0000000000..ad14a0a639
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output
@@ -0,0 +1,10 @@
+
+DEAL:0:0::16
+DEAL:0:0::accumulated: 
+DEAL:0:0::16
+DEAL:0:0::76
+DEAL:0:0::accumulated: 
+DEAL:0:0::76
+DEAL:0:0::16
+DEAL:0:0::accumulated: 
+DEAL:0:0::16
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output
new file mode 100644
index 0000000000..24d5d3f17a
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output
@@ -0,0 +1,20 @@
+
+DEAL:0:0::16
+DEAL:0:0::accumulated: 
+DEAL:0:0::16
+DEAL:0:0::61
+DEAL:0:0::accumulated: 
+DEAL:0:0::76
+DEAL:0:0::16
+DEAL:0:0::accumulated: 
+DEAL:0:0::16
+
+DEAL:1:1::21
+DEAL:1:1::64
+DEAL:1:1::15
+
+
+DEAL:2:2::9
+DEAL:2:2::27
+DEAL:2:2::9
+
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output
new file mode 100644
index 0000000000..4a18c5f509
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output
@@ -0,0 +1,30 @@
+
+DEAL:0:0::16
+DEAL:0:0::accumulated: 
+DEAL:0:0::16
+DEAL:0:0::49
+DEAL:0:0::accumulated: 
+DEAL:0:0::76
+DEAL:0:0::16
+DEAL:0:0::accumulated: 
+DEAL:0:0::16
+
+DEAL:1:1::15
+DEAL:1:1::69
+DEAL:1:1::15
+
+
+DEAL:2:2::0
+DEAL:2:2::46
+DEAL:2:2::15
+
+
+DEAL:3:3::15
+DEAL:3:3::45
+DEAL:3:3::9
+
+
+DEAL:4:4::9
+DEAL:4:4::9
+DEAL:4:4::0
+
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc
new file mode 100644
index 0000000000..e19e4ec036
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc
@@ -0,0 +1,143 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+//
+
+// same as dof_tools_32 test DoFTools::extract_dofs_with_support_on() but
+// for a slightly different configuration
+
+#include "../tests.h"
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+template<int dim>
+bool
+pred_left(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) < 0.49);
+}
+
+template<int dim>
+bool
+pred_right(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) > 0.51);
+}
+
+
+template<int dim>
+bool
+pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) < 0.49 &&
+          cell->center()(1) < 0.49) ||
+         (cell->center()(0) > 0.49 &&
+          cell->center()(1) > 0.49);
+}
+
+
+
+template <int dim>
+void test (const bool left = true)
+{
+
+  // Setup system
+  Triangulation<dim> triangulation;
+
+  GridGenerator::hyper_rectangle (triangulation,
+                                  Point<dim>(0,0),
+                                  Point<dim>(1,1));
+
+  triangulation.refine_global(1);
+
+  // Extra refinement to generate hanging nodes
+  for (typename Triangulation<dim>::active_cell_iterator
+       cell = triangulation.begin_active();
+       cell != triangulation.end(); ++cell)
+    if (pred_r<dim>(cell))
+      cell->set_refine_flag ();
+
+  triangulation.prepare_coarsening_and_refinement();
+  triangulation.execute_coarsening_and_refinement ();
+
+  DoFHandler<dim> dh (triangulation);
+
+  FE_Q<dim> fe(2);
+  dh.distribute_dofs (fe);
+
+  ConstraintMatrix cm;
+  DoFTools::make_hanging_node_constraints(dh,cm);
+
+  IndexSet support = left ?
+                     DoFTools::extract_dofs_with_support_on(dh, pred_left<dim>, cm) :
+                     DoFTools::extract_dofs_with_support_on(dh, pred_right<dim>, cm);
+  support.print(deallog);
+
+  // print grid and DoFs for visual inspection
+  if (false)
+    {
+      std::cout << "-------------------- " << Utilities::int_to_string(dim) << std::endl;
+      cm.print(std::cout);
+
+      std::map<types::global_dof_index, Point<dim> > support_points;
+      MappingQ1<dim> mapping;
+      DoFTools::map_dofs_to_support_points(mapping, dh, support_points);
+
+      const std::string filename =
+        "grid" + Utilities::int_to_string(dim) + ".gp";
+      std::ofstream f(filename.c_str());
+
+      f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl
+        << "set output \"grid" << Utilities::int_to_string(dim) << ".png\"" << std::endl
+        << "set size square" << std::endl
+        << "set view equal xy" << std::endl
+        << "unset xtics" << std::endl
+        << "unset ytics" << std::endl
+        << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" << std::endl;
+      GridOut().write_gnuplot (triangulation, f);
+      f << "e" << std::endl;
+
+      DoFTools::write_gnuplot_dof_support_point_info(f,
+                                                     support_points);
+      f << "e" << std::endl;
+    }
+
+  dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+  initlog();
+
+  test<2>(true);
+  test<2>(false);
+
+  return 0;
+}
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output
new file mode 100644
index 0000000000..76ee14d0e0
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output
@@ -0,0 +1,3 @@
+
+DEAL::{[9,10], 12, [14,38]}
+DEAL::{1, 3, [5,8], [39,60]}
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc
new file mode 100644
index 0000000000..7126ab0558
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc
@@ -0,0 +1,327 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+//
+
+// test DoFTools::extract_dofs_with_support_on() for calculation of the RHS
+// function in parallel
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+std::string output_name(const unsigned int subdomain)
+{
+  return "output_" +
+         Utilities::int_to_string(subdomain) +
+         ".vtu";
+}
+
+template<int dim>
+bool
+pred_d(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) < 0.49);
+}
+
+template<int dim>
+bool
+pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  return (cell->center()(0) < 0.49 &&
+          cell->center()(1) < 0.49) ||
+         (cell->center()(0) > 0.49 &&
+          cell->center()(1) > 0.49);
+}
+
+
+
+template <int dim>
+void test ()
+{
+
+  // Setup system
+  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+
+  GridGenerator::hyper_rectangle (triangulation,
+                                  Point<dim>(0,0),
+                                  Point<dim>(1,1));
+
+  triangulation.refine_global(1);
+
+  // Extra refinement to generate hanging nodes
+  for (typename Triangulation<dim>::active_cell_iterator
+       cell = triangulation.begin_active();
+       cell != triangulation.end(); ++cell)
+    if (cell->is_locally_owned() && pred_r<dim>(cell))
+      cell->set_refine_flag ();
+
+  triangulation.prepare_coarsening_and_refinement();
+  triangulation.execute_coarsening_and_refinement ();
+
+  DoFHandler<dim> dh (triangulation);
+
+  FE_Q<dim> fe(2);
+  dh.distribute_dofs (fe);
+
+  const IndexSet &locally_owned_set = dh.locally_owned_dofs();
+  IndexSet locally_relevant_set;
+  DoFTools::extract_locally_relevant_dofs (dh,
+                                           locally_relevant_set);
+
+  ConstraintMatrix cm;
+  cm.reinit (locally_relevant_set);
+  DoFTools::make_hanging_node_constraints(dh,cm);
+  cm.close ();
+
+  std::vector<types::global_dof_index> local_dof_indices(fe.dofs_per_cell);
+
+  // get support on the predicate
+  IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d<dim>, cm);
+  IndexSet local_support = support & locally_owned_set;
+
+  // rhs vectors:
+  LinearAlgebra::distributed::Vector<double> sparse_rhs;
+  LinearAlgebra::distributed::Vector<double> rhs;
+  sparse_rhs.reinit(locally_owned_set, locally_relevant_set, MPI_COMM_WORLD);
+  rhs.reinit(       locally_owned_set, locally_relevant_set, MPI_COMM_WORLD);
+
+  sparse_rhs.zero_out_ghosts();
+  rhs.zero_out_ghosts();
+
+  // assemble RHS which has a local support:
+  const std::function< double(const Point<dim> &)> rhs_func
+  = [=] (const Point<dim> &p) -> double
+  {
+    return p[0] > 0.5 ? 0. : 0.5 - p[0];
+  };
+
+  Vector<double> local_rhs(fe.dofs_per_cell);
+  QGauss<dim> quadrature(3);
+  FEValues<dim> fe_values(fe, quadrature,
+                          update_values | update_JxW_values | update_quadrature_points);
+  for (typename DoFHandler<dim>::active_cell_iterator
+       cell = dh.begin_active();
+       cell != dh.end(); ++cell)
+    if (cell->is_locally_owned() && pred_d<dim>(cell))
+      {
+        fe_values.reinit(cell);
+        cell->get_dof_indices (local_dof_indices);
+
+        const std::vector<Point<dim>> &q_points = fe_values.get_quadrature_points();
+
+        local_rhs = 0.;
+        for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+          for (unsigned int q = 0; q < quadrature.size(); ++q)
+            local_rhs[i] += fe_values.shape_value (i, q) *
+                            rhs_func(q_points[q]) *
+                            fe_values.JxW (q);
+
+        cm.distribute_local_to_global (local_rhs,
+                                       local_dof_indices,
+                                       rhs);
+
+        // copy-paste of CM distribute_local_to_global and
+        // add is_element() checks:
+        auto local_vector_begin   = local_rhs.begin();
+        const auto local_vector_end = local_rhs.end();
+        auto local_indices_begin = local_dof_indices.begin();
+        const std::vector<std::pair<types::global_dof_index,double> > *line_ptr;
+        for ( ; local_vector_begin != local_vector_end;
+              ++local_vector_begin, ++local_indices_begin)
+          {
+            line_ptr = cm.get_constraint_entries(*local_indices_begin);
+            if (line_ptr==NULL) // unconstrained
+              {
+                if (support.is_element(*local_indices_begin))
+                  sparse_rhs(*local_indices_begin) += *local_vector_begin;
+              }
+            else
+              {
+                const unsigned int line_size = line_ptr->size();
+                for (unsigned int j=0; j<line_size; ++j)
+                  if (support.is_element((*line_ptr)[j].first))
+                    sparse_rhs((*line_ptr)[j].first)
+                    += *local_vector_begin * (*line_ptr)[j].second;
+              }
+          }
+        /*
+        cm.distribute_local_to_global (local_rhs,
+                                       local_dof_indices,
+                                       sparse_rhs);
+        */
+      }
+
+  rhs.compress(VectorOperation::add);
+  sparse_rhs.compress(VectorOperation::add);
+
+  rhs.update_ghost_values();
+  sparse_rhs.update_ghost_values();
+
+  rhs.add(-1., sparse_rhs);
+
+  // print grid and DoFs for visual inspection
+  if (false)
+    {
+      const unsigned int this_mpi_process = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+      const unsigned int n_mpi_processes = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+      for (unsigned int i = 0; i < n_mpi_processes; ++i)
+        {
+          MPI_Barrier(MPI_COMM_WORLD);
+          if (i == this_mpi_process)
+            {
+              std::cout << "-------------------- " << this_mpi_process << std::endl;
+              cm.print(std::cout);
+
+              std::cout << "local support:" << std::endl;
+              local_support.print(std::cout);
+              std::cout << "support:" << std::endl;
+              support.print(std::cout);
+              std::cout << "locally owned:" << std::endl;
+              locally_owned_set.print(std::cout);
+              std::cout << "locally relevant set:" << std::endl;
+              locally_relevant_set.print(std::cout);
+            }
+        }
+
+      {
+        std::map<types::global_dof_index, Point<dim> > support_points;
+        MappingQ1<dim> mapping;
+        DoFTools::map_dofs_to_support_points(mapping, dh, support_points);
+
+        const std::string filename =
+          "grid" + dealii::Utilities::int_to_string(n_mpi_processes) + dealii::Utilities::int_to_string(this_mpi_process);
+        std::ofstream f((filename+".gp").c_str());
+
+        f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl
+          << "set output \"" << filename << ".png\"" << std::endl
+          << "set size square" << std::endl
+          << "set view equal xy" << std::endl
+          << "unset xtics" << std::endl
+          << "unset ytics" << std::endl
+          << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" << std::endl;
+        GridOut().write_gnuplot (triangulation, f);
+        f << "e" << std::endl;
+
+        DoFTools::write_gnuplot_dof_support_point_info(f,
+                                                       support_points);
+        f << "e" << std::endl;
+
+      }
+
+      {
+        DataOut<dim> data_out;
+        data_out.attach_dof_handler (dh);
+
+        std::vector<LinearAlgebra::distributed::Vector<double> > shape_functions(dh.n_dofs());
+        for (unsigned int i = 0; i < dh.n_dofs(); ++i)
+          {
+            LinearAlgebra::distributed::Vector<double> sl(locally_owned_set, MPI_COMM_WORLD);
+            sl = 0.;
+            if (locally_owned_set.is_element(i))
+              sl[i] = 1.0;
+            cm.distribute(sl);
+
+            LinearAlgebra::distributed::Vector<double> &s = shape_functions[i];
+            s.reinit(locally_owned_set,
+                     locally_relevant_set,
+                     MPI_COMM_WORLD);
+            s = 0.;
+            s = sl;
+
+            data_out.add_data_vector (s,
+                                      std::string("N_") +
+                                      dealii::Utilities::int_to_string(i));
+          }
+
+        Vector<float> subdomain (triangulation.n_active_cells());
+        for (unsigned int i=0; i<subdomain.size(); ++i)
+          subdomain(i) = triangulation.locally_owned_subdomain();
+        data_out.add_data_vector (subdomain, "subdomain");
+        data_out.build_patches ();
+
+        const std::string filename = output_name(triangulation.locally_owned_subdomain());
+
+        std::ofstream output (filename.c_str());
+        data_out.write_vtu (output);
+
+        if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+          {
+            std::vector<std::string> filenames;
+            for (unsigned int i=0; i<dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
+              filenames.push_back (output_name(i));
+
+            const std::string master_name = "output.pvtu";
+            std::ofstream pvtu_master (master_name.c_str());
+            data_out.write_pvtu_record (pvtu_master, filenames);
+          }
+      }
+    }
+
+  for (unsigned int i = 0; i < locally_owned_set.n_elements(); ++i)
+    {
+      const unsigned int ind = locally_owned_set.nth_index_in_set(i);
+      const double v = rhs[ind];
+      AssertThrow(std::abs(v) < 1e-12,
+                  ExcMessage(
+                    "Element " + std::to_string(ind) +
+                    " has an error " + std::to_string(v) +
+                    " on the process " + std::to_string(dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD))
+                  ));
+    }
+
+  if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    deallog << "Ok" << std::endl;
+
+
+  dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile,/*do not print job id*/false);
+  deallog.depth_console(0);
+  deallog.precision(4);
+
+  dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  test<2>();
+
+  return 0;
+}
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output
new file mode 100644
index 0000000000..04c93dbf6b
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output
@@ -0,0 +1 @@
+DEAL::Ok
diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output
new file mode 100644
index 0000000000..04c93dbf6b
--- /dev/null
+++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output
@@ -0,0 +1 @@
+DEAL::Ok