Merge pull request #14057 from kronbichler/fix_reinit_mf

Fix bug in MatrixFree::reinit for MG levels with empty process

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(),
@@ -1451,16 +1450,18 @@ namespace internal
         constexpr unsigned int max_children_per_cell =
           GeometryInfo<dim>::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<std::vector<
           std::pair<unsigned int,
                     std::array<unsigned int, max_children_per_cell>>>>
           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)

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 <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+do_test(const unsigned int n_refinements)
+{
+  parallel::distributed::Triangulation<dim> tria(
+    MPI_COMM_WORLD,
+    Triangulation<dim>::limit_level_difference_at_vertices,
+    parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(n_refinements);
+
+  FE_Q<dim>       fe(1);
+  DoFHandler<dim> 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<double> constraints;
+  constraints.close();
+
+  MappingQ1<dim> mapping;
+
+  for (unsigned int level = 0; level <= tria.n_global_levels(); ++level)
+    {
+      typename MatrixFree<dim, double>::AdditionalData data;
+      data.initialize_mapping = false;
+      if (level < tria.n_global_levels())
+        data.mg_level = level;
+
+      MatrixFree<dim, double> 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);
+}

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