Use faster data structure to compute matrix-free renumbering

diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc
index 438c9beac50fefcd4b0fdbad0db1e21054d7d0da..a924e551d9a102ae6f6dd736be82952819eef0e3 100644 (file)
--- a/source/dofs/dof_renumbering.cc
+++ b/source/dofs/dof_renumbering.cc
@@ -2679,10 +2679,10 @@ namespace DoFRenumbering
                                     const IndexSet &              owned_dofs,
                                     std::vector<unsigned int> &   result,
                                     unsigned int &counter_dof_numbers) {
-        if (owned_dofs.is_element(dof_index))
+        const types::global_dof_index local_dof_index =
+          owned_dofs.index_within_set(dof_index);
+        if (local_dof_index != numbers::invalid_dof_index)
           {
-            const unsigned int local_dof_index =
-              owned_dofs.index_within_set(dof_index);
             AssertIndexRange(local_dof_index, result.size());
             if (result[local_dof_index] == numbers::invalid_unsigned_int)
               result[local_dof_index] = counter_dof_numbers++;
@@ -2690,7 +2690,7 @@ namespace DoFRenumbering
       };
 
       unsigned int                         counter_dof_numbers = 0;
-      std::vector<unsigned int>            local_dofs, dofs_extracted;
+      std::vector<unsigned int>            dofs_extracted;
       std::vector<types::global_dof_index> dof_indices(
         matrix_free.get_dof_handler(component).get_fe().dofs_per_cell);
 
@@ -2705,19 +2705,21 @@ namespace DoFRenumbering
       // that we pick the unconstrained indices in the hierarchical order for
       // part (a) as this makes it easy to identify the DoFs on the individual
       // entities, whereas we pick the numbers with constraints eliminated for
-      // part (b). For the latter, we finally need to sort the indices and
-      // remove duplicates to really only count the DoFs in each range of cell
-      // batches once.
+      // part (b). For the latter, we keep track of each DoF's interaction
+      // with different ranges of cell batches, i.e., call-backs into the
+      // operation_on_cell_range() function.
+      const unsigned int        n_owned_dofs = owned_dofs.n_elements();
       std::vector<unsigned int> dof_numbers_mf_order(
-        owned_dofs.n_elements(), dealii::numbers::invalid_unsigned_int);
-      std::vector<unsigned char> touch_count(owned_dofs.n_elements());
+        n_owned_dofs, dealii::numbers::invalid_unsigned_int);
+      std::vector<unsigned int> last_touch_by_cell_batch_range(
+        n_owned_dofs, dealii::numbers::invalid_unsigned_int);
+      std::vector<unsigned char> touch_count(n_owned_dofs);
 
       const auto operation_on_cell_range =
         [&](const MatrixFree<dim, Number, VectorizedArrayType> &data,
             unsigned int &,
             const unsigned int &,
             const std::pair<unsigned int, unsigned int> &cell_range) {
-          local_dofs.clear();
           for (unsigned int cell = cell_range.first; cell < cell_range.second;
                ++cell)
             {
@@ -2759,22 +2761,22 @@ namespace DoFRenumbering
                     }
                 }
 
-              // part (b): compute the touch count in the current cell batch
-              // range
+              // part (b): increment the touch count of a dof appearing in the
+              // current cell batch if it was last touched by another than the
+              // present cell batch range (we track them via storing the last
+              // cell batch range that touched a particular dof)
               data.get_dof_info(component).get_dof_indices_on_cell_batch(
                 dofs_extracted, cell);
-              local_dofs.insert(local_dofs.end(),
-                                dofs_extracted.begin(),
-                                dofs_extracted.end());
+              for (unsigned int dof_index : dofs_extracted)
+                if (dof_index < n_owned_dofs &&
+                    last_touch_by_cell_batch_range[dof_index] !=
+                      cell_range.first)
+                  {
+                    ++touch_count[dof_index];
+                    last_touch_by_cell_batch_range[dof_index] =
+                      cell_range.first;
+                  }
             }
-
-          std::sort(local_dofs.begin(), local_dofs.end());
-          local_dofs.erase(std::unique(local_dofs.begin(), local_dofs.end()),
-                           local_dofs.end());
-
-          for (unsigned int i : local_dofs)
-            if (i < touch_count.size())
-              touch_count[i]++;
         };
 
       // Finally run the matrix-free loop.
@@ -2786,7 +2788,7 @@ namespace DoFRenumbering
       matrix_free.template cell_loop<unsigned int, unsigned int>(
         operation_on_cell_range, counter_dof_numbers, counter_dof_numbers);
 
-      AssertDimension(counter_dof_numbers, owned_dofs.n_elements());
+      AssertDimension(counter_dof_numbers, n_owned_dofs);
 
       return std::make_pair(dof_numbers_mf_order, touch_count);
     }