Do not set active_fe_indices on non-locally owned cells.

On the other hand, this piece of code tries to do exactly this: after ghost
exchange, we try to set the active_fe_index on a ghost cell to the value
we received from the proper owner of that cell. We work around this by
avoiding cell->set_active_fe_index() on such cells, and manipulate the
underlying data structure directly.

Because these data structures are private, we can no longer do this in a
segregated function in an anonymous namespace, but need to do this in a
function that lives in the 'Implementation' struct.

diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc
index 623742f127519e76b6d85a77b193bff825a20644..a37f1798fb62867d8a5409d9fab1a0ec9469661e 100644 (file)
--- a/source/hp/dof_handler.cc
+++ b/source/hp/dof_handler.cc
@@ -863,6 +863,65 @@ namespace internal
 
           return std::min(max_couplings,dof_handler.n_dofs());
         }
+
+
+        /**
+         * 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>
+        static
+        void communicate_active_fe_indices (dealii::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. we would like to call cell->set_active_fe_index(),
+              // but that function does not allow setting these indices on
+              // non-locally_owned cells. so we have to work around the
+              // issue a little bit by accessing the underlying data
+              // structures directly
+              for (auto cell : dof_handler.active_cell_iterators())
+                if (cell->is_ghost())
+                  dof_handler.levels[cell->level()]->
+                  set_active_fe_index(cell->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());
+            }
+        }
       };
     }
   }
@@ -1199,64 +1258,6 @@ 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)
   {
@@ -1266,7 +1267,7 @@ namespace hp
 
     // 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);
+    dealii::internal::hp::DoFHandler::Implementation::communicate_active_fe_indices (*this);
 
     // If an underlying shared::Tria allows artificial cells,
     // then save the current set of subdomain ids, and set