Simplify DoFTools::make_flux_sparsity_pattern and test more cases

diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc
index d4f6f9babc21eeb88a210eab2bf7018e5f72bede..c7e392b1eb58a86ef0ec15657710e8dc14d782bb 100644 (file)
--- a/source/dofs/dof_tools_sparsity.cc
+++ b/source/dofs/dof_tools_sparsity.cc
@@ -983,21 +983,36 @@ namespace DoFTools
         std::vector<types::global_dof_index> dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof));
         std::vector<types::global_dof_index> dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof));
 
-        std::vector<Table<2,Coupling> > int_dof_mask (fe.size());
-
-        int_dof_mask = dof_couplings_from_component_couplings (fe, int_mask);
-
-        // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
-        // to constraints.add_entries_local_to_global()
-        std::vector<Table<2,bool> > bool_int_dof_mask (fe.size());
+        const unsigned int n_components = fe.n_components();
+        AssertDimension (int_mask.size(0), n_components);
+        AssertDimension (int_mask.size(1), n_components);
+        AssertDimension (flux_mask.size(0), n_components);
+        AssertDimension (flux_mask.size(1), n_components);
+
+        // note that we also need to set the respective entries if flux_mask
+        // says so. this is necessary since we need to consider all degrees of
+        // freedom on a cell for interior faces.
+        Table<2,Coupling> int_and_flux_mask (n_components, n_components);
+        for (unsigned int c1=0; c1<n_components; ++c1)
+          for (unsigned int c2=0; c2<n_components; ++c2)
+            if (int_mask(c1,c2) != none || flux_mask(c1,c2) != none)
+              int_and_flux_mask(c1,c2) = always;
+
+        std::vector<Table<2,Coupling> > int_and_flux_dof_mask
+          = dof_couplings_from_component_couplings (fe, int_and_flux_mask);
+
+        // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we
+        // can pass it to constraints.add_entries_local_to_global()
+        std::vector<Table<2,bool> > bool_int_and_flux_dof_mask (fe.size());
         for (unsigned int f=0; f<fe.size(); ++f)
           {
-            bool_int_dof_mask[f].reinit (TableIndices<2>(fe[f].dofs_per_cell,fe[f].dofs_per_cell));
-            bool_int_dof_mask[f].fill (false);
+            bool_int_and_flux_dof_mask[f].reinit (TableIndices<2>(fe[f].dofs_per_cell,
+                                                                  fe[f].dofs_per_cell));
+            bool_int_and_flux_dof_mask[f].fill (false);
             for (unsigned int i=0; i<fe[f].dofs_per_cell; ++i)
               for (unsigned int j=0; j<fe[f].dofs_per_cell; ++j)
-                if (int_dof_mask[f](i,j) != none)
-                  bool_int_dof_mask[f](i,j) = true;
+                if (int_and_flux_dof_mask[f](i,j) != none)
+                  bool_int_and_flux_dof_mask[f](i,j) = true;
           }
 
 
@@ -1014,12 +1029,14 @@ namespace DoFTools
               dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell);
               cell->get_dof_indices (dofs_on_this_cell);
 
-              // make sparsity pattern for this cell
+              // make sparsity pattern for this cell also taking into account
+              // the couplings due to face contributions on the same cell
               constraints.add_entries_local_to_global (dofs_on_this_cell,
                                                        sparsity,
                                                        keep_constrained_dofs,
-                                                       bool_int_dof_mask[cell->active_fe_index()]);
-              // Loop over all interior neighbors
+                                                       bool_int_and_flux_dof_mask[cell->active_fe_index()]);
+
+              // Loop over interior faces
               for (unsigned int face = 0;
                    face < GeometryInfo<dim>::faces_per_cell;
                    ++face)
@@ -1029,35 +1046,6 @@ namespace DoFTools
 
                   const bool periodic_neighbor = cell->has_periodic_neighbor (face);
 
-                  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
-                    {
-                      const unsigned int ii
-                        = (cell->get_fe().is_primitive(i) ?
-                           cell->get_fe().system_to_component_index(i).first
-                           :
-                           cell->get_fe().get_nonzero_components(i).first_selected_component()
-                          );
-
-                      Assert (ii < cell->get_fe().n_components(), ExcInternalError());
-
-                      for (unsigned int j=0; j<cell->get_fe().dofs_per_cell; ++j)
-                        {
-                          const unsigned int jj
-                            = (cell->get_fe().is_primitive(j) ?
-                               cell->get_fe().system_to_component_index(j).first
-                               :
-                               cell->get_fe().get_nonzero_components(j).first_selected_component()
-                              );
-
-                          Assert (jj < cell->get_fe().n_components(), ExcInternalError());
-
-                          if ((flux_mask(ii,jj) == always)
-                              ||
-                              (flux_mask(ii,jj) == nonzero))
-                            sparsity.add (dofs_on_this_cell[i],
-                                          dofs_on_this_cell[j]);
-                        }
-                    }
                   if ((!cell->at_boundary(face)) || periodic_neighbor)
                     {
                       typename dealii::hp::DoFHandler<dim,spacedim>::level_cell_iterator

diff --git a/tests/hp/flux_sparsity_02.cc b/tests/hp/flux_sparsity_02.cc
index f9757fac9ffc5c5ba6d65122a70adb3fa476a71d..aab2d77f908aba83df441012fd55ac67f50fed7a 100644 (file)
--- a/tests/hp/flux_sparsity_02.cc
+++ b/tests/hp/flux_sparsity_02.cc
@@ -73,14 +73,75 @@ check ()
   dof_handler.begin_active()->set_active_fe_index(1);
   dof_handler.distribute_dofs(fe_collection);
 
-  DynamicSparsityPattern dsp(dof_handler.n_dofs(),
-                             dof_handler.n_dofs());
-
-  DoFTools::make_flux_sparsity_pattern (dof_handler, dsp);
-  dsp.compress();
-
-  // Print sparsity pattern, we expect that all dofs couple with each other.
-  dsp.print(deallog.get_file_stream());
+  deallog << "Dimension " << dim << std::endl;
+
+  {
+    deallog << "cell and face coupling" << std::endl;
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(),
+                               dof_handler.n_dofs());
+
+    Table<2, DoFTools::Coupling> cell_coupling(1,1);
+    Table<2, DoFTools::Coupling> face_coupling(1,1);
+    cell_coupling(0,0) = DoFTools::always;
+    face_coupling(0,0) = DoFTools::always;
+
+    DoFTools::make_flux_sparsity_pattern (dof_handler, dsp, cell_coupling, face_coupling);
+    dsp.compress();
+
+    // Print sparsity pattern, we expect that all dofs couple with each other.
+    dsp.print(deallog.get_file_stream());
+  }
+
+  {
+    deallog << "cell coupling" << std::endl;
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(),
+                               dof_handler.n_dofs());
+
+    Table<2, DoFTools::Coupling> cell_coupling(1,1);
+    Table<2, DoFTools::Coupling> face_coupling(1,1);
+    cell_coupling(0,0) = DoFTools::always;
+    face_coupling(0,0) = DoFTools::none;
+
+    DoFTools::make_flux_sparsity_pattern (dof_handler, dsp, cell_coupling, face_coupling);
+    dsp.compress();
+
+    // Print sparsity pattern, we expect no couplings across the face.
+    dsp.print(deallog.get_file_stream());
+  }
+
+  {
+    deallog << "face coupling" << std::endl;
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(),
+                               dof_handler.n_dofs());
+
+    Table<2, DoFTools::Coupling> cell_coupling(1,1);
+    Table<2, DoFTools::Coupling> face_coupling(1,1);
+    cell_coupling(0,0) = DoFTools::none;
+    face_coupling(0,0) = DoFTools::always;
+
+    DoFTools::make_flux_sparsity_pattern (dof_handler, dsp, cell_coupling, face_coupling);
+    dsp.compress();
+
+    // Print sparsity pattern, we expect that all dofs couple with each other.
+    dsp.print(deallog.get_file_stream());
+  }
+
+  {
+    deallog << "no coupling" << std::endl;
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(),
+                               dof_handler.n_dofs());
+
+    Table<2, DoFTools::Coupling> cell_coupling(1,1);
+    Table<2, DoFTools::Coupling> face_coupling(1,1);
+    cell_coupling(0,0) = DoFTools::none;
+    face_coupling(0,0) = DoFTools::none;
+
+    DoFTools::make_flux_sparsity_pattern (dof_handler, dsp, cell_coupling, face_coupling);
+    dsp.compress();
+
+    // Print sparsity pattern, we expect no couplings.
+    dsp.print(deallog.get_file_stream());
+  }
 }
 
 

diff --git a/tests/hp/flux_sparsity_02.output b/tests/hp/flux_sparsity_02.output
index 26fea0c98501ce8101b83b2f72dbefdd25de0058..98b1f3efa4fa47d890ee46d367fd0edad24500ec 100644 (file)
--- a/tests/hp/flux_sparsity_02.output
+++ b/tests/hp/flux_sparsity_02.output
@@ -1,9 +1,51 @@
 
+DEAL:2d::Dimension 2
+DEAL:2d::cell and face coupling
 [0,0,1,2,3,4]
 [1,0,1,2,3,4]
 [2,0,1,2,3,4]
 [3,0,1,2,3,4]
 [4,0,1,2,3,4]
+DEAL:2d::cell coupling
+[0,0,1,2,3]
+[1,0,1,2,3]
+[2,0,1,2,3]
+[3,0,1,2,3]
+[4,4]
+DEAL:2d::face coupling
+[0,0,1,2,3,4]
+[1,0,1,2,3,4]
+[2,0,1,2,3,4]
+[3,0,1,2,3,4]
+[4,0,1,2,3,4]
+DEAL:2d::no coupling
+[0]
+[1]
+[2]
+[3]
+[4]
+DEAL:3d::Dimension 3
+DEAL:3d::cell and face coupling
+[0,0,1,2,3,4,5,6,7,8]
+[1,0,1,2,3,4,5,6,7,8]
+[2,0,1,2,3,4,5,6,7,8]
+[3,0,1,2,3,4,5,6,7,8]
+[4,0,1,2,3,4,5,6,7,8]
+[5,0,1,2,3,4,5,6,7,8]
+[6,0,1,2,3,4,5,6,7,8]
+[7,0,1,2,3,4,5,6,7,8]
+[8,0,1,2,3,4,5,6,7,8]
+DEAL:3d::cell coupling
+[0,0,1,2,3,4,5,6,7]
+[1,0,1,2,3,4,5,6,7]
+[2,0,1,2,3,4,5,6,7]
+[3,0,1,2,3,4,5,6,7]
+[4,0,1,2,3,4,5,6,7]
+[5,0,1,2,3,4,5,6,7]
+[6,0,1,2,3,4,5,6,7]
+[7,0,1,2,3,4,5,6,7]
+[8,8]
+DEAL:3d::face coupling
 [0,0,1,2,3,4,5,6,7,8]
 [1,0,1,2,3,4,5,6,7,8]
 [2,0,1,2,3,4,5,6,7,8]
@@ -13,3 +55,13 @@
 [6,0,1,2,3,4,5,6,7,8]
 [7,0,1,2,3,4,5,6,7,8]
 [8,0,1,2,3,4,5,6,7,8]
+DEAL:3d::no coupling
+[0]
+[1]
+[2]
+[3]
+[4]
+[5]
+[6]
+[7]
+[8]