From 23ab0eb0c011b9814eda63eed7f2ae95cf3914f3 Mon Sep 17 00:00:00 2001
From: David Wells <wellsd2@rpi.edu>
Date: Wed, 7 Dec 2016 10:01:57 -0500
Subject: [PATCH] Cache end_active() when assigning hp dofs.

This fixes a performance bug where distributing dofs on some grids would take an
amount of time quadratic in the number of dofs when it should be linear.
---
 doc/news/changes.h       |  6 ++++++
 source/hp/dof_handler.cc | 44 ++++++++++++++++++++++------------------
 2 files changed, 30 insertions(+), 20 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 116b67a09e..9faa80be66 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -235,6 +235,12 @@ inconvenience this causes.
 
 <ol>
 
+<li> Improved: The run time for the method hp::DoFHandler::distribute_dofs was
+     quadratic in the total number of dofs for some grids. This has been fixed.
+     <br>
+     (David Wells, 2016/12/07)
+</li>
+
 <li> New: The inverse of a rank-2 SymmetricTensor can now be directly computed
      with SymmetricTensor::invert() instead of having to use the
      Tensor::invert() function.
diff --git a/source/hp/dof_handler.cc b/source/hp/dof_handler.cc
index c0bcf19470..18d7c11403 100644
--- a/source/hp/dof_handler.cc
+++ b/source/hp/dof_handler.cc
@@ -607,9 +607,10 @@ namespace internal
 
               types::global_dof_index next_free_dof = 0;
               types::global_dof_index cache_size = 0;
-              for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
-                   cell=dof_handler.begin_active(level);
-                   cell!=dof_handler.end_active(level); ++cell)
+              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
+              cell=dof_handler.begin_active(level),
+              endc=dof_handler.end_active(level);
+              for (; cell!=endc; ++cell)
                 if (!cell->has_children())
                   {
                     dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
@@ -638,9 +639,10 @@ namespace internal
           for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
             {
               types::global_dof_index counter = 0;
-              for (typename HpDoFHandler<dim,spacedim>::cell_iterator
-                   cell=dof_handler.begin_active(level);
-                   cell!=dof_handler.end_active(level); ++cell)
+              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
+              cell=dof_handler.begin_active(level),
+              endc=dof_handler.end_active(level);
+              for (; cell!=endc; ++cell)
                 if (!cell->has_children())
                   counter += cell->get_fe().dofs_per_line;
 
@@ -705,7 +707,6 @@ namespace internal
             dof_handler.faces = new internal::hp::DoFIndicesOnFaces<2>;
           }
 
-
           // QUAD (CELL) DOFs
 
           // count how much space we need
@@ -735,9 +736,10 @@ namespace internal
 
               types::global_dof_index next_free_dof = 0;
               types::global_dof_index cache_size = 0;
-              for (typename HpDoFHandler<dim, spacedim>::active_cell_iterator
-                   cell=dof_handler.begin_active(level);
-                   cell!=dof_handler.end_active(level); ++cell)
+              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
+              cell=dof_handler.begin_active(level),
+              endc=dof_handler.end_active(level);
+              for (; cell!=endc; ++cell)
                 if (!cell->has_children())
                   {
                     dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
@@ -766,9 +768,10 @@ namespace internal
           for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
             {
               types::global_dof_index counter = 0;
-              for (typename HpDoFHandler<dim,spacedim>::cell_iterator
-                   cell=dof_handler.begin_active(level);
-                   cell!=dof_handler.end_active(level); ++cell)
+              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
+              cell=dof_handler.begin_active(level),
+              endc=dof_handler.end_active(level);
+              for (; cell!=endc; ++cell)
                 if (!cell->has_children())
                   counter += cell->get_fe().dofs_per_quad;
 
@@ -1081,7 +1084,6 @@ namespace internal
             dof_handler.faces = new internal::hp::DoFIndicesOnFaces<3>;
           }
 
-
           // HEX (CELL) DOFs
 
           // count how much space we need
@@ -1111,9 +1113,10 @@ namespace internal
 
               types::global_dof_index next_free_dof = 0;
               types::global_dof_index cache_size = 0;
-              for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
-                   cell=dof_handler.begin_active(level);
-                   cell!=dof_handler.end_active(level); ++cell)
+              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
+              cell=dof_handler.begin_active(level),
+              endc=dof_handler.end_active(level);
+              for (; cell!=endc; ++cell)
                 if (!cell->has_children())
                   {
                     dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
@@ -1142,9 +1145,10 @@ namespace internal
           for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
             {
               types::global_dof_index counter = 0;
-              for (typename HpDoFHandler<dim,spacedim>::cell_iterator
-                   cell=dof_handler.begin_active(level);
-                   cell!=dof_handler.end_active(level); ++cell)
+              typename HpDoFHandler<dim,spacedim>::active_cell_iterator
+              cell=dof_handler.begin_active(level),
+              endc=dof_handler.end_active(level);
+              for (; cell!=endc; ++cell)
                 if (!cell->has_children())
                   counter += cell->get_fe().dofs_per_hex;
 
-- 
2.39.5