From 4600f351066524ed3f529c788d3aefce24a47cd5 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Fri, 15 Apr 2022 19:41:10 -0400
Subject: [PATCH] Make DoFRenumbering::support_point_wise() about twice as
 fast.

1. There must be an equal number of DoFs per component, so they do not need to
   be calculated.
2. IndexSet::add_indices() works best when indices are sorted, so provide them
   one component at a time.
3. Clean up the output files to match assertions.
---
 source/dofs/dof_renumbering.cc                | 46 +++++++++--
 ..._01.mpirun=2.with_p4est=true.debug.output} |  0
 ...01.mpirun=2.with_p4est=true.release.output | 51 ++++++++++++
 ..._01.mpirun=3.with_p4est=true.debug.output} |  0
 ...01.mpirun=3.with_p4est=true.release.output | 77 +++++++++++++++++++
 ...numbering_01.with_p4est=true.debug.output} |  0
 ...umbering_01.with_p4est=true.release.output | 25 ++++++
 7 files changed, 192 insertions(+), 7 deletions(-)
 rename tests/dofs/{nodal_renumbering_01.mpirun=2.with_p4est=true.output => nodal_renumbering_01.mpirun=2.with_p4est=true.debug.output} (100%)
 create mode 100644 tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.release.output
 rename tests/dofs/{nodal_renumbering_01.mpirun=3.with_p4est=true.output => nodal_renumbering_01.mpirun=3.with_p4est=true.debug.output} (100%)
 create mode 100644 tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.release.output
 rename tests/dofs/{nodal_renumbering_01.with_p4est=true.output => nodal_renumbering_01.with_p4est=true.debug.output} (100%)
 create mode 100644 tests/dofs/nodal_renumbering_01.with_p4est=true.release.output

diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc
index 4567a3139d..a371d9a270 100644
--- a/source/dofs/dof_renumbering.cc
+++ b/source/dofs/dof_renumbering.cc
@@ -2296,8 +2296,21 @@ namespace DoFRenumbering
                                           std::vector<unsigned int>(),
                                           false);
 
-    const std::vector<types::global_dof_index> dofs_per_component =
-      DoFTools::count_dofs_per_fe_component(dof_handler, true);
+#ifdef DEBUG
+    {
+      const std::vector<types::global_dof_index> dofs_per_component =
+        DoFTools::count_dofs_per_fe_component(dof_handler, true);
+      for (const auto &dpc : dofs_per_component)
+        Assert(dofs_per_component[0] == dpc, ExcNotImplemented());
+    }
+#endif
+    const unsigned int n_components =
+      dof_handler.get_fe_collection().n_components();
+    Assert(dof_handler.n_dofs() % n_components == 0, ExcInternalError());
+    const types::global_dof_index dofs_per_component =
+      dof_handler.n_dofs() / n_components;
+    const types::global_dof_index local_dofs_per_component =
+      dof_handler.n_locally_owned_dofs() / n_components;
 
     // At this point we have no more communication to do - simplify things by
     // returning early if possible
@@ -2312,10 +2325,29 @@ namespace DoFRenumbering
     // This index set equals what dof_handler.locally_owned_dofs() would be if
     // we executed the componentwise renumbering.
     IndexSet component_renumbered_dofs(dof_handler.n_dofs());
-    component_renumbered_dofs.add_indices(component_renumbering.begin(),
-                                          component_renumbering.end());
-    for (const auto &dpc : dofs_per_component)
-      AssertThrow(dofs_per_component[0] == dpc, ExcNotImplemented());
+    // DoFs in each component are now consecutive, which IndexSet::add_indices()
+    // can exploit by avoiding calls to sort. Make use of that by adding DoFs
+    // one component at a time:
+    std::vector<types::global_dof_index> component_dofs(
+      local_dofs_per_component);
+    for (unsigned int component = 0; component < n_components; ++component)
+      {
+        for (std::size_t i = 0; i < local_dofs_per_component; ++i)
+          component_dofs[i] =
+            component_renumbering[n_components * i + component];
+        component_renumbered_dofs.add_indices(component_dofs.begin(),
+                                              component_dofs.end());
+      }
+    component_renumbered_dofs.compress();
+#ifdef DEBUG
+    {
+      IndexSet component_renumbered_dofs2(dof_handler.n_dofs());
+      component_renumbered_dofs2.add_indices(component_renumbering.begin(),
+                                             component_renumbering.end());
+      Assert(component_renumbered_dofs2 == component_renumbered_dofs,
+             ExcInternalError());
+    }
+#endif
     for (const FiniteElement<dim, spacedim> &fe :
          dof_handler.get_fe_collection())
       {
@@ -2382,7 +2414,7 @@ namespace DoFRenumbering
                       const auto local_index =
                         component_renumbered_dofs.index_within_set(
                           component_renumbered_cell_dofs[i] +
-                          dofs_per_component[0] * component);
+                          dofs_per_component * component);
 
                       if (component_to_nodal[local_index] ==
                           numbers::invalid_dof_index)
diff --git a/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.output b/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.debug.output
similarity index 100%
rename from tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.output
rename to tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.debug.output
diff --git a/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.release.output b/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.release.output
new file mode 100644
index 0000000000..cb3976230d
--- /dev/null
+++ b/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.release.output
@@ -0,0 +1,51 @@
+
+DEAL:0::new case with locally owned dofs = {[0,7]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::new case with locally owned dofs = {[0,99]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::new case with locally owned dofs = {[0,4033]}
+DEAL:0::
+DEAL:0::difference norm = 1.59444e-15
+DEAL:0::new case with locally owned dofs = {[0,3]}
+DEAL:0::
+DEAL:0::new case with locally owned dofs = {[0,230]}
+DEAL:0::
+DEAL:0::new case with locally owned dofs = {[0,1721]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:0::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:0::The violated condition was:
+DEAL:0::fe.dofs_per_cell == 0 || fe.has_support_points()
+DEAL:0::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:0::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:0::The violated condition was:
+DEAL:0::fe.dofs_per_cell == 0 || fe.has_support_points()
+
+DEAL:1::new case with locally owned dofs = {[8,9]}
+DEAL:1::
+DEAL:1::difference norm = 0.00000
+DEAL:1::new case with locally owned dofs = {[100,249]}
+DEAL:1::
+DEAL:1::difference norm = 0.00000
+DEAL:1::new case with locally owned dofs = {[4034,7965]}
+DEAL:1::
+DEAL:1::difference norm = 1.59444e-15
+DEAL:1::new case with locally owned dofs = {4}
+DEAL:1::
+DEAL:1::new case with locally owned dofs = {[231,440]}
+DEAL:1::
+DEAL:1::new case with locally owned dofs = {[1722,3361]}
+DEAL:1::
+DEAL:1::difference norm = 0.00000
+DEAL:1::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:1::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:1::The violated condition was:
+DEAL:1::fe.dofs_per_cell == 0 || fe.has_support_points()
+DEAL:1::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:1::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:1::The violated condition was:
+DEAL:1::fe.dofs_per_cell == 0 || fe.has_support_points()
+
diff --git a/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.output b/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.debug.output
similarity index 100%
rename from tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.output
rename to tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.debug.output
diff --git a/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.release.output b/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.release.output
new file mode 100644
index 0000000000..6b809e17cb
--- /dev/null
+++ b/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.release.output
@@ -0,0 +1,77 @@
+
+DEAL:0::new case with locally owned dofs = {}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::new case with locally owned dofs = {[0,49]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::new case with locally owned dofs = {[0,2731]}
+DEAL:0::
+DEAL:0::difference norm = 1.67272e-15
+DEAL:0::new case with locally owned dofs = {}
+DEAL:0::
+DEAL:0::new case with locally owned dofs = {[0,120]}
+DEAL:0::
+DEAL:0::new case with locally owned dofs = {[0,1101]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:0::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:0::The violated condition was:
+DEAL:0::fe.dofs_per_cell == 0 || fe.has_support_points()
+DEAL:0::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:0::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:0::The violated condition was:
+DEAL:0::fe.dofs_per_cell == 0 || fe.has_support_points()
+
+DEAL:1::new case with locally owned dofs = {[0,7]}
+DEAL:1::
+DEAL:1::difference norm = 0.00000
+DEAL:1::new case with locally owned dofs = {[50,149]}
+DEAL:1::
+DEAL:1::difference norm = 0.00000
+DEAL:1::new case with locally owned dofs = {[2732,5305]}
+DEAL:1::
+DEAL:1::difference norm = 1.67272e-15
+DEAL:1::new case with locally owned dofs = {[0,3]}
+DEAL:1::
+DEAL:1::new case with locally owned dofs = {[121,340]}
+DEAL:1::
+DEAL:1::new case with locally owned dofs = {[1102,2361]}
+DEAL:1::
+DEAL:1::difference norm = 0.00000
+DEAL:1::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:1::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:1::The violated condition was:
+DEAL:1::fe.dofs_per_cell == 0 || fe.has_support_points()
+DEAL:1::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:1::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:1::The violated condition was:
+DEAL:1::fe.dofs_per_cell == 0 || fe.has_support_points()
+
+
+DEAL:2::new case with locally owned dofs = {[8,9]}
+DEAL:2::
+DEAL:2::difference norm = 0.00000
+DEAL:2::new case with locally owned dofs = {[150,249]}
+DEAL:2::
+DEAL:2::difference norm = 0.00000
+DEAL:2::new case with locally owned dofs = {[5306,7935]}
+DEAL:2::
+DEAL:2::difference norm = 1.67272e-15
+DEAL:2::new case with locally owned dofs = {4}
+DEAL:2::
+DEAL:2::new case with locally owned dofs = {[341,440]}
+DEAL:2::
+DEAL:2::new case with locally owned dofs = {[2362,3361]}
+DEAL:2::
+DEAL:2::difference norm = 0.00000
+DEAL:2::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:2::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:2::The violated condition was:
+DEAL:2::fe.dofs_per_cell == 0 || fe.has_support_points()
+DEAL:2::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:2::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:2::The violated condition was:
+DEAL:2::fe.dofs_per_cell == 0 || fe.has_support_points()
+
diff --git a/tests/dofs/nodal_renumbering_01.with_p4est=true.output b/tests/dofs/nodal_renumbering_01.with_p4est=true.debug.output
similarity index 100%
rename from tests/dofs/nodal_renumbering_01.with_p4est=true.output
rename to tests/dofs/nodal_renumbering_01.with_p4est=true.debug.output
diff --git a/tests/dofs/nodal_renumbering_01.with_p4est=true.release.output b/tests/dofs/nodal_renumbering_01.with_p4est=true.release.output
new file mode 100644
index 0000000000..d3c3ac0e59
--- /dev/null
+++ b/tests/dofs/nodal_renumbering_01.with_p4est=true.release.output
@@ -0,0 +1,25 @@
+
+DEAL:0::new case with locally owned dofs = {[0,9]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::new case with locally owned dofs = {[0,249]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::new case with locally owned dofs = {[0,7941]}
+DEAL:0::
+DEAL:0::difference norm = 1.62791e-15
+DEAL:0::new case with locally owned dofs = {[0,4]}
+DEAL:0::
+DEAL:0::new case with locally owned dofs = {[0,440]}
+DEAL:0::
+DEAL:0::new case with locally owned dofs = {[0,3361]}
+DEAL:0::
+DEAL:0::difference norm = 0.00000
+DEAL:0::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:0::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:0::The violated condition was:
+DEAL:0::fe.dofs_per_cell == 0 || fe.has_support_points()
+DEAL:0::FE_NedelecSZ nodal renumbering failed successfully with message:
+DEAL:0::void dealii::DoFRenumbering::compute_support_point_wise(std::vector<unsigned int>&, const dealii::DoFHandler<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
+DEAL:0::The violated condition was:
+DEAL:0::fe.dofs_per_cell == 0 || fe.has_support_points()
-- 
2.39.5