fix make_flux_sparsity_pattern for FENothing

diff --git a/doc/news/changes/minor/20240219SchreterRitthaler b/doc/news/changes/minor/20240219SchreterRitthaler
new file mode 100644 (file)
index 0000000..bf1f898
--- /dev/null
+++ b/doc/news/changes/minor/20240219SchreterRitthaler
@@ -0,0 +1,3 @@
+Fixed: Fix make_flux_sparsity_pattern for a FECollection containing FENothing elements. 
+<br>
+(Magdalena Schreter-Fleischhacker, Andreas Ritthaler, 2024/02/19)

diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h
index c6d24ec4969c2f0a5fb5520d32e868ee81a5cf32..359196623fd792a952a1f721f1b05a9e11ab4095 100644 (file)
--- a/include/deal.II/lac/affine_constraints.templates.h
+++ b/include/deal.II/lac/affine_constraints.templates.h
@@ -4714,12 +4714,17 @@ AffineConstraints<number>::add_entries_local_to_global(
   const size_type n_local_rows = row_indices.size();
   const size_type n_local_cols = col_indices.size();
 
+  // Early return if the length of row and column indices is zero, relevant
+  // for the usage with FENothing.
+  if (n_local_cols == 0 && n_local_rows == 0)
+    return;
+
   typename internal::AffineConstraints::ScratchDataAccessor<number>
     scratch_data(this->scratch_data);
   std::vector<std::pair<size_type, size_type>> &cell_entries =
     scratch_data->new_entries;
   cell_entries.resize(0);
-  cell_entries.reserve(row_indices.size() * col_indices.size());
+  cell_entries.reserve(n_local_rows * n_local_cols);
 
   // if constrained entries should be kept, need to add rows and columns of
   // those to the sparsity pattern

diff --git a/tests/hp/flux_sparsity_03.cc b/tests/hp/flux_sparsity_03.cc
new file mode 100644 (file)
index 0000000..73e3028
--- /dev/null
+++ b/tests/hp/flux_sparsity_03.cc
@@ -0,0 +1,144 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2020 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.
+//
+// ---------------------------------------------------------------------
+
+// A test that checks make_flux_sparsity_pattern
+// with a FECollection containing FE_Q and FE_Nothing elements.
+
+//   x__________x__________ _________
+//   |          |          |         |
+//   |    FE_Q  |   FE_N   |   FE_N  |     x: constrained DoFs
+//   |          |          |         |
+//   |__________|__________|_________|
+//   x          x
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+template <int dim>
+void
+check()
+{
+  // create triangulation
+  Triangulation<dim>        triangulation;
+  std::vector<unsigned int> subdivisions(dim, 1U);
+  subdivisions[0] = 3;
+  subdivisions[1] = 1;
+  if (dim == 3)
+    subdivisions[2] = 1;
+  Point<dim> p1, p2;
+  switch (dim)
+    {
+      case 2:
+        p1[0] = p1[1] = 0.0;
+        p2[0]         = 3.0;
+        p2[1]         = 1.0;
+        break;
+      case 3:
+        p1[0] = p1[1] = p1[2] = 0.0;
+        p2[0]                 = 3.0;
+        p2[1] = p2[2] = 1.0;
+        break;
+      default:
+        DEAL_II_NOT_IMPLEMENTED();
+    }
+
+  GridGenerator::subdivided_hyper_rectangle(triangulation,
+                                            subdivisions,
+                                            p1,
+                                            p2);
+
+  // create FE Collection and insert two FE objects
+  // Q1 and one FE object Nothing
+  hp::FECollection<dim> fe_collection;
+
+  fe_collection.push_back(FE_Q<dim>(1));
+  fe_collection.push_back(FE_Nothing<dim>());
+
+  // set-up DoFHandler
+  DoFHandler<dim> dof_handler(triangulation);
+
+  // set active fe index
+  unsigned int i = 0;
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      if (i == 2)
+        {
+          cell->set_active_fe_index(0);
+          ++i;
+        }
+      else
+        {
+          cell->set_active_fe_index(1);
+          ++i;
+        }
+    }
+
+  dof_handler.distribute_dofs(fe_collection);
+
+  {
+    // set constraints
+    AffineConstraints<double> constraints;
+
+    for (unsigned i = 0; i < Utilities::pow<unsigned int>(2, dim); ++i)
+      constraints.add_constraint(i, {}, i);
+
+    constraints.close();
+
+    // setup sparsity pattern
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+    DoFTools::make_flux_sparsity_pattern(dof_handler, dsp, constraints);
+    dsp.compress();
+
+    // print sparsity pattern
+    deallog << "nonzero matrix elements: " << dsp.n_nonzero_elements()
+            << std::endl;
+    dsp.print(deallog.get_file_stream());
+  }
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("2d");
+  check<2>();
+  deallog.pop();
+  deallog.push("3d");
+  check<3>();
+  deallog.pop();
+}

diff --git a/tests/hp/flux_sparsity_03.output b/tests/hp/flux_sparsity_03.output
new file mode 100644 (file)
index 0000000..adbee02
--- /dev/null
+++ b/tests/hp/flux_sparsity_03.output
@@ -0,0 +1,15 @@
+
+DEAL:2d::nonzero matrix elements: 16
+[0,0,1,2,3]
+[1,0,1,2,3]
+[2,0,1,2,3]
+[3,0,1,2,3]
+DEAL:3d::nonzero matrix elements: 64
+[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]