From bbcebc42821fe98b950e9b4fdd4fa0ba75de9695 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 12 Sep 2017 13:42:59 -0600
Subject: [PATCH] Limit renumbering hp DoFs to non-artificial cells.

---
 source/dofs/dof_handler_policy.cc | 176 +++++++++++++++---------------
 1 file changed, 90 insertions(+), 86 deletions(-)

diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc
index f6cd3dacd2..ff966fe300 100644
--- a/source/dofs/dof_handler_policy.cc
+++ b/source/dofs/dof_handler_policy.cc
@@ -1614,22 +1614,23 @@ namespace internal
                             hp::DoFHandler<dim,spacedim>               &dof_handler)
         {
           for (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
-               cell = dof_handler.begin_active();
+               cell=dof_handler.begin_active();
                cell!=dof_handler.end(); ++cell)
-            {
-              const unsigned int fe_index = cell->active_fe_index ();
+            if (!cell->is_artificial())
+              {
+                const unsigned int fe_index = cell->active_fe_index ();
 
-              for (unsigned int d=0; d<dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>(); ++d)
-                {
-                  const types::global_dof_index old_dof_index = cell->dof_index(d,fe_index);
-                  if (old_dof_index != numbers::invalid_dof_index)
-                    cell->set_dof_index (d,
-                                         (indices.size() == 0)?
-                                         (new_numbers[old_dof_index]) :
-                                         (new_numbers[indices.index_within_set(old_dof_index)]),
-                                         fe_index);
-                }
-            }
+                for (unsigned int d=0; d<dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>(); ++d)
+                  {
+                    const types::global_dof_index old_dof_index = cell->dof_index(d,fe_index);
+                    if (old_dof_index != numbers::invalid_dof_index)
+                      cell->set_dof_index (d,
+                                           (indices.size() == 0)?
+                                           (new_numbers[old_dof_index]) :
+                                           (new_numbers[indices.index_within_set(old_dof_index)]),
+                                           fe_index);
+                  }
+              }
         }
 
 
@@ -1668,32 +1669,33 @@ namespace internal
 
             for (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
                  cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
-              for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
-                if (cell->line(l)->user_flag_set() == false)
-                  {
-                    const typename hp::DoFHandler<dim,spacedim>::line_iterator line = cell->line(l);
-                    line->set_user_flag();
-
-                    const unsigned int n_active_fe_indices
-                      = line->n_active_fe_indices ();
+              if (!cell->is_artificial())
+                for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
+                  if (cell->line(l)->user_flag_set() == false)
+                    {
+                      const typename hp::DoFHandler<dim,spacedim>::line_iterator line = cell->line(l);
+                      line->set_user_flag();
 
-                    for (unsigned int f=0; f<n_active_fe_indices; ++f)
-                      {
-                        const unsigned int fe_index
-                          = line->nth_active_fe_index (f);
+                      const unsigned int n_active_fe_indices
+                        = line->n_active_fe_indices ();
 
-                        for (unsigned int d=0; d<dof_handler.get_fe(fe_index).dofs_per_line; ++d)
-                          {
-                            const types::global_dof_index old_dof_index = line->dof_index(d,fe_index);
-                            if (old_dof_index != numbers::invalid_dof_index)
-                              line->set_dof_index (d,
-                                                   (indices.size() == 0)?
-                                                   (new_numbers[old_dof_index]) :
-                                                   (new_numbers[indices.index_within_set(old_dof_index)]),
-                                                   fe_index);
-                          }
-                      }
-                  }
+                      for (unsigned int f=0; f<n_active_fe_indices; ++f)
+                        {
+                          const unsigned int fe_index
+                            = line->nth_active_fe_index (f);
+
+                          for (unsigned int d=0; d<dof_handler.get_fe(fe_index).dofs_per_line; ++d)
+                            {
+                              const types::global_dof_index old_dof_index = line->dof_index(d,fe_index);
+                              if (old_dof_index != numbers::invalid_dof_index)
+                                line->set_dof_index (d,
+                                                     (indices.size() == 0)?
+                                                     (new_numbers[old_dof_index]) :
+                                                     (new_numbers[indices.index_within_set(old_dof_index)]),
+                                                     fe_index);
+                            }
+                        }
+                    }
 
             // at the end, restore the user
             // flags for the lines
@@ -1725,32 +1727,33 @@ namespace internal
 
             for (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
                  cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
-              for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
-                if (cell->line(l)->user_flag_set() == false)
-                  {
-                    const typename hp::DoFHandler<dim,spacedim>::line_iterator line = cell->line(l);
-                    line->set_user_flag();
-
-                    const unsigned int n_active_fe_indices
-                      = line->n_active_fe_indices ();
+              if (!cell->is_artificial())
+                for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
+                  if (cell->line(l)->user_flag_set() == false)
+                    {
+                      const typename hp::DoFHandler<dim,spacedim>::line_iterator line = cell->line(l);
+                      line->set_user_flag();
 
-                    for (unsigned int f=0; f<n_active_fe_indices; ++f)
-                      {
-                        const unsigned int fe_index
-                          = line->nth_active_fe_index (f);
+                      const unsigned int n_active_fe_indices
+                        = line->n_active_fe_indices ();
 
-                        for (unsigned int d=0; d<dof_handler.get_fe(fe_index).dofs_per_line; ++d)
-                          {
-                            const types::global_dof_index old_dof_index = line->dof_index(d,fe_index);
-                            if (old_dof_index != numbers::invalid_dof_index)
-                              line->set_dof_index (d,
-                                                   (indices.size() == 0)?
-                                                   (new_numbers[old_dof_index]) :
-                                                   (new_numbers[indices.index_within_set(old_dof_index)]),
-                                                   fe_index);
-                          }
-                      }
-                  }
+                      for (unsigned int f=0; f<n_active_fe_indices; ++f)
+                        {
+                          const unsigned int fe_index
+                            = line->nth_active_fe_index (f);
+
+                          for (unsigned int d=0; d<dof_handler.get_fe(fe_index).dofs_per_line; ++d)
+                            {
+                              const types::global_dof_index old_dof_index = line->dof_index(d,fe_index);
+                              if (old_dof_index != numbers::invalid_dof_index)
+                                line->set_dof_index (d,
+                                                     (indices.size() == 0)?
+                                                     (new_numbers[old_dof_index]) :
+                                                     (new_numbers[indices.index_within_set(old_dof_index)]),
+                                                     fe_index);
+                            }
+                        }
+                    }
 
             // at the end, restore the user
             // flags for the lines
@@ -1768,32 +1771,33 @@ namespace internal
 
             for (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
                  cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
-              for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
-                if (cell->quad(q)->user_flag_set() == false)
-                  {
-                    const typename hp::DoFHandler<dim,spacedim>::quad_iterator quad = cell->quad(q);
-                    quad->set_user_flag();
-
-                    const unsigned int n_active_fe_indices
-                      = quad->n_active_fe_indices ();
+              if (!cell->is_artificial())
+                for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
+                  if (cell->quad(q)->user_flag_set() == false)
+                    {
+                      const typename hp::DoFHandler<dim,spacedim>::quad_iterator quad = cell->quad(q);
+                      quad->set_user_flag();
 
-                    for (unsigned int f=0; f<n_active_fe_indices; ++f)
-                      {
-                        const unsigned int fe_index
-                          = quad->nth_active_fe_index (f);
+                      const unsigned int n_active_fe_indices
+                        = quad->n_active_fe_indices ();
 
-                        for (unsigned int d=0; d<dof_handler.get_fe(fe_index).dofs_per_quad; ++d)
-                          {
-                            const types::global_dof_index old_dof_index = quad->dof_index(d,fe_index);
-                            if (old_dof_index != numbers::invalid_dof_index)
-                              quad->set_dof_index (d,
-                                                   (indices.size() == 0)?
-                                                   (new_numbers[old_dof_index]) :
-                                                   (new_numbers[indices.index_within_set(old_dof_index)]),
-                                                   fe_index);
-                          }
-                      }
-                  }
+                      for (unsigned int f=0; f<n_active_fe_indices; ++f)
+                        {
+                          const unsigned int fe_index
+                            = quad->nth_active_fe_index (f);
+
+                          for (unsigned int d=0; d<dof_handler.get_fe(fe_index).dofs_per_quad; ++d)
+                            {
+                              const types::global_dof_index old_dof_index = quad->dof_index(d,fe_index);
+                              if (old_dof_index != numbers::invalid_dof_index)
+                                quad->set_dof_index (d,
+                                                     (indices.size() == 0)?
+                                                     (new_numbers[old_dof_index]) :
+                                                     (new_numbers[indices.index_within_set(old_dof_index)]),
+                                                     fe_index);
+                            }
+                        }
+                    }
 
             // at the end, restore the user flags for the quads
             const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
-- 
2.39.5