From 9b037871bbf3b6163d90e2172c762b7736027a9f Mon Sep 17 00:00:00 2001
From: David Wells <wellsd2@rpi.edu>
Date: Wed, 2 Mar 2016 20:14:24 -0500
Subject: [PATCH] Redo the logic in make_flux_sparsity_pattern.

As was reported on the mailing list, this function did not work in 1D
because it relied on setting the user flag on a face. This is not
possible because a 'face' in 1D ultimately resolves to the class
TriaAccessor<0, 1, 1>, which does not support user flags.
---
 source/dofs/dof_tools_sparsity.cc      | 25 ++++++++++++++++++-------
 source/dofs/dof_tools_sparsity.inst.in |  4 +---
 2 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc
index a2d927315d..5a98b4e698 100644
--- a/source/dofs/dof_tools_sparsity.cc
+++ b/source/dofs/dof_tools_sparsity.cc
@@ -754,8 +754,6 @@ namespace DoFTools
                 {
                   const typename DoFHandlerType::face_iterator
                   cell_face = cell->face (face_n);
-                  if (cell_face->user_flag_set ())
-                    continue;
 
                   if (cell->at_boundary (face_n))
                     {
@@ -782,12 +780,27 @@ namespace DoFTools
                     {
                       typename DoFHandlerType::level_cell_iterator
                       neighbor = cell->neighbor (face_n);
-                      // Refinement edges are taken care of by coarser cells
-                      if (cell->neighbor_is_coarser (face_n))
+                      // If the cells are on the same level then only add to
+                      // the sparsity pattern if the current cell is 'greater'
+                      // in the total ordering.
+                      if (neighbor->level () == cell->level ()
+                          &&
+                          neighbor->index () > cell->index ())
                         continue;
+                      // If we are more refined then the neighbor, then we
+                      // will automatically find the active neighbor cell when
+                      // we call 'neighbor (face_n)' above. The opposite is
+                      // not true; if the neighbor is more refined then the
+                      // call 'neighbor (face_n)' will *not* return an active
+                      // cell. Hence, only add things to the sparsity pattern
+                      // if the neighbor is coarser than the current cell.
+                      if (neighbor->level () != cell->level ()
+                          &&
+                          !cell->neighbor_is_coarser (face_n))
+                        continue; // (the neighbor is finer)
 
                       const unsigned int
-                      neighbor_face_n = cell->neighbor_of_neighbor (face_n);
+                      neighbor_face_n = cell->neighbor_face_no (face_n);
 
                       if (cell_face->has_children ())
                         {
@@ -864,7 +877,6 @@ namespace DoFTools
                                         }
                                     }
                                 }
-                              sub_neighbor->face(neighbor_face_n)->set_user_flag ();
                             }
                         }
                       else
@@ -933,7 +945,6 @@ namespace DoFTools
                                     }
                                 }
                             }
-                          neighbor->face (neighbor_face_n)->set_user_flag ();
                         }
                     }
                 }
diff --git a/source/dofs/dof_tools_sparsity.inst.in b/source/dofs/dof_tools_sparsity.inst.in
index 1fadf26751..440a1890b2 100644
--- a/source/dofs/dof_tools_sparsity.inst.in
+++ b/source/dofs/dof_tools_sparsity.inst.in
@@ -182,9 +182,6 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS)
      const ConstraintMatrix &constraints,
      const bool, const unsigned int);
 
-
-#if deal_II_dimension > 1
-
     template void
     DoFTools::make_flux_sparsity_pattern<DoFHandler<deal_II_dimension>,SP>
     (const DoFHandler<deal_II_dimension> &dof,
@@ -192,6 +189,7 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS)
      const Table<2,Coupling>&,
      const Table<2,Coupling>&);
 
+#if deal_II_dimension > 1
     template void
     DoFTools::make_flux_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,SP>
     (const hp::DoFHandler<deal_II_dimension> &dof,
-- 
2.39.5