From 8c928cfd88fbce76039628162a55a182947a9ef6 Mon Sep 17 00:00:00 2001
From: Marc Fehling <mafehling.git@gmail.com>
Date: Tue, 2 Mar 2021 21:13:25 -0700
Subject: [PATCH] Optimize DoFHandler::prepare_for_coarsening_and_refinement().

---
 source/dofs/dof_handler.cc | 107 +++++++++++++++++--------------------
 1 file changed, 48 insertions(+), 59 deletions(-)

diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc
index 7f9aa9839e..8d30cbbeff 100644
--- a/source/dofs/dof_handler.cc
+++ b/source/dofs/dof_handler.cc
@@ -2677,7 +2677,6 @@ DoFHandler<dim, spacedim>::prepare_coarsening_and_refinement(
   // establish hierarchy
   //
   // - create bimap between hierarchy levels and FE indices
-  // - classify cells based on their level
 
   // map from FE index to level in hierarchy
   // FE indices that are not covered in the hierarchy are not in the map
@@ -2690,77 +2689,67 @@ DoFHandler<dim, spacedim>::prepare_coarsening_and_refinement(
   for (unsigned int i = 0; i < fe_index_for_hierarchy_level.size(); ++i)
     hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[i]] = i;
 
-  // map from level in hierarchy to cells on that particular level
-  // cells that are not covered in the hierarchy are not in the map
-  std::vector<std::list<active_cell_iterator>> cells_on_hierarchy_level(
-    fe_index_for_hierarchy_level.size());
-  if (fe_index_for_hierarchy_level.size() == fe_collection.size())
-    {
-      // all FE indices are part of hierarchy
-      for (const auto &cell : active_cell_iterators())
-        cells_on_hierarchy_level
-          [hierarchy_level_for_fe_index[cell->future_fe_index()]]
-            .push_back(cell);
-    }
-  else
-    {
-      // only some FE indices are part of hierarchy
-      typename decltype(hierarchy_level_for_fe_index)::iterator iter;
-      for (const auto &cell : active_cell_iterators())
-        {
-          iter = hierarchy_level_for_fe_index.find(cell->future_fe_index());
-          if (iter != hierarchy_level_for_fe_index.end())
-            cells_on_hierarchy_level[iter->second].push_back(cell);
-        }
-    }
-
   //
   // limit level difference of neighboring cells
   //
   // - always raise levels to match criterion, never lower them
   //   (hence iterating from highest to lowest level)
   // - update future FE indices
-  // - move cells in level representation correspondigly
-
-  bool changed_fe_indices = false;
-  for (unsigned int level = fe_index_for_hierarchy_level.size() - 1;
-       level > max_difference;
-       --level)
-    for (const auto &cell : cells_on_hierarchy_level[level])
-      for (unsigned int f = 0; f < cell->n_faces(); ++f)
-        if (cell->face(f)->at_boundary() == false)
-          {
-            const auto neighbor = cell->neighbor(f);
 
-            const auto neighbor_fe_and_level =
-              hierarchy_level_for_fe_index.find(neighbor->future_fe_index());
+  bool fe_indices_changed = false;
+  bool fe_indices_changed_in_cycle;
+  do
+    {
+      fe_indices_changed_in_cycle = false;
+
+      for (const auto &cell : active_cell_iterators())
+        {
+          const auto cell_fe_and_level =
+            hierarchy_level_for_fe_index.find(cell->future_fe_index());
 
-            // ignore neighbors that are not part of the hierarchy
-            if (neighbor_fe_and_level == hierarchy_level_for_fe_index.end())
-              continue;
+          // ignore cells that are not part of the hierarchy
+          if (cell_fe_and_level == hierarchy_level_for_fe_index.end())
+            continue;
 
-            const auto neighbor_level = neighbor_fe_and_level->second;
+          const unsigned int cell_level = cell_fe_and_level->second;
 
-            if ((level - max_difference) > neighbor_level)
+          // ignore lowest levels of the hierarchy that always fulfill the
+          // max_difference criterion
+          if (cell_level <= max_difference)
+            continue;
+
+          for (unsigned int f = 0; f < cell->n_faces(); ++f)
+            if (cell->face(f)->at_boundary() == false)
               {
-                // remove neighbor from old level representation
-                auto &     old_cells = cells_on_hierarchy_level[neighbor_level];
-                const auto iter =
-                  std::find(old_cells.begin(), old_cells.end(), neighbor);
-                old_cells.erase(iter);
-
-                // move neighbor to new level representation
-                const unsigned int new_level = level - max_difference;
-                cells_on_hierarchy_level[new_level].push_back(neighbor);
-
-                // update future FE index
-                neighbor->set_future_fe_index(
-                  fe_index_for_hierarchy_level[new_level]);
-                changed_fe_indices = true;
+                const auto neighbor = cell->neighbor(f);
+
+                const auto neighbor_fe_and_level =
+                  hierarchy_level_for_fe_index.find(
+                    neighbor->future_fe_index());
+
+                // ignore neighbors that are not part of the hierarchy
+                if (neighbor_fe_and_level == hierarchy_level_for_fe_index.end())
+                  continue;
+
+                const unsigned int neighbor_level =
+                  neighbor_fe_and_level->second;
+
+                if ((cell_level - max_difference) > neighbor_level)
+                  {
+                    // update future FE index
+                    const unsigned int new_level = cell_level - max_difference;
+                    neighbor->set_future_fe_index(
+                      fe_index_for_hierarchy_level[new_level]);
+
+                    fe_indices_changed_in_cycle = true;
+                    fe_indices_changed          = true;
+                  }
               }
-          }
+        }
+    }
+  while (fe_indices_changed_in_cycle);
 
-  return changed_fe_indices;
+  return fe_indices_changed;
 }
 
 
-- 
2.39.5