Merge pull request #14868 from vyushut/flux_sparsity_tests

make_flux_sparsity_pattern() revision

source/dofs/dof_tools_sparsity.cc

@@ -852,9 +852,6 @@ namespace DoFTools
                         }
                       else
                         {
-                          if (!face_has_flux_coupling(cell, face_n))
-                            continue;
-
                           typename DoFHandler<dim,
                                               spacedim>::level_cell_iterator
                             neighbor =
@@ -888,6 +885,10 @@ namespace DoFTools
                               neighbor->is_locally_owned())
                             continue; // (the neighbor is finer)
 
+                          if (!face_has_flux_coupling(cell, face_n))
+                            continue;
+
+
                           const unsigned int neighbor_face_n =
                             periodic_neighbor ?
                               cell->periodic_neighbor_face_no(face_n) :
@@ -1197,9 +1198,6 @@ namespace DoFTools
                             level_cell_iterator neighbor =
                               cell->neighbor_or_periodic_neighbor(face);
 
-                          if (!face_has_flux_coupling(cell, face))
-                            continue;
-
                           // Like the non-hp-case: If the cells are on the same
                           // level (and both are active, locally-owned cells)
                           // then only add to the sparsity pattern if the
@@ -1229,6 +1227,10 @@ namespace DoFTools
                               neighbor->is_locally_owned())
                             continue; // (the neighbor is finer)
 
+
+                          if (!face_has_flux_coupling(cell, face))
+                            continue;
+
                           // In 1d, go straight to the cell behind this
                           // particular cell's most terminal cell. This makes us
                           // skip the if (neighbor->has_children()) section

tests/sparsity/flux_sparsity_pattern_visiting_once.cc

@@ -0,0 +1,105 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2021 - 2022 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// @sect3{This test considers two elements which share a common, regular face.
+// The logic of make_flux_sparsity_pattern() loops over all faces and filters
+// some of them according to a predicate face_has_flux_coupling(). We check
+// that the only internal face in the current setup is visited exactly once and
+// that the predicate face_has_flux_coupling() is evaluated once for this face.}
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+Triangulation<dim>
+make_two_elements()
+{
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, -2.0, 2.0);
+  triangulation.begin_active()->set_refine_flag(
+    RefinementCase<dim>::cut_axis(0));
+  triangulation.execute_coarsening_and_refinement();
+  return triangulation;
+}
+
+
+template <int dim>
+bool
+is_face_on_OY(const typename DoFHandler<dim>::active_cell_iterator &cell,
+              const unsigned int                                    face_index)
+{
+  deallog
+    << "This sentence should appear once when the corresponding face is visited only once on cell "
+    << cell->index() << std::endl;
+  return (std::abs(cell->face(face_index)->center()(0)) < 0.01);
+}
+
+template <int dim>
+void
+create_and_output_flux_sparsity_with_filter(
+  DoFHandler<dim> &                       dof_handler,
+  std::function<bool(const typename DoFHandler<dim>::active_cell_iterator,
+                     const unsigned int)> filter)
+{
+  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+  DoFTools::make_flux_sparsity_pattern(
+    dof_handler,
+    dsp,
+    AffineConstraints<double>(),
+    true,
+    Table<2, DoFTools::Coupling>(1, 1, new auto(DoFTools::always)),
+    Table<2, DoFTools::Coupling>(1, 1, new auto(DoFTools::always)),
+    numbers::invalid_subdomain_id,
+    [&](const auto &cell, const unsigned int face_index) {
+      return filter(cell, face_index);
+    });
+  dsp.print(deallog.get_file_stream());
+}
+
+template <int dim>
+void
+check()
+{
+  // create a square with two elements
+  const Triangulation<dim> tria = make_two_elements<dim>();
+  // Generate Q1 dofs for the mesh grid
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(FE_Q<dim>(1));
+  // create and output the sparsity pattern by using a filter
+  create_and_output_flux_sparsity_with_filter<dim>(dof_handler,
+                                                   is_face_on_OY<dim>);
+}
+
+int
+main()
+{
+  initlog();
+  deallog.push("2d");
+  check<2>();
+  deallog.pop();
+  deallog.push("3d");
+  check<3>();
+  deallog.pop();
+}

tests/sparsity/flux_sparsity_pattern_visiting_once.output

@@ -0,0 +1,21 @@
+
+DEAL:2d::This sentence should appear once when the corresponding face is visited only once on cell 1
+[0,0,1,2,3,4,5]
+[1,0,1,2,3,4,5]
+[2,0,1,2,3,4,5]
+[3,0,1,2,3,4,5]
+[4,0,1,2,3,4,5]
+[5,0,1,2,3,4,5]
+DEAL:3d::This sentence should appear once when the corresponding face is visited only once on cell 1
+[0,0,1,2,3,4,5,6,7,8,9,10,11]
+[1,0,1,2,3,4,5,6,7,8,9,10,11]
+[2,0,1,2,3,4,5,6,7,8,9,10,11]
+[3,0,1,2,3,4,5,6,7,8,9,10,11]
+[4,0,1,2,3,4,5,6,7,8,9,10,11]
+[5,0,1,2,3,4,5,6,7,8,9,10,11]
+[6,0,1,2,3,4,5,6,7,8,9,10,11]
+[7,0,1,2,3,4,5,6,7,8,9,10,11]
+[8,0,1,2,3,4,5,6,7,8,9,10,11]
+[9,0,1,2,3,4,5,6,7,8,9,10,11]
+[10,0,1,2,3,4,5,6,7,8,9,10,11]
+[11,0,1,2,3,4,5,6,7,8,9,10,11]