From 61336fa989f0ca367a635ef71dfccac6307d20b4 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Sat, 1 Jul 2017 01:20:12 -0600
Subject: [PATCH] Move the DoF renumbering functionality.

Specifically, move it out of the DoFHandler class and into the policy class
where it can use shared infrastructure.
---
 include/deal.II/hp/dof_handler.h  |  23 --
 source/dofs/dof_handler.cc        |   4 +-
 source/dofs/dof_handler_policy.cc | 268 +++++++++++++++++++++--
 source/hp/dof_handler.cc          | 340 +++---------------------------
 4 files changed, 279 insertions(+), 356 deletions(-)

diff --git a/include/deal.II/hp/dof_handler.h b/include/deal.II/hp/dof_handler.h
index 3ab91b2bc3..1e72595cee 100644
--- a/include/deal.II/hp/dof_handler.h
+++ b/include/deal.II/hp/dof_handler.h
@@ -810,29 +810,6 @@ namespace hp
     void
     compute_quad_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const;
 
-    /**
-     * Renumber the objects with the given and all lower structural
-     * dimensions, i.e. renumber vertices by giving a template argument of
-     * zero to the int2type argument, lines and vertices with one, etc.
-     *
-     * Note that in contrast to the public renumber_dofs() function, these
-     * internal functions do not ensure that the new DoFs are contiguously
-     * numbered. The function may therefore also be used to assign different
-     * DoFs the same number, for example to unify hp DoFs corresponding to
-     * different finite elements but co-located on the same entity.
-     */
-    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                                 dealii::internal::int2type<0>);
-
-    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                                 dealii::internal::int2type<1>);
-
-    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                                 dealii::internal::int2type<2>);
-
-    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                                 dealii::internal::int2type<3>);
-
     /**
      * Space to store the DoF numbers for the different levels. Analogous to
      * the <tt>levels[]</tt> tree of the Triangulation objects.
diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc
index b54fab23f2..11376c30be 100644
--- a/source/dofs/dof_handler.cc
+++ b/source/dofs/dof_handler.cc
@@ -653,9 +653,7 @@ DoFHandler<dim,spacedim>::DoFHandler (const Triangulation<dim,spacedim> &tria)
   faces(nullptr),
   mg_faces (nullptr)
 {
-  // decide whether we need a
-  // sequential or a parallel
-  // distributed policy
+  // decide whether we need a sequential or a parallel distributed policy
   if (dynamic_cast<const parallel::shared::Triangulation< dim, spacedim>*>
       (&tria)
       != nullptr)
diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc
index 5be7b4676f..47d9205bb2 100644
--- a/source/dofs/dof_handler_policy.cc
+++ b/source/dofs/dof_handler_policy.cc
@@ -761,6 +761,255 @@ namespace internal
 
 
 
+        template <int dim, int spacedim>
+        static
+        void
+        renumber_vertex_dofs (const std::vector<types::global_dof_index> &new_numbers,
+                              const IndexSet                             &indices,
+                              hp::DoFHandler<dim,spacedim>               &dof_handler,
+                              const bool                                  check_validity)
+        {
+          for (unsigned int vertex_index=0; vertex_index<dof_handler.get_triangulation().n_vertices();
+               ++vertex_index)
+            {
+              const unsigned int n_active_fe_indices
+                = dealii::internal::DoFAccessor::Implementation::
+                  n_active_vertex_fe_indices (dof_handler, vertex_index);
+
+              for (unsigned int f=0; f<n_active_fe_indices; ++f)
+                {
+                  const unsigned int fe_index
+                    = dealii::internal::DoFAccessor::Implementation::
+                      nth_active_vertex_fe_index (dof_handler, vertex_index, f);
+
+                  for (unsigned int d=0; d<dof_handler.get_fe()[fe_index].dofs_per_vertex; ++d)
+                    {
+                      const types::global_dof_index old_dof_index
+                        = dealii::internal::DoFAccessor::Implementation::
+                          get_vertex_dof_index(dof_handler,
+                                               vertex_index,
+                                               fe_index,
+                                               d);
+                      if (old_dof_index != numbers::invalid_dof_index)
+                        dealii::internal::DoFAccessor::Implementation::
+                        set_vertex_dof_index (dof_handler,
+                                              vertex_index,
+                                              fe_index,
+                                              d,
+                                              (indices.size() == 0)?
+                                              (new_numbers[old_dof_index]) :
+                                              (new_numbers[indices.index_within_set(old_dof_index)]));
+                      else if (check_validity)
+                        // if index is invalid_dof_index: check if this one
+                        // really is unused
+                        Assert (dof_handler.get_triangulation()
+                                .vertex_used(vertex_index)
+                                == false,
+                                ExcInternalError ());
+                    }
+                }
+            }
+        }
+
+
+
+        template <int dim, int spacedim>
+        static
+        void
+        renumber_cell_dofs (const std::vector<types::global_dof_index> &new_numbers,
+                            const IndexSet                             &indices,
+                            hp::DoFHandler<dim,spacedim>               &dof_handler)
+        {
+          for (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
+               cell = dof_handler.begin_active();
+               cell!=dof_handler.end(); ++cell)
+            {
+              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);
+                }
+            }
+        }
+
+
+
+        template <int spacedim>
+        static
+        void
+        renumber_face_dofs (const std::vector<types::global_dof_index> &/*new_numbers*/,
+                            const IndexSet                             &/*indices*/,
+                            hp::DoFHandler<1,spacedim>                 &/*dof_handler*/)
+        {
+          // nothing to do in 1d since there are no separate faces -- we've already
+          // taken care of this when dealing with the vertices
+        }
+
+
+
+        template <int spacedim>
+        static
+        void
+        renumber_face_dofs (const std::vector<types::global_dof_index> &new_numbers,
+                            const IndexSet                             &indices,
+                            hp::DoFHandler<2,spacedim>                 &dof_handler)
+        {
+          const unsigned int dim = 2;
+
+          // deal with DoFs on lines
+          {
+            // save user flags on lines so we can use them to mark lines
+            // we've already treated
+            std::vector<bool> saved_line_user_flags;
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .save_user_flags_line (saved_line_user_flags);
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .clear_user_flags_line ();
+
+            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 ();
+
+                    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
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .load_user_flags_line (saved_line_user_flags);
+          }
+        }
+
+
+
+        template <int spacedim>
+        static
+        void
+        renumber_face_dofs (const std::vector<types::global_dof_index> &new_numbers,
+                            const IndexSet                             &indices,
+                            hp::DoFHandler<3,spacedim>                 &dof_handler)
+        {
+          const unsigned int dim = 3;
+
+          // deal with DoFs on lines
+          {
+            // save user flags on lines so we can use them to mark lines
+            // we've already treated
+            std::vector<bool> saved_line_user_flags;
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .save_user_flags_line (saved_line_user_flags);
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .clear_user_flags_line ();
+
+            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 ();
+
+                    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
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .load_user_flags_line (saved_line_user_flags);
+          }
+
+          // then deal with dofs on quads
+          {
+            std::vector<bool> saved_quad_user_flags;
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .save_user_flags_quad (saved_quad_user_flags);
+            const_cast<dealii::Triangulation<dim,spacedim>&>(dof_handler.get_triangulation())
+            .clear_user_flags_quad ();
+
+            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 ();
+
+                    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())
+            .load_user_flags_quad (saved_quad_user_flags);
+          }
+        }
+
+
+
 
         /**
          * Implementation of DoFHandler::renumber_dofs()
@@ -773,15 +1022,15 @@ namespace internal
          * (The IndexSet argument is not used in 1d because we only need
          * it for parallel meshes and 1d doesn't support that right now.)
          */
-        template <int dim, int spacedim>
+        template <typename DoFHandlerType>
         static
         void
         renumber_dofs (const std::vector<types::global_dof_index> &new_numbers,
                        const IndexSet                             &indices,
-                       DoFHandler<dim,spacedim>                   &dof_handler,
+                       DoFHandlerType                             &dof_handler,
                        const bool                                  check_validity)
         {
-          if (dim == 1)
+          if (DoFHandlerType::dimension == 1)
             Assert (indices == IndexSet(0), ExcNotImplemented());
 
           // renumber DoF indices on vertices, cells, and faces. this
@@ -808,19 +1057,6 @@ namespace internal
 
 
 
-        template <int dim, int spacedim>
-        static
-        void
-        renumber_dofs (const std::vector<types::global_dof_index> &/*new_numbers*/,
-                       const IndexSet                             &/*indices*/,
-                       hp::DoFHandler<dim,spacedim>               &/*dof_handler*/,
-                       const bool                                  /*check_validity*/)
-        {
-          Assert (false, ExcNotImplemented());
-        }
-
-
-
         /* --------------------- renumber_mg_dofs functionality ---------------- */
 
         /**
diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc
index 55b874de1e..4ec421b208 100644
--- a/source/hp/dof_handler.cc
+++ b/source/hp/dof_handler.cc
@@ -1672,6 +1672,19 @@ namespace hp
             ExcMessage ("The given triangulation is parallel distributed but "
                         "this class does not currently support this."));
 
+    // decide whether we need a sequential or a parallel distributed policy
+//    if (dynamic_cast<const parallel::shared::Triangulation< dim, spacedim>*>
+//        (&tria)
+//        != nullptr)
+//      policy.reset (new internal::DoFHandler::Policy::ParallelShared<dim,spacedim>(*this));
+//    else if (dynamic_cast<const parallel::distributed::Triangulation< dim, spacedim >*>
+//             (&tria)
+//             == nullptr)
+    policy.reset (new internal::DoFHandler::Policy::Sequential<DoFHandler<dim,spacedim> >(*this));
+//    else
+//      policy.reset (new internal::DoFHandler::Policy::ParallelDistributed<dim,spacedim>(*this));
+
+
     create_active_fe_table ();
 
     tria_listeners.push_back
@@ -2676,47 +2689,21 @@ namespace hp
                 ExcInternalError());
       }
 
-    // finally, do the renumbering
-    // and set the number of actually
+    // finally, do the renumbering and set the number of actually
     // used dof indices
-    renumber_dofs_internal (new_dof_indices, dealii::internal::int2type<dim>());
+    policy->renumber_dofs (new_dof_indices);
 
-    // now set the elements of the
-    // number cache appropriately
-    number_cache.n_global_dofs        = next_free_dof;
-    number_cache.n_locally_owned_dofs = number_cache.n_global_dofs;
+    // now set the elements of the number cache appropriately
+    number_cache = dealii::internal::DoFHandler::NumberCache (next_free_dof);
 
-    if (dynamic_cast<const parallel::shared::Triangulation< dim, spacedim >*>
-        (&this->get_triangulation())
-        == nullptr)
-      {
-        number_cache.locally_owned_dofs
-          = IndexSet (number_cache.n_global_dofs);
-        number_cache.locally_owned_dofs.add_range (0,
-                                                   number_cache.n_global_dofs);
-        Assert (number_cache.n_global_dofs < std::numeric_limits<unsigned int>::max (),
-                ExcMessage ("Global number of degrees of freedom is too large."));
-        number_cache.n_locally_owned_dofs_per_processor
-          = std::vector<types::global_dof_index> (1,
-                                                  (types::global_dof_index) number_cache.n_global_dofs);
-      }
-    else
-      {
-        AssertThrow(false, ExcNotImplemented() );
-        //number_cache.locally_owned_dofs = dealii::DoFTools::locally_owned_dofs_with_subdomain(this,tria->locally_owned_subdomain() );
-        //TODO: update n_locally_owned_dofs_per_processor as well
-      }
-
-    number_cache.locally_owned_dofs_per_processor
-      = std::vector<IndexSet> (1,
-                               number_cache.locally_owned_dofs);
-
-    // update the cache used for cell dof indices and compress the data on the levels. do
-    // the latter on separate tasks to gain parallelism, starting with the highest
-    // level (there is most to do there, so start it first)
-    for (active_cell_iterator cell = begin_active();
-         cell != end(); ++cell)
-      cell->update_cell_dof_indices_cache ();
+    Assert ((dynamic_cast<const parallel::shared::Triangulation< dim, spacedim >*>
+             (&this->get_triangulation())
+             == nullptr),
+            ExcNotImplemented());
+    Assert ((dynamic_cast<const parallel::distributed::Triangulation< dim, spacedim >*>
+             (&this->get_triangulation())
+             == nullptr),
+            ExcNotImplemented());
 
     {
       Threads::TaskGroup<> tg;
@@ -2774,12 +2761,7 @@ namespace hp
     }
 
     // do the renumbering
-    renumber_dofs_internal (new_numbers, dealii::internal::int2type<dim>());
-
-    // update the cache used for cell dof indices
-    for (active_cell_iterator cell = begin_active();
-         cell != end(); ++cell)
-      cell->update_cell_dof_indices_cache ();
+    number_cache = policy->renumber_dofs(new_numbers);
 
     // now re-compress the dof indices
     {
@@ -2793,276 +2775,6 @@ namespace hp
 
 
 
-  template <int dim, int spacedim>
-  void
-  DoFHandler<dim,spacedim>::
-  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                          dealii::internal::int2type<0>)
-  {
-    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-
-    for (unsigned int vertex_index=0; vertex_index<get_triangulation().n_vertices();
-         ++vertex_index)
-      {
-        const unsigned int n_active_fe_indices
-          = dealii::internal::DoFAccessor::Implementation::
-            n_active_vertex_fe_indices (*this, vertex_index);
-
-        for (unsigned int f=0; f<n_active_fe_indices; ++f)
-          {
-            const unsigned int fe_index
-              = dealii::internal::DoFAccessor::Implementation::
-                nth_active_vertex_fe_index (*this, vertex_index, f);
-
-            for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_vertex; ++d)
-              {
-                const types::global_dof_index vertex_dof_index
-                  = dealii::internal::DoFAccessor::Implementation::
-                    get_vertex_dof_index(*this,
-                                         vertex_index,
-                                         fe_index,
-                                         d);
-                dealii::internal::DoFAccessor::Implementation::
-                set_vertex_dof_index (*this,
-                                      vertex_index,
-                                      fe_index,
-                                      d,
-                                      new_numbers[vertex_dof_index]);
-              }
-          }
-      }
-  }
-
-
-
-  template <int dim, int spacedim>
-  void
-  DoFHandler<dim,spacedim>::
-  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                          dealii::internal::int2type<1>)
-  {
-    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-
-    renumber_dofs_internal (new_numbers, internal::int2type<0>());
-
-    // save user flags on lines so we
-    // can use them to mark lines
-    // we've already treated
-    std::vector<bool> saved_line_user_flags;
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .save_user_flags_line (saved_line_user_flags);
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .clear_user_flags_line ();
-
-    for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
-      for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
-        if (cell->line(l)->user_flag_set() == false)
-          {
-            const line_iterator line = cell->line(l);
-            line->set_user_flag();
-
-            const unsigned int n_active_fe_indices
-              = line->n_active_fe_indices ();
-
-            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<(*finite_elements)[fe_index].dofs_per_line; ++d)
-                  line->set_dof_index (d,
-                                       new_numbers[line->dof_index(d,fe_index)],
-                                       fe_index);
-              }
-          }
-
-    // at the end, restore the user
-    // flags for the lines
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .load_user_flags_line (saved_line_user_flags);
-  }
-
-
-
-//TODO: Merge the following three functions -- they are identical
-  template <>
-  void
-  DoFHandler<2,2>::
-  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                          dealii::internal::int2type<2>)
-  {
-    const unsigned int dim = 2;
-    const unsigned int spacedim = 2;
-
-    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-
-    renumber_dofs_internal (new_numbers, internal::int2type<1>());
-
-    // save user flags on quads so we
-    // can use them to mark quads
-    // we've already treated
-    std::vector<bool> saved_quad_user_flags;
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .save_user_flags_quad (saved_quad_user_flags);
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .clear_user_flags_quad ();
-
-    for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
-      for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
-        if (cell->quad(q)->user_flag_set() == false)
-          {
-            const quad_iterator quad = cell->quad(q);
-            quad->set_user_flag();
-
-            const unsigned int n_active_fe_indices
-              = quad->n_active_fe_indices ();
-
-            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<(*finite_elements)[fe_index].dofs_per_quad; ++d)
-                  quad->set_dof_index (d,
-                                       new_numbers[quad->dof_index(d,fe_index)],
-                                       fe_index);
-              }
-          }
-
-    // at the end, restore the user
-    // flags for the quads
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .load_user_flags_quad (saved_quad_user_flags);
-  }
-
-
-
-  template <>
-  void
-  DoFHandler<2,3>::
-  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                          dealii::internal::int2type<2>)
-  {
-    const unsigned int dim = 2;
-    const unsigned int spacedim = 3;
-
-    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-
-    renumber_dofs_internal (new_numbers, internal::int2type<1>());
-
-    // save user flags on quads so we
-    // can use them to mark quads
-    // we've already treated
-    std::vector<bool> saved_quad_user_flags;
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .save_user_flags_quad (saved_quad_user_flags);
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .clear_user_flags_quad ();
-
-    for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
-      for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
-        if (cell->quad(q)->user_flag_set() == false)
-          {
-            const quad_iterator quad = cell->quad(q);
-            quad->set_user_flag();
-
-            const unsigned int n_active_fe_indices
-              = quad->n_active_fe_indices ();
-
-            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<(*finite_elements)[fe_index].dofs_per_quad; ++d)
-                  quad->set_dof_index (d,
-                                       new_numbers[quad->dof_index(d,fe_index)],
-                                       fe_index);
-              }
-          }
-
-    // at the end, restore the user
-    // flags for the quads
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .load_user_flags_quad (saved_quad_user_flags);
-  }
-
-
-  template <>
-  void
-  DoFHandler<3,3>::
-  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                          dealii::internal::int2type<2>)
-  {
-    const unsigned int dim = 3;
-    const unsigned int spacedim = 3;
-
-    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-
-    renumber_dofs_internal (new_numbers, internal::int2type<1>());
-
-    // save user flags on quads so we
-    // can use them to mark quads
-    // we've already treated
-    std::vector<bool> saved_quad_user_flags;
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .save_user_flags_quad (saved_quad_user_flags);
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .clear_user_flags_quad ();
-
-    for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
-      for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
-        if (cell->quad(q)->user_flag_set() == false)
-          {
-            const quad_iterator quad = cell->quad(q);
-            quad->set_user_flag();
-
-            const unsigned int n_active_fe_indices
-              = quad->n_active_fe_indices ();
-
-            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<(*finite_elements)[fe_index].dofs_per_quad; ++d)
-                  quad->set_dof_index (d,
-                                       new_numbers[quad->dof_index(d,fe_index)],
-                                       fe_index);
-              }
-          }
-
-    // at the end, restore the user
-    // flags for the quads
-    const_cast<dealii::Triangulation<dim,spacedim>&>(*tria)
-    .load_user_flags_quad (saved_quad_user_flags);
-  }
-
-
-  template <>
-  void
-  DoFHandler<3,3>::
-  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
-                          dealii::internal::int2type<3>)
-  {
-    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-
-    renumber_dofs_internal (new_numbers, internal::int2type<2>());
-
-    // we're in 3d, so hexes are also cells. by design, we only visit each cell once
-    for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
-      {
-        const unsigned int fe_index = cell->active_fe_index ();
-
-        for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_hex; ++d)
-          cell->set_dof_index (d,
-                               new_numbers[cell->dof_index(d,fe_index)],
-                               fe_index);
-      }
-  }
-
-
-
   template <int dim, int spacedim>
   unsigned int
   DoFHandler<dim, spacedim>::max_couplings_between_dofs () const
-- 
2.39.5