From: Peter Munch <peterrmuench@gmail.com>
Date: Fri, 8 Apr 2022 12:05:29 +0000 (+0200)
Subject: Fix MGTransferGlobalCoarsening::interpolate_to_mg()
X-Git-Tag: v9.4.0-rc1~324^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=19cf04c0d3a756f6d050029e00a7b3eb5c2bd9da;p=dealii.git

Fix MGTransferGlobalCoarsening::interpolate_to_mg()
---

diff --git a/include/deal.II/matrix_free/constraint_info.h b/include/deal.II/matrix_free/constraint_info.h
index dd5c97a483..fdf56f4429 100644
--- a/include/deal.II/matrix_free/constraint_info.h
+++ b/include/deal.II/matrix_free/constraint_info.h
@@ -125,8 +125,11 @@ namespace internal
       this->plain_dof_indices_per_cell.resize(n_cells);
       this->constraint_indicator_per_cell.resize(n_cells);
 
-      if (use_fast_hanging_node_algorithm &&
-          dof_handler.get_triangulation().has_hanging_nodes())
+      // note: has_hanging_nodes() is a global operatrion
+      const bool has_hanging_nodes =
+        dof_handler.get_triangulation().has_hanging_nodes();
+
+      if (use_fast_hanging_node_algorithm && has_hanging_nodes)
         {
           hanging_nodes = std::make_unique<HangingNodes<dim>>(
             dof_handler.get_triangulation());
diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
index 62a10e09f2..3a14658609 100644
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
@@ -1481,7 +1481,8 @@ namespace internal
 
       // indices
       {
-        transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back());
+        transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back(),
+                                               numbers::invalid_unsigned_int);
 
         std::vector<types::global_dof_index> local_dof_indices(
           transfer.schemes[0].n_dofs_per_cell_coarse);
@@ -1525,7 +1526,9 @@ namespace internal
           dof_handler_coarse,
           transfer.schemes[0].n_coarse_cells +
             transfer.schemes[1].n_coarse_cells,
-          use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse));
+          constraints_coarse.n_constraints() > 0 &&
+            use_fast_hanging_node_algorithm(dof_handler_coarse,
+                                            mg_level_coarse));
 
         process_cells(
           [&](const auto &cell_coarse, const auto &cell_fine) {
@@ -1574,9 +1577,20 @@ namespace internal
               for (unsigned int i = 0;
                    i < transfer.schemes[1].n_dofs_per_cell_coarse;
                    ++i)
-                level_dof_indices_fine_1[cell_local_children_indices[c][i]] =
-                  transfer.partitioner_fine->global_to_local(
+                {
+                  const auto index = transfer.partitioner_fine->global_to_local(
                     local_dof_indices[lexicographic_numbering_fine[i]]);
+
+                  Assert(level_dof_indices_fine_1
+                               [cell_local_children_indices[c][i]] ==
+                             numbers::invalid_unsigned_int ||
+                           level_dof_indices_fine_1
+                               [cell_local_children_indices[c][i]] == index,
+                         ExcInternalError());
+
+                  level_dof_indices_fine_1[cell_local_children_indices[c][i]] =
+                    index;
+                }
             }
 
             // move pointers (only once at the end)
@@ -1653,7 +1667,7 @@ namespace internal
                     for (unsigned int j = 0; j < fe->n_dofs_per_cell(); ++j)
                       transfer.schemes[1]
                         .restriction_matrix_1d[i * n_child_dofs_1d + j +
-                                               c * shift] =
+                                               c * shift] +=
                         matrix(renumbering[i], renumbering[j]);
                 }
             }
@@ -1694,7 +1708,7 @@ namespace internal
                       transfer.schemes[1].restriction_matrix
                         [i * n_dofs_per_cell *
                            GeometryInfo<dim>::max_children_per_cell +
-                         j + c * n_dofs_per_cell] = matrix(i, j);
+                         j + c * n_dofs_per_cell] += matrix(i, j);
                 }
             }
           }
@@ -2014,7 +2028,9 @@ namespace internal
         transfer.constraint_info.reinit(
           dof_handler_coarse,
           cell_no.back(),
-          use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse));
+          constraints_coarse.n_constraints() > 0 &&
+            use_fast_hanging_node_algorithm(dof_handler_coarse,
+                                            mg_level_coarse));
 
         process_cells([&](const auto &cell_coarse, const auto &cell_fine) {
           const auto fe_pair_no =
diff --git a/tests/multigrid-global-coarsening/interpolate_01.cc b/tests/multigrid-global-coarsening/interpolate_01.cc
new file mode 100644
index 0000000000..9b9de8ec0f
--- /dev/null
+++ b/tests/multigrid-global-coarsening/interpolate_01.cc
@@ -0,0 +1,201 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+/**
+ * Test MGTransferGlobalCoarsening::interpolate_to_mg() for FE_Q and FE_DGQ.
+ */
+
+#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_system.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/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 <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+create_quadrant(Triangulation<dim> &tria, const unsigned int n_refinements)
+{
+  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>
+std::shared_ptr<Utilities::MPI::Partitioner>
+create_partitioner(const DoFHandler<dim> &dof_handler)
+{
+  return std::make_shared<Utilities::MPI::Partitioner>(
+    dof_handler.locally_owned_dofs(),
+    DoFTools::extract_locally_active_dofs(dof_handler),
+    dof_handler.get_communicator());
+}
+
+
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+  RightHandSideFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    (void)component;
+
+    return p[0];
+  }
+};
+
+
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int n_refinements,
+     const unsigned int fe_degree_fine,
+     const bool         do_dg)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  const unsigned int min_level = 0;
+  const unsigned int max_level = n_refinements;
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  create_quadrant(tria, n_refinements);
+
+  const auto trias =
+    MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence(tria);
+
+  MGLevelObject<DoFHandler<dim>> dof_handlers(min_level, max_level);
+  MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> transfers(min_level,
+                                                               max_level);
+
+  for (auto l = min_level; l <= max_level; ++l)
+    {
+      auto &dof_handler = dof_handlers[l];
+
+      dof_handler.reinit(*trias[l]);
+
+      if (do_dg)
+        dof_handler.distribute_dofs(FE_DGQ<dim>{fe_degree_fine});
+      else
+        dof_handler.distribute_dofs(FE_Q<dim>{fe_degree_fine});
+    }
+
+  for (unsigned int l = min_level; l < max_level; ++l)
+    transfers[l + 1].reinit(dof_handlers[l + 1], dof_handlers[l]);
+
+  MGTransferGlobalCoarsening<dim, VectorType> transfer(
+    transfers, [&](const auto l, auto &vec) {
+      vec.reinit(create_partitioner(dof_handlers[l]));
+    });
+
+  VectorType vec;
+  vec.reinit(create_partitioner(dof_handlers[max_level]));
+
+  RightHandSideFunction<dim> fu;
+
+  VectorTools::interpolate(dof_handlers[max_level], fu, vec);
+
+  MGLevelObject<VectorType> results(min_level, max_level);
+
+  transfer.interpolate_to_mg(dof_handlers[max_level], results, vec);
+
+  for (unsigned int l = min_level; l <= max_level; ++l)
+    {
+      results[l].update_ghost_values();
+      Vector<float> norm_per_cell(trias[l]->n_active_cells());
+      VectorTools::integrate_difference(dof_handlers[l],
+                                        results[l],
+                                        fu,
+                                        norm_per_cell,
+                                        QGauss<dim>(fe_degree_fine + 2),
+                                        VectorTools::L2_norm);
+      const double error_L2_norm =
+        VectorTools::compute_global_error(*trias[l],
+                                          norm_per_cell,
+                                          VectorTools::L2_norm);
+
+
+      deallog << l << " " << results[l].l2_norm() << " " << error_L2_norm
+              << std::endl;
+    }
+
+  deallog << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  initlog();
+
+  deallog.precision(8);
+
+  for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements)
+    for (unsigned int degree = 1; degree <= 4; ++degree)
+      {
+        test<2>(n_refinements, degree, true);
+        test<2>(n_refinements, degree, false);
+      }
+}
diff --git a/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output
new file mode 100644
index 0000000000..ec37f09d9c
--- /dev/null
+++ b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=1.output
@@ -0,0 +1,121 @@
+
+DEAL::0 2.0000000 8.1041682e-16
+DEAL::1 2.8284271 3.6480780e-16
+DEAL::2 3.4641016 0.0000000
+DEAL::
+DEAL::0 2.0000000 0.0000000
+DEAL::1 2.4494897 0.0000000
+DEAL::2 2.7838822 0.0000000
+DEAL::
+DEAL::0 2.4494897 4.3469376e-16
+DEAL::1 3.8729833 1.9846630e-16
+DEAL::2 4.8989795 1.2272899e-16
+DEAL::
+DEAL::0 2.4494897 1.1453553e-16
+DEAL::1 3.5355339 9.9806058e-17
+DEAL::2 4.2866070 1.1488487e-16
+DEAL::
+DEAL::0 3.0983867 6.6306930e-16
+DEAL::1 5.0596443 2.7087997e-16
+DEAL::2 6.4498062 2.3247867e-16
+DEAL::
+DEAL::0 3.0983867 3.0349490e-16
+DEAL::1 4.7328638 2.0099536e-16
+DEAL::2 5.8694122 2.1659599e-16
+DEAL::
+DEAL::0 3.7796447 1.1027361e-15
+DEAL::1 6.2678317 3.9148216e-16
+DEAL::2 8.0178373 3.1129099e-16
+DEAL::
+DEAL::0 3.7796447 2.3002617e-16
+DEAL::1 5.9461873 2.9830841e-16
+DEAL::2 7.4558223 3.2071973e-16
+DEAL::
+DEAL::0 2.0000000 1.1994317e-15
+DEAL::1 2.8284271 8.0596747e-16
+DEAL::2 3.4641016 6.0282265e-16
+DEAL::3 6.3245553 0.0000000
+DEAL::
+DEAL::0 2.0000000 0.0000000
+DEAL::1 2.4494897 0.0000000
+DEAL::2 2.7838822 0.0000000
+DEAL::3 4.2573466 0.0000000
+DEAL::
+DEAL::0 2.4494897 6.6147021e-16
+DEAL::1 3.8729833 4.1436381e-16
+DEAL::2 4.8989795 3.4330648e-16
+DEAL::3 9.3273791 1.4568899e-16
+DEAL::
+DEAL::0 2.4494897 1.1453553e-16
+DEAL::1 3.5355339 9.9806058e-17
+DEAL::2 4.2866070 1.1488487e-16
+DEAL::3 7.3271754 1.3584636e-16
+DEAL::
+DEAL::0 3.0983867 9.8309852e-16
+DEAL::1 5.0596443 5.6517728e-16
+DEAL::2 6.4498062 3.8488399e-16
+DEAL::3 12.393547 2.3245655e-16
+DEAL::
+DEAL::0 3.0983867 4.3584763e-16
+DEAL::1 4.7328638 2.7306050e-16
+DEAL::2 5.8694122 2.3526684e-16
+DEAL::3 10.419933 2.2576644e-16
+DEAL::
+DEAL::0 3.7796447 1.4273981e-15
+DEAL::1 6.2678317 7.7585762e-16
+DEAL::2 8.0178373 7.4056572e-16
+DEAL::3 15.468863 3.7969912e-16
+DEAL::
+DEAL::0 3.7796447 2.8778492e-16
+DEAL::1 5.9461873 3.3252973e-16
+DEAL::2 7.4558223 3.4411320e-16
+DEAL::3 13.509587 3.6867808e-16
+DEAL::
+DEAL::0 2.0000000 1.4131724e-15
+DEAL::1 2.8284271 1.0606982e-15
+DEAL::2 3.4641016 8.5294693e-16
+DEAL::3 6.3245553 4.3812845e-16
+DEAL::4 10.723805 0.0000000
+DEAL::
+DEAL::0 2.0000000 0.0000000
+DEAL::1 2.4494897 0.0000000
+DEAL::2 2.7838822 0.0000000
+DEAL::3 4.2573466 0.0000000
+DEAL::4 6.4128777 0.0000000
+DEAL::
+DEAL::0 2.4494897 8.1236828e-16
+DEAL::1 3.8729833 5.3562441e-16
+DEAL::2 4.8989795 4.2374356e-16
+DEAL::3 9.3273791 2.4185602e-16
+DEAL::4 15.992186 1.4728562e-16
+DEAL::
+DEAL::0 2.4494897 1.1453553e-16
+DEAL::1 3.5355339 9.9806058e-17
+DEAL::2 4.2866070 1.1488487e-16
+DEAL::3 7.3271754 1.3584636e-16
+DEAL::4 11.805719 1.3167380e-16
+DEAL::
+DEAL::0 3.0983867 1.4470137e-15
+DEAL::1 5.0596443 8.9797785e-16
+DEAL::2 6.4498062 6.6972640e-16
+DEAL::3 12.393547 3.5145888e-16
+DEAL::4 21.297887 2.3268951e-16
+DEAL::
+DEAL::0 3.0983867 6.3814779e-16
+DEAL::1 4.7328638 4.5675034e-16
+DEAL::2 5.8694122 3.7189323e-16
+DEAL::3 10.419933 2.4820397e-16
+DEAL::4 17.160274 2.2156718e-16
+DEAL::
+DEAL::0 3.7796447 1.4543742e-15
+DEAL::1 6.2678317 8.5375061e-16
+DEAL::2 8.0178373 8.4327874e-16
+DEAL::3 15.468863 5.4672827e-16
+DEAL::4 26.608940 3.8089812e-16
+DEAL::
+DEAL::0 3.7796447 3.6858565e-16
+DEAL::1 5.9461873 3.5838440e-16
+DEAL::2 7.4558223 3.6691834e-16
+DEAL::3 13.509587 3.6351703e-16
+DEAL::4 22.497222 3.6587400e-16
+DEAL::
diff --git a/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=4.output b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=4.output
new file mode 100644
index 0000000000..ec37f09d9c
--- /dev/null
+++ b/tests/multigrid-global-coarsening/interpolate_01.with_p4est=true.mpirun=4.output
@@ -0,0 +1,121 @@
+
+DEAL::0 2.0000000 8.1041682e-16
+DEAL::1 2.8284271 3.6480780e-16
+DEAL::2 3.4641016 0.0000000
+DEAL::
+DEAL::0 2.0000000 0.0000000
+DEAL::1 2.4494897 0.0000000
+DEAL::2 2.7838822 0.0000000
+DEAL::
+DEAL::0 2.4494897 4.3469376e-16
+DEAL::1 3.8729833 1.9846630e-16
+DEAL::2 4.8989795 1.2272899e-16
+DEAL::
+DEAL::0 2.4494897 1.1453553e-16
+DEAL::1 3.5355339 9.9806058e-17
+DEAL::2 4.2866070 1.1488487e-16
+DEAL::
+DEAL::0 3.0983867 6.6306930e-16
+DEAL::1 5.0596443 2.7087997e-16
+DEAL::2 6.4498062 2.3247867e-16
+DEAL::
+DEAL::0 3.0983867 3.0349490e-16
+DEAL::1 4.7328638 2.0099536e-16
+DEAL::2 5.8694122 2.1659599e-16
+DEAL::
+DEAL::0 3.7796447 1.1027361e-15
+DEAL::1 6.2678317 3.9148216e-16
+DEAL::2 8.0178373 3.1129099e-16
+DEAL::
+DEAL::0 3.7796447 2.3002617e-16
+DEAL::1 5.9461873 2.9830841e-16
+DEAL::2 7.4558223 3.2071973e-16
+DEAL::
+DEAL::0 2.0000000 1.1994317e-15
+DEAL::1 2.8284271 8.0596747e-16
+DEAL::2 3.4641016 6.0282265e-16
+DEAL::3 6.3245553 0.0000000
+DEAL::
+DEAL::0 2.0000000 0.0000000
+DEAL::1 2.4494897 0.0000000
+DEAL::2 2.7838822 0.0000000
+DEAL::3 4.2573466 0.0000000
+DEAL::
+DEAL::0 2.4494897 6.6147021e-16
+DEAL::1 3.8729833 4.1436381e-16
+DEAL::2 4.8989795 3.4330648e-16
+DEAL::3 9.3273791 1.4568899e-16
+DEAL::
+DEAL::0 2.4494897 1.1453553e-16
+DEAL::1 3.5355339 9.9806058e-17
+DEAL::2 4.2866070 1.1488487e-16
+DEAL::3 7.3271754 1.3584636e-16
+DEAL::
+DEAL::0 3.0983867 9.8309852e-16
+DEAL::1 5.0596443 5.6517728e-16
+DEAL::2 6.4498062 3.8488399e-16
+DEAL::3 12.393547 2.3245655e-16
+DEAL::
+DEAL::0 3.0983867 4.3584763e-16
+DEAL::1 4.7328638 2.7306050e-16
+DEAL::2 5.8694122 2.3526684e-16
+DEAL::3 10.419933 2.2576644e-16
+DEAL::
+DEAL::0 3.7796447 1.4273981e-15
+DEAL::1 6.2678317 7.7585762e-16
+DEAL::2 8.0178373 7.4056572e-16
+DEAL::3 15.468863 3.7969912e-16
+DEAL::
+DEAL::0 3.7796447 2.8778492e-16
+DEAL::1 5.9461873 3.3252973e-16
+DEAL::2 7.4558223 3.4411320e-16
+DEAL::3 13.509587 3.6867808e-16
+DEAL::
+DEAL::0 2.0000000 1.4131724e-15
+DEAL::1 2.8284271 1.0606982e-15
+DEAL::2 3.4641016 8.5294693e-16
+DEAL::3 6.3245553 4.3812845e-16
+DEAL::4 10.723805 0.0000000
+DEAL::
+DEAL::0 2.0000000 0.0000000
+DEAL::1 2.4494897 0.0000000
+DEAL::2 2.7838822 0.0000000
+DEAL::3 4.2573466 0.0000000
+DEAL::4 6.4128777 0.0000000
+DEAL::
+DEAL::0 2.4494897 8.1236828e-16
+DEAL::1 3.8729833 5.3562441e-16
+DEAL::2 4.8989795 4.2374356e-16
+DEAL::3 9.3273791 2.4185602e-16
+DEAL::4 15.992186 1.4728562e-16
+DEAL::
+DEAL::0 2.4494897 1.1453553e-16
+DEAL::1 3.5355339 9.9806058e-17
+DEAL::2 4.2866070 1.1488487e-16
+DEAL::3 7.3271754 1.3584636e-16
+DEAL::4 11.805719 1.3167380e-16
+DEAL::
+DEAL::0 3.0983867 1.4470137e-15
+DEAL::1 5.0596443 8.9797785e-16
+DEAL::2 6.4498062 6.6972640e-16
+DEAL::3 12.393547 3.5145888e-16
+DEAL::4 21.297887 2.3268951e-16
+DEAL::
+DEAL::0 3.0983867 6.3814779e-16
+DEAL::1 4.7328638 4.5675034e-16
+DEAL::2 5.8694122 3.7189323e-16
+DEAL::3 10.419933 2.4820397e-16
+DEAL::4 17.160274 2.2156718e-16
+DEAL::
+DEAL::0 3.7796447 1.4543742e-15
+DEAL::1 6.2678317 8.5375061e-16
+DEAL::2 8.0178373 8.4327874e-16
+DEAL::3 15.468863 5.4672827e-16
+DEAL::4 26.608940 3.8089812e-16
+DEAL::
+DEAL::0 3.7796447 3.6858565e-16
+DEAL::1 5.9461873 3.5838440e-16
+DEAL::2 7.4558223 3.6691834e-16
+DEAL::3 13.509587 3.6351703e-16
+DEAL::4 22.497222 3.6587400e-16
+DEAL::