From 5ba5f153f5cc7de0341cd7a606fbeda9f0e9b8e7 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 7 Aug 2017 13:51:54 -0600
Subject: [PATCH] Exchange active_fe_index between processors.

Every processor should only set active_fe_index on their locally owned
cells, and then get those on ghost cells through parallel exchange.
---
 source/dofs/dof_handler_policy.cc |  1 +
 source/hp/dof_handler.cc          | 62 +++++++++++++++++++++++++++++++
 2 files changed, 63 insertions(+)

diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc
index 811208a226..e270f834dd 100644
--- a/source/dofs/dof_handler_policy.cc
+++ b/source/dofs/dof_handler_policy.cc
@@ -3513,6 +3513,7 @@ namespace internal
         const unsigned int
         n_cpus = Utilities::MPI::n_mpi_processes (triangulation->get_communicator());
 
+
         /*
            The following algorithm has a number of stages that are all documented
            in the paper that describes the parallel::distributed functionality:
diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc
index 450f8a62dc..a1d12159a4 100644
--- a/source/hp/dof_handler.cc
+++ b/source/hp/dof_handler.cc
@@ -1399,6 +1399,64 @@ namespace hp
 
 
 
+  namespace
+  {
+    /**
+     * Given a hp::DoFHandler object, make sure that the active_fe_indices that
+     * a user has set for locally owned cells are communicated to all ghost
+     * cells as well.
+     */
+    template <int dim, int spacedim>
+    void communicate_active_fe_indices (hp::DoFHandler<dim,spacedim> &dof_handler)
+    {
+      if (const parallel::shared::Triangulation<dim, spacedim> *tr =
+            dynamic_cast<const parallel::shared::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation()))
+        {
+          // we have a shared triangulation. in this case, every processor
+          // knows about all cells, but every processor only has knowledge
+          // about the active_fe_index on the cells it owns.
+          //
+          // we can create a complete set of active_fe_indices by letting
+          // every processor create a vector of indices for all cells,
+          // filling only those on the cells it owns and setting the indices
+          // on the other cells to zero. then we add all of these vectors
+          // up, and because every vector entry has exactly one processor
+          // that owns it, the sum is correct
+          std::vector<unsigned int> active_fe_indices (tr->n_active_cells(),
+                                                       (unsigned int)0);
+          for (const auto &cell : dof_handler.active_cell_iterators())
+            if (cell->is_locally_owned())
+              active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
+
+          Utilities::MPI::sum (active_fe_indices,
+                               tr->get_communicator(),
+                               active_fe_indices);
+
+          // now go back and fill the active_fe_index on ghost cells
+          for (auto cell : dof_handler.active_cell_iterators())
+            if (cell->is_ghost())
+              cell->set_active_fe_index(active_fe_indices[cell->active_cell_index()]);
+        }
+      else if (const parallel::distributed::Triangulation<dim, spacedim> *tr =
+                 dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation()))
+        {
+          Assert (false, ExcNotImplemented());
+        }
+      else
+        {
+          // a sequential triangulation. there is nothing we need to do here
+          Assert ((dynamic_cast<const parallel::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation())
+                   == nullptr),
+                  ExcInternalError());
+        }
+    }
+  }
+
+
+
+
+
+
   template <int dim, int spacedim>
   void DoFHandler<dim,spacedim>::distribute_dofs (const hp::FECollection<dim,spacedim> &ff)
   {
@@ -1406,6 +1464,10 @@ namespace hp
 
     finite_elements = &ff;
 
+    // at the beginning, make sure every processor knows the
+    // active_fe_indices on both its own cells and all ghost cells
+    communicate_active_fe_indices (*this);
+
     // This call ensures that the active_fe_indices vectors are
     // initialized correctly.
     create_active_fe_table ();
-- 
2.39.5