Move the main part of distributing hp DoFs to the policy.

diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc
index a298d25c1e464598e061def3ed63000379b71a00..b18b4c92c81340aae25cf8e9f9314e9655509019 100644 (file)
--- a/source/dofs/dof_handler_policy.cc
+++ b/source/dofs/dof_handler_policy.cc
@@ -289,14 +289,41 @@ namespace internal
         template <int spacedim>
         static
         types::global_dof_index
-        distribute_dofs_on_cell (const hp::DoFHandler<1,spacedim>                                &dof_handler,
+        distribute_dofs_on_cell (const hp::DoFHandler<1,spacedim> &,
                                  const typename hp::DoFHandler<1,spacedim>::active_cell_iterator &cell,
                                  types::global_dof_index                                          next_free_dof)
         {
-          (void)dof_handler;
-          (void)cell;
-          (void)next_free_dof;
-          return 0;
+          const unsigned int dim = 1;
+
+          const FiniteElement<dim,spacedim> &fe       = cell->get_fe();
+          const unsigned int                 fe_index = cell->active_fe_index ();
+
+          // number dofs on vertices. to do so, check whether dofs for
+          // this vertex have been distributed and for the present fe
+          // (only check the first dof), and if this isn't the case
+          // distribute new ones there
+          if (fe.dofs_per_vertex > 0)
+            for (unsigned int vertex=0; vertex<GeometryInfo<1>::vertices_per_cell; ++vertex)
+              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
+                  numbers::invalid_dof_index)
+                for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
+                  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
+
+          // finally for the line. this one shouldn't be numbered yet
+          if (fe.dofs_per_line > 0)
+            {
+              Assert ((cell->dof_index(0, fe_index) ==
+                       numbers::invalid_dof_index),
+                      ExcInternalError());
+
+              for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
+                cell->set_dof_index (d, next_free_dof, fe_index);
+            }
+
+          // note that this cell has been processed
+          cell->set_user_flag ();
+
+          return next_free_dof;
         }
 
 
@@ -304,14 +331,57 @@ namespace internal
         template <int spacedim>
         static
         types::global_dof_index
-        distribute_dofs_on_cell (const hp::DoFHandler<2,spacedim>                                &dof_handler,
+        distribute_dofs_on_cell (const hp::DoFHandler<2,spacedim> &,
                                  const typename hp::DoFHandler<2,spacedim>::active_cell_iterator &cell,
                                  types::global_dof_index                                          next_free_dof)
         {
-          (void)dof_handler;
-          (void)cell;
-          (void)next_free_dof;
-          return 0;
+          const unsigned int dim = 2;
+
+          const FiniteElement<dim,spacedim> &fe       = cell->get_fe();
+          const unsigned int                 fe_index = cell->active_fe_index ();
+
+          // number dofs on vertices. to do so, check whether dofs for
+          // this vertex have been distributed and for the present fe
+          // (only check the first dof), and if this isn't the case
+          // distribute new ones there
+          if (fe.dofs_per_vertex > 0)
+            for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
+              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
+                  numbers::invalid_dof_index)
+                for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
+                  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
+
+          // next the sides. do the same as above: check whether the
+          // line is already numbered for the present fe_index, and if
+          // not do it
+          if (fe.dofs_per_line > 0)
+            for (unsigned int l=0; l<GeometryInfo<2>::lines_per_cell; ++l)
+              {
+                typename hp::DoFHandler<dim,spacedim>::line_iterator
+                line = cell->line(l);
+
+                if (line->dof_index(0,fe_index) ==
+                    numbers::invalid_dof_index)
+                  for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
+                    line->set_dof_index (d, next_free_dof, fe_index);
+              }
+
+
+          // finally for the quad. this one shouldn't be numbered yet
+          if (fe.dofs_per_quad > 0)
+            {
+              Assert ((cell->dof_index(0, fe_index) ==
+                       numbers::invalid_dof_index),
+                      ExcInternalError());
+
+              for (unsigned int d=0; d<fe.dofs_per_quad; ++d, ++next_free_dof)
+                cell->set_dof_index (d, next_free_dof, fe_index);
+            }
+
+          // note that this cell has been processed
+          cell->set_user_flag ();
+
+          return next_free_dof;
         }
 
 
@@ -319,14 +389,70 @@ namespace internal
         template <int spacedim>
         static
         types::global_dof_index
-        distribute_dofs_on_cell (const hp::DoFHandler<3,spacedim>                                &dof_handler,
+        distribute_dofs_on_cell (const hp::DoFHandler<3,spacedim> &,
                                  const typename hp::DoFHandler<3,spacedim>::active_cell_iterator &cell,
                                  types::global_dof_index                                          next_free_dof)
         {
-          (void)dof_handler;
-          (void)cell;
-          (void)next_free_dof;
-          return 0;
+          const unsigned int dim = 3;
+
+          const FiniteElement<dim,spacedim> &fe       = cell->get_fe();
+          const unsigned int                 fe_index = cell->active_fe_index ();
+
+          // number dofs on vertices. to do so, check whether dofs for
+          // this vertex have been distributed and for the present fe
+          // (only check the first dof), and if this isn't the case
+          // distribute new ones there
+          if (fe.dofs_per_vertex > 0)
+            for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
+              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
+                  numbers::invalid_dof_index)
+                for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
+                  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
+
+          // next the four lines. do the same as above: check whether
+          // the line is already numbered for the present fe_index,
+          // and if not do it
+          if (fe.dofs_per_line > 0)
+            for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
+              {
+                typename hp::DoFHandler<dim,spacedim>::line_iterator
+                line = cell->line(l);
+
+                if (line->dof_index(0,fe_index) ==
+                    numbers::invalid_dof_index)
+                  for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
+                    line->set_dof_index (d, next_free_dof, fe_index);
+              }
+
+          // same for quads
+          if (fe.dofs_per_quad > 0)
+            for (unsigned int q=0; q<GeometryInfo<3>::quads_per_cell; ++q)
+              {
+                typename hp::DoFHandler<dim,spacedim>::quad_iterator
+                quad = cell->quad(q);
+
+                if (quad->dof_index(0,fe_index) ==
+                    numbers::invalid_dof_index)
+                  for (unsigned int d=0; d<fe.dofs_per_quad; ++d, ++next_free_dof)
+                    quad->set_dof_index (d, next_free_dof, fe_index);
+              }
+
+
+          // finally for the hex. this one shouldn't be numbered yet
+          if (fe.dofs_per_hex > 0)
+            {
+              Assert ((cell->dof_index(0, fe_index) ==
+                       numbers::invalid_dof_index),
+                      ExcInternalError());
+
+              for (unsigned int d=0; d<fe.dofs_per_hex; ++d, ++next_free_dof)
+                cell->set_dof_index (d, next_free_dof, fe_index);
+            }
+
+          // note that this cell has been processed
+          cell->set_user_flag ();
+
+          return next_free_dof;
         }
 
 

diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc
index 54893a95586d69fe388b7eff6bc39aa8b13a943a..ff56aeff390f7d8d3039ed35186e7fdcf547b111 100644 (file)
--- a/source/hp/dof_handler.cc
+++ b/source/hp/dof_handler.cc
@@ -328,204 +328,6 @@ namespace internal
 
 
 
-        /**
-         * Distribute dofs on the given cell,
-         * with new dofs starting with index
-         * @p next_free_dof. Return the next
-         * unused index number. The finite
-         * element used is the one given to
-         * @p distribute_dofs, which is copied
-         * to @p selected_fe.
-         *
-         * This function is excluded from the
-         * @p distribute_dofs function since
-         * it can not be implemented dimension
-         * independent.
-         */
-        template <int spacedim>
-        static
-        types::global_dof_index
-        distribute_dofs_on_cell (const typename dealii::hp::DoFHandler<1,spacedim>::active_cell_iterator &cell,
-                                 types::global_dof_index                                                 next_free_dof)
-        {
-          const unsigned int dim = 1;
-
-          const FiniteElement<dim,spacedim> &fe       = cell->get_fe();
-          const unsigned int                 fe_index = cell->active_fe_index ();
-
-          // number dofs on vertices. to do
-          // so, check whether dofs for
-          // this vertex have been
-          // distributed and for the
-          // present fe (only check the
-          // first dof), and if this isn't
-          // the case distribute new ones
-          // there
-          if (fe.dofs_per_vertex > 0)
-            for (unsigned int vertex=0; vertex<GeometryInfo<1>::vertices_per_cell; ++vertex)
-              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
-                  numbers::invalid_dof_index)
-                for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
-                  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
-
-          // finally for the line. this one
-          // shouldn't be numbered yet
-          if (fe.dofs_per_line > 0)
-            {
-              Assert ((cell->dof_index(0, fe_index) ==
-                       numbers::invalid_dof_index),
-                      ExcInternalError());
-
-              for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
-                cell->set_dof_index (d, next_free_dof, fe_index);
-            }
-
-          // note that this cell has been processed
-          cell->set_user_flag ();
-
-          return next_free_dof;
-        }
-
-
-        template <int spacedim>
-        static
-        types::global_dof_index
-        distribute_dofs_on_cell (const typename dealii::hp::DoFHandler<2,spacedim>::active_cell_iterator &cell,
-                                 types::global_dof_index                                                 next_free_dof)
-        {
-          const unsigned int dim = 2;
-
-          const FiniteElement<dim,spacedim> &fe       = cell->get_fe();
-          const unsigned int                 fe_index = cell->active_fe_index ();
-
-          // number dofs on vertices. to do
-          // so, check whether dofs for
-          // this vertex have been
-          // distributed and for the
-          // present fe (only check the
-          // first dof), and if this isn't
-          // the case distribute new ones
-          // there
-          if (fe.dofs_per_vertex > 0)
-            for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
-              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
-                  numbers::invalid_dof_index)
-                for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
-                  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
-
-          // next the sides. do the
-          // same as above: check whether
-          // the line is already numbered
-          // for the present fe_index, and
-          // if not do it
-          if (fe.dofs_per_line > 0)
-            for (unsigned int l=0; l<GeometryInfo<2>::lines_per_cell; ++l)
-              {
-                typename HpDoFHandler<dim,spacedim>::line_iterator
-                line = cell->line(l);
-
-                if (line->dof_index(0,fe_index) ==
-                    numbers::invalid_dof_index)
-                  for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
-                    line->set_dof_index (d, next_free_dof, fe_index);
-              }
-
-
-          // finally for the quad. this one
-          // shouldn't be numbered yet
-          if (fe.dofs_per_quad > 0)
-            {
-              Assert ((cell->dof_index(0, fe_index) ==
-                       numbers::invalid_dof_index),
-                      ExcInternalError());
-
-              for (unsigned int d=0; d<fe.dofs_per_quad; ++d, ++next_free_dof)
-                cell->set_dof_index (d, next_free_dof, fe_index);
-            }
-
-          // note that this cell has been processed
-          cell->set_user_flag ();
-
-          return next_free_dof;
-        }
-
-
-        template <int spacedim>
-        static
-        types::global_dof_index
-        distribute_dofs_on_cell (const typename dealii::hp::DoFHandler<3,spacedim>::active_cell_iterator &cell,
-                                 types::global_dof_index                                                 next_free_dof)
-        {
-          const unsigned int dim = 3;
-
-          const FiniteElement<dim,spacedim> &fe       = cell->get_fe();
-          const unsigned int                 fe_index = cell->active_fe_index ();
-
-          // number dofs on vertices. to do
-          // so, check whether dofs for
-          // this vertex have been
-          // distributed and for the
-          // present fe (only check the
-          // first dof), and if this isn't
-          // the case distribute new ones
-          // there
-          if (fe.dofs_per_vertex > 0)
-            for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
-              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
-                  numbers::invalid_dof_index)
-                for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
-                  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
-
-          // next the four lines. do the
-          // same as above: check whether
-          // the line is already numbered
-          // for the present fe_index, and
-          // if not do it
-          if (fe.dofs_per_line > 0)
-            for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
-              {
-                typename HpDoFHandler<dim,spacedim>::line_iterator
-                line = cell->line(l);
-
-                if (line->dof_index(0,fe_index) ==
-                    numbers::invalid_dof_index)
-                  for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
-                    line->set_dof_index (d, next_free_dof, fe_index);
-              }
-
-          // same for quads
-          if (fe.dofs_per_quad > 0)
-            for (unsigned int q=0; q<GeometryInfo<3>::quads_per_cell; ++q)
-              {
-                typename HpDoFHandler<dim,spacedim>::quad_iterator
-                quad = cell->quad(q);
-
-                if (quad->dof_index(0,fe_index) ==
-                    numbers::invalid_dof_index)
-                  for (unsigned int d=0; d<fe.dofs_per_quad; ++d, ++next_free_dof)
-                    quad->set_dof_index (d, next_free_dof, fe_index);
-              }
-
-
-          // finally for the hex. this one
-          // shouldn't be numbered yet
-          if (fe.dofs_per_hex > 0)
-            {
-              Assert ((cell->dof_index(0, fe_index) ==
-                       numbers::invalid_dof_index),
-                      ExcInternalError());
-
-              for (unsigned int d=0; d<fe.dofs_per_hex; ++d, ++next_free_dof)
-                cell->set_dof_index (d, next_free_dof, fe_index);
-            }
-
-          // note that this cell has been processed
-          cell->set_user_flag ();
-
-          return next_free_dof;
-        }
-
-
         /**
          * Reserve enough space in the
          * <tt>levels[]</tt> objects to store the
@@ -2626,18 +2428,7 @@ namespace hp
 
     // Step 1: distribute DoFs on all
     // active entities
-    {
-      types::global_dof_index next_free_dof = 0;
-      active_cell_iterator cell = begin_active(),
-                           endc = end();
-
-      for (; cell != endc; ++cell)
-        next_free_dof
-          = dealii::internal::hp::DoFHandler::Implementation::distribute_dofs_on_cell<spacedim> (cell,
-              next_free_dof);
-
-      number_cache.n_global_dofs = next_free_dof;
-    }
+    number_cache = policy->distribute_dofs ();
 
 
     /////////////////////////////////