From 5fe7213ec48e2032448511d532c8368c92850347 Mon Sep 17 00:00:00 2001
From: fsonner <florian.sonner@iwr.uni-heidelberg.de>
Date: Fri, 24 Apr 2015 16:56:20 +0200
Subject: [PATCH] Fix block-wise dof renumbering for multigrid

Re-adds a version of DoFRenumbering::block_wise working on a single
multigrid level. For this an internal helper function is changed to
work with non-active cells.
---
 include/deal.II/dofs/dof_renumbering.h | 12 ++++-
 source/dofs/dof_renumbering.cc         | 66 ++++++++++++++++----------
 2 files changed, 52 insertions(+), 26 deletions(-)

diff --git a/include/deal.II/dofs/dof_renumbering.h b/include/deal.II/dofs/dof_renumbering.h
index 74dc7ddbe9..be688f48dc 100644
--- a/include/deal.II/dofs/dof_renumbering.h
+++ b/include/deal.II/dofs/dof_renumbering.h
@@ -657,6 +657,15 @@ namespace DoFRenumbering
   void
   block_wise (DoFHandler<dim,spacedim> &dof_handler);
 
+  /**
+   * Sort the degrees of freedom by vector block. It does the same thing as the
+   * above function, only that it does this for one single level of a multi-
+   * level discretization. The non-multigrid part of the the DoFHandler is not
+   * touched.
+   */
+  template <int dim, int spacedim>
+  void
+  block_wise (DoFHandler<dim,spacedim> &dof_handler, unsigned int level);
 
   /**
    * Sort the degrees of freedom by block. It does the same thing as the above
@@ -687,7 +696,8 @@ namespace DoFRenumbering
   types::global_dof_index
   compute_block_wise (std::vector<types::global_dof_index> &new_dof_indices,
                       const ITERATOR &start,
-                      const ENDITERATOR &end);
+                      const ENDITERATOR &end,
+                      bool is_level_operation);
 
   /**
    * @}
diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc
index 1ff221061c..e317bec837 100644
--- a/source/dofs/dof_renumbering.cc
+++ b/source/dofs/dof_renumbering.cc
@@ -841,7 +841,7 @@ namespace DoFRenumbering
     const types::global_dof_index result =
       compute_block_wise<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,
       typename DoFHandler<dim,spacedim>::level_cell_iterator>
-      (renumbering, start, end);
+      (renumbering, start, end, false);
     if (result == 0)
       return;
 
@@ -879,7 +879,7 @@ namespace DoFRenumbering
     const types::global_dof_index result =
       compute_block_wise<dim, dim, typename hp::DoFHandler<dim>::active_cell_iterator,
       typename hp::DoFHandler<dim>::level_cell_iterator>(renumbering,
-                                                         start, end);
+                                                         start, end, false);
 
     if (result == 0)
       return;
@@ -892,25 +892,25 @@ namespace DoFRenumbering
 
 
 
-  template <int dim>
+  template <int dim, int spacedim>
   void
-  block_wise (DoFHandler<dim> &dof_handler, const unsigned int level)
+  block_wise (DoFHandler<dim,spacedim> &dof_handler, const unsigned int level)
   {
     Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
            ExcNotInitialized());
 
     std::vector<types::global_dof_index> renumbering (dof_handler.n_dofs(level),
-                                                      DoFHandler<dim>::invalid_dof_index);
+                                                      DoFHandler<dim, spacedim>::invalid_dof_index);
 
-    typename DoFHandler<dim>::level_cell_iterator
+    typename DoFHandler<dim, spacedim>::level_cell_iterator
     start =dof_handler.begin(level);
-    typename DoFHandler<dim>::level_cell_iterator
+    typename DoFHandler<dim, spacedim>::level_cell_iterator
     end = dof_handler.end(level);
 
     const types::global_dof_index result =
-      compute_block_wise<dim, dim, typename DoFHandler<dim>::level_cell_iterator,
-      typename DoFHandler<dim>::level_cell_iterator>(
-        renumbering, start, end);
+      compute_block_wise<dim, spacedim, typename DoFHandler<dim, spacedim>::level_cell_iterator,
+      typename DoFHandler<dim, spacedim>::level_cell_iterator>(
+        renumbering, start, end, true);
 
     if (result == 0) return;
 
@@ -927,7 +927,8 @@ namespace DoFRenumbering
   types::global_dof_index
   compute_block_wise (std::vector<types::global_dof_index> &new_indices,
                       const ITERATOR    &start,
-                      const ENDITERATOR &end)
+                      const ENDITERATOR &end,
+                      bool is_level_operation)
   {
     const hp::FECollection<dim,spacedim>
     fe_collection (start->get_dof_handler().get_fe ());
@@ -971,20 +972,35 @@ namespace DoFRenumbering
     std::vector<std::vector<types::global_dof_index> >
     block_to_dof_map (fe_collection.n_blocks());
     for (ITERATOR cell=start; cell!=end; ++cell)
-      if (cell->is_locally_owned())
-        {
-          // on each cell: get dof indices
-          // and insert them into the global
-          // list using their component
-          const unsigned int fe_index = cell->active_fe_index();
-          const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
-          local_dof_indices.resize (dofs_per_cell);
-          cell->get_active_or_mg_dof_indices (local_dof_indices);
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            if (start->get_dof_handler().locally_owned_dofs().is_element(local_dof_indices[i]))
-              block_to_dof_map[block_list[fe_index][i]].
-              push_back (local_dof_indices[i]);
-        }
+      {
+        if (is_level_operation)
+          {
+            //we are dealing with mg dofs, skip foreign level cells:
+            if ((start->get_dof_handler().get_tria().locally_owned_subdomain() != numbers::invalid_subdomain_id)
+                && (cell->level_subdomain_id()
+                    != start->get_dof_handler().get_tria().locally_owned_subdomain()))
+              continue;
+          }
+        else
+          {
+            //we are dealing with active dofs, skip the loop if not locally
+            // owned:
+            if (!cell->active() || !cell->is_locally_owned())
+              continue;
+          }
+
+        // on each cell: get dof indices
+        // and insert them into the global
+        // list using their component
+        const unsigned int fe_index = cell->active_fe_index();
+        const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
+        local_dof_indices.resize (dofs_per_cell);
+        cell->get_active_or_mg_dof_indices (local_dof_indices);
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          if (start->get_dof_handler().locally_owned_dofs().is_element(local_dof_indices[i]))
+            block_to_dof_map[block_list[fe_index][i]].
+            push_back (local_dof_indices[i]);
+    }
 
     // now we've got all indices sorted
     // into buckets labeled by their
-- 
2.39.5