From: Martin Kronbichler Date: Sun, 26 Jun 2022 18:24:05 +0000 (+0200) Subject: Fix bug in MatrixFree::reinit for MG levels with empty process X-Git-Tag: v9.5.0-rc1~1136^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c13f5ac19b05b8a79547609a51b8b9223a4b1c9c;p=dealii.git Fix bug in MatrixFree::reinit for MG levels with empty process --- diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 8ff2e3566c..f8bc962846 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1089,8 +1089,7 @@ namespace internal if (use_fast_hanging_node_algorithm) { - const auto &reference_cells = - dof_handler[0]->get_triangulation().get_reference_cells(); + const auto &reference_cells = tria.get_reference_cells(); use_fast_hanging_node_algorithm = std::all_of(reference_cells.begin(), reference_cells.end(), @@ -1463,16 +1462,18 @@ namespace internal constexpr unsigned int max_children_per_cell = GeometryInfo::max_children_per_cell; const unsigned int n_levels = - mg_level == numbers::invalid_unsigned_int ? - dof_handler[0]->get_triangulation().n_levels() - 1 : - mg_level; + mg_level == numbers::invalid_unsigned_int ? tria.n_levels() - 1 : + mg_level; std::vector>>> cell_parents(n_levels); + + // Set up data structures, making sure that the current process has + // cells on that level in the case of MG for (unsigned int level = 0; level < n_levels; ++level) - cell_parents[level].resize( - dof_handler[0]->get_triangulation().n_raw_cells(level)); + if (tria.n_levels() > level) + cell_parents[level].resize(tria.n_raw_cells(level)); for (unsigned int c = 0; c < cell_level_index_end_local; ++c) if (cell_level_index[c].first > 0) diff --git a/tests/matrix_free/reinit_matrix_free_mg.cc b/tests/matrix_free/reinit_matrix_free_mg.cc new file mode 100644 index 0000000000..169b0ef795 --- /dev/null +++ b/tests/matrix_free/reinit_matrix_free_mg.cc @@ -0,0 +1,89 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +// Verify that MatrixFree::reinit works without errors also when we have +// processes without any cells on a larger mesh (we need 17 processes to have +// a 8x8 mesh with one process not having a cell on level 1). + + +#include + +#include + +#include + +#include + +#include + +#include +#include + +#include "../tests.h" + + + +template +void +do_test(const unsigned int n_refinements) +{ + parallel::distributed::Triangulation tria( + MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + + GridGenerator::hyper_cube(tria); + tria.refine_global(n_refinements); + + FE_Q fe(1); + DoFHandler dof_handler(tria); + + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + deallog << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + AffineConstraints constraints; + constraints.close(); + + MappingQ1 mapping; + + for (unsigned int level = 0; level <= tria.n_global_levels(); ++level) + { + typename MatrixFree::AdditionalData data; + data.initialize_mapping = false; + if (level < tria.n_global_levels()) + data.mg_level = level; + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof_handler, constraints, QGauss<1>(2), data); + + deallog << "Level " << data.mg_level << " OK" << std::endl; + } +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + mpi_initlog(); + + do_test<2>(2); + do_test<2>(3); + do_test<2>(4); +} diff --git a/tests/matrix_free/reinit_matrix_free_mg.with_p4est=ture.mpirun=17.output b/tests/matrix_free/reinit_matrix_free_mg.with_p4est=ture.mpirun=17.output new file mode 100644 index 0000000000..0691792d37 --- /dev/null +++ b/tests/matrix_free/reinit_matrix_free_mg.with_p4est=ture.mpirun=17.output @@ -0,0 +1,19 @@ + +DEAL::Number of degrees of freedom: 25 +DEAL::Level 0 OK +DEAL::Level 1 OK +DEAL::Level 2 OK +DEAL::Level 4294967295 OK +DEAL::Number of degrees of freedom: 81 +DEAL::Level 0 OK +DEAL::Level 1 OK +DEAL::Level 2 OK +DEAL::Level 3 OK +DEAL::Level 4294967295 OK +DEAL::Number of degrees of freedom: 289 +DEAL::Level 0 OK +DEAL::Level 1 OK +DEAL::Level 2 OK +DEAL::Level 3 OK +DEAL::Level 4 OK +DEAL::Level 4294967295 OK