Add constrained entries to partitioner_coarse in NonNestedMG

include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h

@@ -1043,49 +1043,6 @@ namespace internal
 
   class MGTwoLevelTransferImplementation
   {
-    template <int dim, typename Number>
-    static std::shared_ptr<const Utilities::MPI::Partitioner>
-    create_coarse_partitioner(
-      const DoFHandler<dim> &                  dof_handler_coarse,
-      const dealii::AffineConstraints<Number> &constraints_coarse,
-      const unsigned int                       mg_level_coarse)
-    {
-      IndexSet locally_relevant_dofs =
-        (mg_level_coarse == numbers::invalid_unsigned_int) ?
-          DoFTools::extract_locally_active_dofs(dof_handler_coarse) :
-          DoFTools::extract_locally_active_level_dofs(dof_handler_coarse,
-                                                      mg_level_coarse);
-
-      std::vector<types::global_dof_index> locally_relevant_dofs_temp;
-
-      for (const auto i : locally_relevant_dofs)
-        {
-          if (locally_relevant_dofs.is_element(i) == false)
-            locally_relevant_dofs_temp.emplace_back(i);
-
-          const auto constraints = constraints_coarse.get_constraint_entries(i);
-
-          if (constraints)
-            for (const auto &p : *constraints)
-              if (locally_relevant_dofs.is_element(p.first) == false)
-                locally_relevant_dofs_temp.emplace_back(p.first);
-        }
-
-      std::sort(locally_relevant_dofs_temp.begin(),
-                locally_relevant_dofs_temp.end());
-      locally_relevant_dofs.add_indices(locally_relevant_dofs_temp.begin(),
-                                        locally_relevant_dofs_temp.end());
-
-      return std::make_shared<Utilities::MPI::Partitioner>(
-        mg_level_coarse == numbers::invalid_unsigned_int ?
-          dof_handler_coarse.locally_owned_dofs() :
-          dof_handler_coarse.locally_owned_mg_dofs(mg_level_coarse),
-        locally_relevant_dofs,
-        dof_handler_coarse.get_communicator());
-    }
-
-
-
     template <int dim, typename Number>
     static void
     compute_weights(
@@ -1212,6 +1169,49 @@ namespace internal
 
 
   public:
+    template <int dim, typename Number>
+    static std::shared_ptr<const Utilities::MPI::Partitioner>
+    create_coarse_partitioner(
+      const DoFHandler<dim> &                  dof_handler_coarse,
+      const dealii::AffineConstraints<Number> &constraints_coarse,
+      const unsigned int                       mg_level_coarse)
+    {
+      IndexSet locally_relevant_dofs =
+        (mg_level_coarse == numbers::invalid_unsigned_int) ?
+          DoFTools::extract_locally_active_dofs(dof_handler_coarse) :
+          DoFTools::extract_locally_active_level_dofs(dof_handler_coarse,
+                                                      mg_level_coarse);
+
+      std::vector<types::global_dof_index> locally_relevant_dofs_temp;
+
+      for (const auto i : locally_relevant_dofs)
+        {
+          if (locally_relevant_dofs.is_element(i) == false)
+            locally_relevant_dofs_temp.emplace_back(i);
+
+          const auto constraints = constraints_coarse.get_constraint_entries(i);
+
+          if (constraints)
+            for (const auto &p : *constraints)
+              if (locally_relevant_dofs.is_element(p.first) == false)
+                locally_relevant_dofs_temp.emplace_back(p.first);
+        }
+
+      std::sort(locally_relevant_dofs_temp.begin(),
+                locally_relevant_dofs_temp.end());
+      locally_relevant_dofs.add_indices(locally_relevant_dofs_temp.begin(),
+                                        locally_relevant_dofs_temp.end());
+
+      return std::make_shared<Utilities::MPI::Partitioner>(
+        mg_level_coarse == numbers::invalid_unsigned_int ?
+          dof_handler_coarse.locally_owned_dofs() :
+          dof_handler_coarse.locally_owned_mg_dofs(mg_level_coarse),
+        locally_relevant_dofs,
+        dof_handler_coarse.get_communicator());
+    }
+
+
+
     template <int dim, typename Number>
     static void
     reinit_geometric_transfer(
@@ -3839,12 +3839,9 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::
 
   // create partitioners and internal vectors
   {
-    IndexSet locally_active_dofs =
-      DoFTools::extract_locally_active_dofs(dof_handler_coarse);
-    this->partitioner_coarse.reset(
-      new Utilities::MPI::Partitioner(dof_handler_coarse.locally_owned_dofs(),
-                                      locally_active_dofs,
-                                      dof_handler_coarse.get_communicator()));
+    this->partitioner_coarse =
+      internal::MGTwoLevelTransferImplementation::create_coarse_partitioner(
+        dof_handler_coarse, constraint_coarse, numbers::invalid_unsigned_int);
 
     this->vec_coarse.reinit(this->partitioner_coarse);
   }

tests/multigrid-global-coarsening/non_nested_multigrid_04.cc

@@ -0,0 +1,209 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+/**
+ * Check that everything works also for locally-refined grids.
+ *
+ */
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+#include <deal.II/multigrid/mg_transfer_global_coarsening.h>
+
+#include "multigrid_util.h"
+
+
+template <int dim>
+void
+create_quadrant(Triangulation<dim> &tria, const unsigned int n_refinements)
+{
+  // Taken from @cite munch2022gc, according to the description in A FLEXIBLE,
+  // PARALLEL, ADAPTIVE GEOMETRIC MULTIGRID METHOD FOR FEM (Clevenger, Heister,
+  // Kanschat, Kronbichler): https://arxiv.org/pdf/1904.03317.pdf
+
+  GridGenerator::hyper_cube(tria, -1.0, +1.0);
+
+  if (n_refinements == 0)
+    return;
+
+  tria.refine_global(1);
+
+  for (unsigned int i = 1; i < n_refinements; i++)
+    {
+      for (auto cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          {
+            bool flag = true;
+            for (int d = 0; d < dim; d++)
+              if (cell->center()[d] > 0.0)
+                flag = false;
+            if (flag)
+              cell->set_refine_flag();
+          }
+      tria.execute_coarsening_and_refinement();
+    }
+
+  AssertDimension(tria.n_global_levels() - 1, n_refinements);
+}
+
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int n_refinements,
+     const unsigned int fe_degree_fine,
+     const bool         do_simplex_mesh)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  const unsigned int min_level = 0;
+  const unsigned int max_level = n_refinements;
+
+  MGLevelObject<std::shared_ptr<parallel::distributed::Triangulation<dim>>>
+                                           triangulations(min_level, max_level);
+  MGLevelObject<DoFHandler<dim>>           dof_handlers(min_level, max_level);
+  MGLevelObject<AffineConstraints<Number>> constraints(min_level, max_level);
+  MGLevelObject<std::unique_ptr<Mapping<dim>>> mappings(min_level, max_level);
+  MGLevelObject<std::shared_ptr<MGTwoLevelTransferNonNested<dim, VectorType>>>
+                                       transfers(min_level, max_level);
+  MGLevelObject<Operator<dim, Number>> operators(min_level, max_level);
+
+
+  // set up levels
+  for (auto l = min_level; l <= max_level; ++l)
+    {
+      triangulations[l] =
+        std::make_shared<parallel::distributed::Triangulation<dim>>(
+          MPI_COMM_WORLD);
+      auto &dof_handler = dof_handlers[l];
+      auto &constraint  = constraints[l];
+      auto &op          = operators[l];
+
+      std::unique_ptr<FiniteElement<dim>> fe;
+      std::unique_ptr<Quadrature<dim>>    quad;
+      std::unique_ptr<Mapping<dim>>       mapping;
+
+      if (do_simplex_mesh)
+        {
+          fe      = std::make_unique<FE_SimplexP<dim>>(fe_degree_fine);
+          quad    = std::make_unique<QGaussSimplex<dim>>(fe_degree_fine + 1);
+          mapping = std::make_unique<MappingFE<dim>>(FE_SimplexP<dim>(1));
+        }
+      else
+        {
+          fe      = std::make_unique<FE_Q<dim>>(fe_degree_fine);
+          quad    = std::make_unique<QGauss<dim>>(fe_degree_fine + 1);
+          mapping = std::make_unique<MappingQ<dim>>(1);
+        }
+
+      mappings[l] = mapping->clone();
+
+      // set up triangulation
+
+      create_quadrant(*triangulations[l], 2 * l + 1);
+
+      // set up dofhandler
+      dof_handler.reinit(*triangulations[l]);
+      dof_handler.distribute_dofs(*fe);
+
+      // set up constraints
+      IndexSet locally_relevant_dofs;
+      DoFTools::extract_locally_relevant_dofs(dof_handler,
+                                              locally_relevant_dofs);
+      constraint.reinit(locally_relevant_dofs);
+      VectorTools::interpolate_boundary_values(
+        *mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraint);
+      DoFTools::make_hanging_node_constraints(dof_handler, constraint);
+      constraint.close();
+
+      // set up operator
+      op.reinit(*mapping, dof_handler, *quad, constraint);
+    }
+
+  // set up transfer operator
+  for (unsigned int l = min_level; l < max_level; ++l)
+    {
+      transfers[l + 1] =
+        std::make_shared<MGTwoLevelTransferNonNested<dim, VectorType>>();
+      transfers[l + 1]->reinit(dof_handlers[l + 1],
+                               dof_handlers[l],
+                               *mappings[l + 1],
+                               *mappings[l],
+                               constraints[l + 1],
+                               constraints[l]);
+    }
+
+  MGTransferGlobalCoarsening<dim, VectorType> transfer(
+    transfers,
+    [&](const auto l, auto &vec) { operators[l].initialize_dof_vector(vec); });
+
+
+  GMGParameters mg_data; // TODO
+
+  VectorType dst, src;
+  operators[max_level].initialize_dof_vector(dst);
+  operators[max_level].initialize_dof_vector(src);
+
+  operators[max_level].rhs(src);
+
+  ReductionControl solver_control(
+    mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false);
+
+  mg_solve(solver_control,
+           dst,
+           src,
+           mg_data,
+           dof_handlers[max_level],
+           operators[max_level],
+           operators,
+           transfer);
+
+  deallog << dim << ' ' << fe_degree_fine << ' ' << n_refinements << ' '
+          << (do_simplex_mesh ? "tri " : "quad") << ' '
+          << solver_control.last_step() << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  deallog.precision(8);
+
+  const unsigned int n_refinements = 3;
+  const unsigned int degree        = 4;
+  test<2>(n_refinements, degree, false /*quadrilateral*/);
+}

tests/multigrid-global-coarsening/non_nested_multigrid_04.with_p4est=true.mpirun=1.output

@@ -0,0 +1,2 @@
+
+DEAL:0::2 4 3 quad 3

tests/multigrid-global-coarsening/non_nested_multigrid_04.with_p4est=true.mpirun=4.output

@@ -0,0 +1,11 @@
+
+DEAL:0::2 4 3 quad 3
+
+DEAL:1::2 4 3 quad 3
+
+
+DEAL:2::2 4 3 quad 3
+
+
+DEAL:3::2 4 3 quad 3
+