Fix MGTransferMatrixFree::interpolate_to_mg() for PBC

diff --git a/include/deal.II/multigrid/mg_transfer_internal.h b/include/deal.II/multigrid/mg_transfer_internal.h
index 7367914349e89468381305cfca1d62b77d379624..246abdfd162a33ff92928c48b551a5a43425b945 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_internal.h
+++ b/include/deal.II/multigrid/mg_transfer_internal.h
@@ -137,6 +137,19 @@ namespace internal
       MGLevelObject<std::shared_ptr<const Utilities::MPI::Partitioner>>
         &vector_partitioners);
 
+
+
+    /**
+     * Helper function for setup_transfer. Checks for identity constrained
+     * dofs and replace with the indices of the dofs to which they are
+     * constrained
+     */
+    void
+    resolve_identity_constraints(
+      const MGConstrainedDoFs *             mg_constrained_dofs,
+      const unsigned int                    level,
+      std::vector<types::global_dof_index> &dof_indices);
+
   } // namespace MGTransfer
 } // namespace internal
 

diff --git a/include/deal.II/multigrid/mg_transfer_matrix_free.h b/include/deal.II/multigrid/mg_transfer_matrix_free.h
index 5d6324be132d8ceddeef4da6f3aaf279ee0b0429..f87fff7d480c13f710208a05ef23fe5480255133 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_matrix_free.h
+++ b/include/deal.II/multigrid/mg_transfer_matrix_free.h
@@ -536,6 +536,7 @@ MGTransferMatrixFree<dim, Number>::interpolate_to_mg(
 
   // do the transfer from level to level-1:
   dst[max_level].update_ghost_values();
+
   for (unsigned int level = max_level; level > min_level; --level)
     {
       // auxiliary vector which always has ghost elements
@@ -569,6 +570,10 @@ MGTransferMatrixFree<dim, Number>::interpolate_to_mg(
             for (unsigned int child = 0; child < cell->n_children(); ++child)
               {
                 cell->child(child)->get_mg_dof_indices(dof_indices);
+
+                internal::MGTransfer::resolve_identity_constraints(
+                  this->mg_constrained_dofs, level, dof_indices);
+
                 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
                   dof_values_fine(i) = (*input)(dof_indices[i]);
                 fe.get_restriction_matrix(child, cell->refinement_case())

diff --git a/source/multigrid/mg_transfer_internal.cc b/source/multigrid/mg_transfer_internal.cc
index 3bed763b933917df9cef2c5f7869063a6952b2ed..d3636be4a1278c248e2a9ff778b34766165097b6 100644 (file)
--- a/source/multigrid/mg_transfer_internal.cc
+++ b/source/multigrid/mg_transfer_internal.cc
@@ -641,39 +641,6 @@ namespace internal
 
 
 
-    namespace
-    {
-      /**
-       * Helper function for setup_transfer. Checks for identity constrained
-       * dofs and replace with the indices of the dofs to which they are
-       * constrained
-       */
-      void
-      replace(const MGConstrainedDoFs *             mg_constrained_dofs,
-              const unsigned int                    level,
-              std::vector<types::global_dof_index> &dof_indices)
-      {
-        if (mg_constrained_dofs != nullptr &&
-            mg_constrained_dofs->get_level_constraints(level).n_constraints() >
-              0)
-          for (auto &ind : dof_indices)
-            if (mg_constrained_dofs->get_level_constraints(level)
-                  .is_identity_constrained(ind))
-              {
-                Assert(mg_constrained_dofs->get_level_constraints(level)
-                           .get_constraint_entries(ind)
-                           ->size() == 1,
-                       ExcInternalError());
-                ind = mg_constrained_dofs->get_level_constraints(level)
-                        .get_constraint_entries(ind)
-                        ->front()
-                        .first;
-              }
-      }
-    } // namespace
-
-
-
     // Sets up most of the internal data structures of the MGTransferMatrixFree
     // class
     template <int dim, typename Number>
@@ -809,7 +776,9 @@ namespace internal
                     continue;
                   cell->child(c)->get_mg_dof_indices(local_dof_indices);
 
-                  replace(mg_constrained_dofs, level, local_dof_indices);
+                  resolve_identity_constraints(mg_constrained_dofs,
+                                               level,
+                                               local_dof_indices);
 
                   const IndexSet &owned_level_dofs =
                     dof_handler.locally_owned_mg_dofs(level);
@@ -880,7 +849,9 @@ namespace internal
                 {
                   cell->get_mg_dof_indices(local_dof_indices);
 
-                  replace(mg_constrained_dofs, level - 1, local_dof_indices);
+                  resolve_identity_constraints(mg_constrained_dofs,
+                                               level - 1,
+                                               local_dof_indices);
 
                   const IndexSet &owned_level_dofs_l0 =
                     dof_handler.locally_owned_mg_dofs(0);
@@ -1030,6 +1001,31 @@ namespace internal
         }
     }
 
+
+
+    void
+    resolve_identity_constraints(
+      const MGConstrainedDoFs *             mg_constrained_dofs,
+      const unsigned int                    level,
+      std::vector<types::global_dof_index> &dof_indices)
+    {
+      if (mg_constrained_dofs != nullptr &&
+          mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0)
+        for (auto &ind : dof_indices)
+          if (mg_constrained_dofs->get_level_constraints(level)
+                .is_identity_constrained(ind))
+            {
+              Assert(mg_constrained_dofs->get_level_constraints(level)
+                         .get_constraint_entries(ind)
+                         ->size() == 1,
+                     ExcInternalError());
+              ind = mg_constrained_dofs->get_level_constraints(level)
+                      .get_constraint_entries(ind)
+                      ->front()
+                      .first;
+            }
+    }
+
   } // namespace MGTransfer
 } // namespace internal
 

diff --git a/tests/multigrid/interpolate_to_mg_01.cc b/tests/multigrid/interpolate_to_mg_01.cc
new file mode 100644 (file)
index 0000000..9d13b93
--- /dev/null
+++ b/tests/multigrid/interpolate_to_mg_01.cc
@@ -0,0 +1,134 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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 function MGTransferMatrixFree::interpolate_to_mg() for periodic
+// boundaries
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/multigrid.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+  RightHandSideFunction(const unsigned int n_components)
+    : Function<dim>(n_components)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    return p[component % dim] * p[component % dim];
+  }
+};
+
+template <int dim>
+void
+test()
+{
+  const unsigned int fe_degree = 1;
+
+  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, -1., 1., true);
+
+  std::vector<GridTools::PeriodicFacePair<TriaIterator<CellAccessor<dim, dim>>>>
+    tria_matched_pairs;
+  GridTools::collect_periodic_faces(tria, 0, 1, 0, tria_matched_pairs);
+  tria.add_periodicity(tria_matched_pairs);
+  tria.refine_global(2);
+
+  const FE_Q<dim> fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+  dof_handler.distribute_mg_dofs();
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+  AffineConstraints<double> constraints;
+  constraints.reinit(locally_relevant_dofs);
+  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+  std::vector<
+    GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator>>
+    dof_matched_pairs;
+  GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, dof_matched_pairs);
+  DoFTools::make_periodicity_constraints<dim, dim>(dof_matched_pairs,
+                                                   constraints);
+  constraints.close();
+
+  LinearAlgebra::distributed::Vector<double> qsol(
+    dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
+
+  VectorTools::interpolate(MappingQ1<dim>(),
+                           dof_handler,
+                           RightHandSideFunction<dim>(1),
+                           qsol);
+
+  MGLevelObject<LinearAlgebra::distributed::Vector<double>> mg_qsol;
+  MGConstrainedDoFs                                         mg_constrained_dofs;
+  MGTransferMatrixFree<dim, double>                         mg_transfer;
+
+  unsigned int n_tria_levels = tria.n_global_levels();
+  mg_qsol.resize(0, n_tria_levels - 1);
+
+  mg_constrained_dofs.initialize(dof_handler);
+  mg_transfer.initialize_constraints(mg_constrained_dofs);
+  mg_transfer.build(dof_handler);
+
+  mg_transfer.interpolate_to_mg(dof_handler, mg_qsol, qsol);
+
+  for (unsigned int i = mg_qsol.min_level(); i <= mg_qsol.max_level(); ++i)
+    deallog << mg_qsol[i].l2_norm() << std::endl;
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  test<2>();
+}

diff --git a/tests/multigrid/interpolate_to_mg_01.with_p4est=true.mpirun=2.output b/tests/multigrid/interpolate_to_mg_01.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..5bc629f
--- /dev/null
+++ b/tests/multigrid/interpolate_to_mg_01.with_p4est=true.mpirun=2.output
@@ -0,0 +1,11 @@
+
+DEAL:0::2.00000
+DEAL:0::2.44949
+DEAL:0::3.25960
+DEAL:0::OK
+
+DEAL:1::2.00000
+DEAL:1::2.44949
+DEAL:1::3.25960
+DEAL:1::OK
+