Check during MF::reinit() if AffineConstraints are closed

diff --git a/doc/news/changes/incompatibilities/20220308Munch b/doc/news/changes/incompatibilities/20220308Munch
new file mode 100644 (file)
index 0000000..b933777
--- /dev/null
+++ b/doc/news/changes/incompatibilities/20220308Munch
@@ -0,0 +1,9 @@
+Changed: To be able to guarantee correctness of
+the application of constraints in MatrixFree, we
+require that users call AffineConstraints::close()
+before calling MatrixFree::reinit() if AffineConstraints
+objects are passed that are not empty. One can check
+if an AffineConstraints object is closed with the
+new function AffineConstraints::is_closed(). 
+<br>
+(Peter Munch, 2022/03/08)

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 51a6b0a12057c0ab4ef41340be29a11d488439bb..15da51cd2e3967eff8fcff53db8e27d69702bf7a 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -372,7 +372,7 @@ void
 MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
   const std::shared_ptr<hp::MappingCollection<dim>> &    mapping,
   const std::vector<const DoFHandler<dim, dim> *> &      dof_handler,
-  const std::vector<const AffineConstraints<number2> *> &constraint,
+  const std::vector<const AffineConstraints<number2> *> &constraints,
   const std::vector<IndexSet> &                          locally_owned_dofs,
   const std::vector<hp::QCollection<q_dim>> &            quad,
   const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
@@ -414,7 +414,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
     {
       clear();
       Assert(dof_handler.size() > 0, ExcMessage("No DoFHandler is given."));
-      AssertDimension(dof_handler.size(), constraint.size());
+      AssertDimension(dof_handler.size(), constraints.size());
       AssertDimension(dof_handler.size(), locally_owned_dofs.size());
 
       task_info.allow_ghosted_vectors_in_loops =
@@ -439,6 +439,16 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
           task_info.n_procs         = 1;
         }
 
+#ifdef DEBUG
+      for (const auto &constraint : constraints)
+        Assert(
+          constraint->is_closed(task_info.communicator),
+          ExcMessage(
+            "You have proveded a non-empty AffineConstraint object that has not "
+            "been closed. Please call AffineConstraints::close() before "
+            "calling MatrixFree::reinit()!"));
+#endif
+
       initialize_dof_handlers(dof_handler, additional_data);
       for (unsigned int no = 0; no < dof_handler.size(); ++no)
         {
@@ -469,7 +479,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
       // constraint_pool_data. It also reorders the way cells are gone through
       // (to separate cells with overlap to other processors from others
       // without).
-      initialize_indices(constraint, locally_owned_dofs, additional_data);
+      initialize_indices(constraints, locally_owned_dofs, additional_data);
     }
 
   // initialize bare structures

diff --git a/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc b/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc
index 6a73585ea9648d6d10372666501aaa8df18fee16..a9f6ae2ae1cdb50770801b52bff5c9d2b61b588f 100644 (file)
--- a/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc
+++ b/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc
@@ -106,6 +106,7 @@ main()
 
   AffineConstraints<double> constraints;
   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+  constraints.close();
 
   const auto print = [](const auto &label, const auto &matrix_free) {
     deallog << label << std::endl;

diff --git a/tests/matrix_free/mixed_mesh_hn_algorithm_01.output b/tests/matrix_free/mixed_mesh_hn_algorithm_01.output
index 384befad9035801ae170d755d9c91aded384b136..540231691334fa1fc1c90d6cab72216aaf4d55e4 100644 (file)
--- a/tests/matrix_free/mixed_mesh_hn_algorithm_01.output
+++ b/tests/matrix_free/mixed_mesh_hn_algorithm_01.output
@@ -1,9 +1,9 @@
 
 DEAL::use_fast_hanging_node_algorithm = true
 DEAL::0   : 12  13  14  
-DEAL::0   : 13  2   2   1   
-DEAL::0   : 14  2   1   1   
-DEAL::0   : 13  2   1   14  
+DEAL::0   : 13  2   1   2   
+DEAL::0   : 14  1   2   1   
+DEAL::0   : 13  1   2   14  
 DEAL::0   : 0   1   2   
 DEAL::0   : 3   1   4   
 DEAL::0   : 5   4   1   
@@ -17,9 +17,9 @@ DEAL::
 DEAL::
 DEAL::use_fast_hanging_node_algorithm = false
 DEAL::0   : 12  13  14  
-DEAL::0   : 13  2   2   1   
-DEAL::0   : 14  2   1   1   
-DEAL::0   : 13  2   1   14  
+DEAL::0   : 13  2   1   2   
+DEAL::0   : 14  1   2   1   
+DEAL::0   : 13  1   2   14  
 DEAL::0   : 0   1   2   
 DEAL::0   : 3   1   4   
 DEAL::0   : 5   4   1   
@@ -27,6 +27,6 @@ DEAL::0   : 1   0   5
 DEAL::0   : 6   5   0   
 DEAL::0   : 7   8   9   10  
 DEAL::0   : 8   12  10  14  
-DEAL::0   : 9   10  3   3   1   
-DEAL::0   : 10  14  3   1   1   
+DEAL::0   : 9   10  3   1   3   
+DEAL::0   : 10  14  1   3   1   
 DEAL::