From 265ef03ac98d90806e68f81a10cfec87339377d9 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Mon, 11 Apr 2016 11:58:04 +0200
Subject: [PATCH] Fix Cuthill-McKee renumbering on large processor counts

---
 doc/news/changes.h             |  9 +++++++++
 source/dofs/dof_renumbering.cc | 21 ++++++++++++++-------
 source/dofs/dof_tools.inst.in  | 18 +++++++++++++-----
 3 files changed, 36 insertions(+), 12 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index fadcb84bdd..1df3d429ec 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -171,6 +171,15 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+ <li> Improved: DoFRenumbering::compute_Cuthill_McKee when used with
+ distributed triangulations contained parts that scaled as the global problem
+ size, rather than the processor-local size. This prevented its use with more
+ than a few hundred cores when hanging node constraints were activated. This
+ has been fixed.
+ <br>
+ (Martin Kronbichler, 2016/04/11)
+ </li>
+
  <li> New: added hessenberg_signal and krylov_space_signal to SolverGMRES.
  These signals allow to retrieve the Hessenberg matrix and the basis vectors
  generated by the Arnoldi algorithm.
diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc
index ae35017783..e076b8ada5 100644
--- a/source/dofs/dof_renumbering.cc
+++ b/source/dofs/dof_renumbering.cc
@@ -392,7 +392,12 @@ namespace DoFRenumbering
     // object will be empty and nothing happens
     ConstraintMatrix constraints;
     if (use_constraints)
-      DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+      {
+        IndexSet relevant_dofs;
+        DoFTools::extract_locally_relevant_dofs(dof_handler, relevant_dofs);
+        constraints.reinit(relevant_dofs);
+        DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+      }
     constraints.close ();
 
     const IndexSet locally_owned = dof_handler.locally_owned_dofs();
@@ -408,8 +413,7 @@ namespace DoFRenumbering
 
     // If the index set is not complete, need to get indices in local index
     // space.
-    if (locally_owned.n_elements() !=
-        locally_owned.size())
+    if (locally_owned.n_elements() != locally_owned.size())
       {
         // Create sparsity pattern from dsp by transferring its indices to
         // processor-local index space and doing Cuthill-McKee there
@@ -419,11 +423,14 @@ namespace DoFRenumbering
         for (unsigned int i=0; i<locally_owned.n_elements(); ++i)
           {
             const types::global_dof_index row = locally_owned.nth_index_in_set(i);
+            const unsigned int row_length = dsp.row_length(row);
             row_entries.clear();
-            for (DynamicSparsityPattern::iterator it =
-                   dsp.begin(row); it != dsp.end(row); ++it)
-              if (it->column() != row && locally_owned.is_element(it->column()))
-                row_entries.push_back(locally_owned.index_within_set(it->column()));
+            for (unsigned int j=0; j<row_length; ++j)
+              {
+                const unsigned int col = dsp.column_number(row, j);
+                if (col != row && locally_owned.is_element(col))
+                  row_entries.push_back(locally_owned.index_within_set(col));
+              }
             sparsity.add_entries(i, row_entries.begin(), row_entries.end(),
                                  true);
           }
diff --git a/source/dofs/dof_tools.inst.in b/source/dofs/dof_tools.inst.in
index feabaa7242..6446b65023 100644
--- a/source/dofs/dof_tools.inst.in
+++ b/source/dofs/dof_tools.inst.in
@@ -32,12 +32,20 @@ for (deal_II_dimension : DIMENSIONS)
 
 #if deal_II_dimension > 1
   template
-  void 
+  void
   DoFTools::extract_locally_relevant_dofs<DoFHandler<deal_II_dimension-1,deal_II_dimension > >
   (const DoFHandler<deal_II_dimension-1,deal_II_dimension > & dof_handler,
    IndexSet & dof_set);
 #endif
 
+#if deal_II_dimension > 2
+  template
+  void
+  DoFTools::extract_locally_relevant_dofs<DoFHandler<deal_II_dimension-2,deal_II_dimension > >
+  (const DoFHandler<deal_II_dimension-2,deal_II_dimension > & dof_handler,
+   IndexSet & dof_set);
+#endif
+
   template
   void DoFTools::make_vertex_patches (SparsityPattern&, const DoFHandler<deal_II_dimension>&,
   unsigned int, bool, bool, bool, bool);
@@ -207,7 +215,7 @@ void
 DoFTools::extract_locally_relevant_dofs<DoFHandler<deal_II_dimension> >
 (const DoFHandler<deal_II_dimension> & dof_handler,
  IndexSet & dof_set);
- 
+
 template
 void
 DoFTools::extract_locally_relevant_dofs<hp::DoFHandler<deal_II_dimension> >
@@ -257,7 +265,7 @@ void
 DoFTools::get_subdomain_association<hp::DoFHandler<deal_II_dimension> >
 (const hp::DoFHandler<deal_II_dimension> &dof_handler,
  std::vector<types::subdomain_id>           &subdomain_association);
- 
+
 template
 std::vector<IndexSet>
 DoFTools::locally_owned_dofs_per_subdomain<DoFHandler<deal_II_dimension> >
@@ -368,7 +376,7 @@ void
 DoFTools::get_subdomain_association<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> >
 (const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &dof_handler,
  std::vector<types::subdomain_id>           &subdomain_association);
- 
+
 template
 std::vector<IndexSet>
 DoFTools::locally_owned_dofs_per_subdomain<DoFHandler<deal_II_dimension,deal_II_dimension+1> >
@@ -422,7 +430,7 @@ void
 DoFTools::get_subdomain_association<hp::DoFHandler<1,3> >
 (const hp::DoFHandler<1,3> &dof_handler,
  std::vector<types::subdomain_id>           &subdomain_association);
- 
+
 template
 std::vector<IndexSet>
 DoFTools::locally_owned_dofs_per_subdomain<DoFHandler<1,3> >
-- 
2.39.5