From b5ab04f912d6259eeb33a815221ef70bb116480f Mon Sep 17 00:00:00 2001
From: Mathias Anselmann <mathias.anselmann@posteo.de>
Date: Mon, 7 Jun 2021 22:48:34 +0200
Subject: [PATCH] changed make_flux_sparsity_pattern to have the ability to add
 constraints

---
 include/deal.II/multigrid/mg_tools.h | 15 +++++++---
 source/multigrid/mg_tools.cc         | 44 +++++++++++++++-------------
 source/multigrid/mg_tools.inst.in    | 23 ++++++++++++---
 3 files changed, 54 insertions(+), 28 deletions(-)

diff --git a/include/deal.II/multigrid/mg_tools.h b/include/deal.II/multigrid/mg_tools.h
index c5cf619e3e..1b0b2d7f05 100644
--- a/include/deal.II/multigrid/mg_tools.h
+++ b/include/deal.II/multigrid/mg_tools.h
@@ -103,11 +103,18 @@ namespace MGTools
    * and
    * @ref DoFTools
    */
-  template <int dim, typename SparsityPatternType, int spacedim>
+  template <int dim,
+            int spacedim,
+            typename SparsityPatternType,
+            typename number = double>
   void
-  make_flux_sparsity_pattern(const DoFHandler<dim, spacedim> &dof_handler,
-                             SparsityPatternType &            sparsity,
-                             const unsigned int               level);
+  make_flux_sparsity_pattern(
+    const DoFHandler<dim, spacedim> &dof_handler,
+    SparsityPatternType &            sparsity,
+    const unsigned int               level,
+    const AffineConstraints<number> &constraints = AffineConstraints<number>(),
+    const bool                       keep_constrained_dofs = true);
+
 
   /**
    * Create sparsity pattern for the fluxes at refinement edges. The matrix
diff --git a/source/multigrid/mg_tools.cc b/source/multigrid/mg_tools.cc
index 2393a661a7..02f4fb12b1 100644
--- a/source/multigrid/mg_tools.cc
+++ b/source/multigrid/mg_tools.cc
@@ -623,11 +623,16 @@ namespace MGTools
 
 
 
-  template <int dim, typename SparsityPatternType, int spacedim>
+  template <int dim,
+            int spacedim,
+            typename SparsityPatternType,
+            typename number>
   void
   make_flux_sparsity_pattern(const DoFHandler<dim, spacedim> &dof,
                              SparsityPatternType &            sparsity,
-                             const unsigned int               level)
+                             const unsigned int               level,
+                             const AffineConstraints<number> &constraints,
+                             const bool keep_constrained_dofs)
   {
     const types::global_dof_index n_dofs = dof.n_dofs(level);
     (void)n_dofs;
@@ -649,9 +654,9 @@ namespace MGTools
 
         cell->get_mg_dof_indices(dofs_on_this_cell);
         // make sparsity pattern for this cell
-        for (unsigned int i = 0; i < dofs_per_cell; ++i)
-          for (unsigned int j = 0; j < dofs_per_cell; ++j)
-            sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
+        constraints.add_entries_local_to_global(dofs_on_this_cell,
+                                                sparsity,
+                                                keep_constrained_dofs);
 
         // Loop over all interior neighbors
         for (const unsigned int face : GeometryInfo<dim>::face_indices())
@@ -674,23 +679,22 @@ namespace MGTools
                 // only add one direction The other is taken care of by
                 // neighbor (except when the neighbor is not owned by the same
                 // processor)
-                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                constraints.add_entries_local_to_global(dofs_on_this_cell,
+                                                        dofs_on_other_cell,
+                                                        sparsity,
+                                                        keep_constrained_dofs);
+
+                if (neighbor->is_locally_owned_on_level() == false)
                   {
-                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
-                      {
-                        sparsity.add(dofs_on_this_cell[i],
-                                     dofs_on_other_cell[j]);
-                      }
+                    constraints.add_entries_local_to_global(
+                      dofs_on_other_cell, sparsity, keep_constrained_dofs);
+
+                    constraints.add_entries_local_to_global(
+                      dofs_on_other_cell,
+                      dofs_on_this_cell,
+                      sparsity,
+                      keep_constrained_dofs);
                   }
-                if (neighbor->is_locally_owned_on_level() == false)
-                  for (unsigned int i = 0; i < dofs_per_cell; ++i)
-                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
-                      {
-                        sparsity.add(dofs_on_other_cell[i],
-                                     dofs_on_other_cell[j]);
-                        sparsity.add(dofs_on_other_cell[i],
-                                     dofs_on_this_cell[j]);
-                      }
               }
           }
       }
diff --git a/source/multigrid/mg_tools.inst.in b/source/multigrid/mg_tools.inst.in
index 3f68d8ca78..f0513a8f5b 100644
--- a/source/multigrid/mg_tools.inst.in
+++ b/source/multigrid/mg_tools.inst.in
@@ -43,6 +43,7 @@ for (PATTERN : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS;
         const unsigned int,
         const AffineConstraints<float> &,
         const bool);
+
 #endif
     \}
   }
@@ -67,10 +68,24 @@ for (PATTERN : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS;
 
 #if deal_II_dimension == deal_II_space_dimension
       template void
-      make_flux_sparsity_pattern<deal_II_dimension>(
-        const DoFHandler<deal_II_dimension> &,
-        PATTERN &,
-        const unsigned int);
+      make_flux_sparsity_pattern<deal_II_dimension,
+                                 deal_II_space_dimension,
+                                 PATTERN,
+                                 double>(const DoFHandler<deal_II_dimension> &,
+                                         PATTERN &,
+                                         const unsigned int,
+                                         const AffineConstraints<double> &,
+                                         const bool);
+
+      template void
+      make_flux_sparsity_pattern<deal_II_dimension,
+                                 deal_II_space_dimension,
+                                 PATTERN,
+                                 float>(const DoFHandler<deal_II_dimension> &,
+                                        PATTERN &,
+                                        const unsigned int,
+                                        const AffineConstraints<float> &,
+                                        const bool);
 
       template void
       make_flux_sparsity_pattern_edge<deal_II_dimension>(
-- 
2.39.5