From ac097b9d2c060fed98ab6bd47bb1d67ab20edfc3 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 7 Mar 2022 22:21:08 +0100 Subject: [PATCH] Check during MF::reinit() if AffineConstraints are closed --- doc/news/changes/incompatibilities/20220308Munch | 9 +++++++++ .../deal.II/matrix_free/matrix_free.templates.h | 16 +++++++++++++--- tests/matrix_free/mixed_mesh_hn_algorithm_01.cc | 1 + .../mixed_mesh_hn_algorithm_01.output | 16 ++++++++-------- 4 files changed, 31 insertions(+), 11 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20220308Munch diff --git a/doc/news/changes/incompatibilities/20220308Munch b/doc/news/changes/incompatibilities/20220308Munch new file mode 100644 index 0000000000..b933777b05 --- /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(). +
+(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 51a6b0a120..15da51cd2e 100644 --- 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::internal_reinit( const std::shared_ptr> & mapping, const std::vector *> & dof_handler, - const std::vector *> &constraint, + const std::vector *> &constraints, const std::vector & locally_owned_dofs, const std::vector> & quad, const typename MatrixFree::AdditionalData @@ -414,7 +414,7 @@ MatrixFree::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::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::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 6a73585ea9..a9f6ae2ae1 100644 --- 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 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 384befad90..5402316913 100644 --- 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:: -- 2.39.5