first part of review

diff --git a/doc/news/changes/major/20230426FederMunchKronbichler b/doc/news/changes/major/20230426FederMunchKronbichler
new file mode 100644 (file)
index 0000000..753fcf8
--- /dev/null
+++ b/doc/news/changes/major/20230426FederMunchKronbichler
@@ -0,0 +1,3 @@
+New: Add a two level transfer operator between non-nested multigrid levels.
+<br>
+(Marco Feder, Peter Munch, Martin Kronbichler, 2023/04/26)

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index b2f32a40ba7a67e2af71303e082c767f1078d703..e02030183a83f5ab385da2235ffb72578f5bfacb 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -3488,10 +3488,10 @@ namespace internal
         constexpr unsigned int         array_size = length > 0 ? length : 1;
         std::array<Number, array_size> shape_values_x;
         std::array<Number, array_size> shape_derivs_x;
-        for (unsigned int i = 0; i < array_size; ++i)
+        for (unsigned int j = 0; j < array_size; ++j)
           {
-            shape_values_x[i] = shapes[i][0][0];
-            shape_derivs_x[i] = shapes[i][1][0];
+            shape_values_x[j] = shapes[j][0][0];
+            shape_derivs_x[j] = shapes[j][1][0];
           }
         for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
           {

diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
index 55c64c1f2245b30bc10be8a656a0683e6e68bef4..c6b2980810df08525bbc2fd0e6db06af01ac198c 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
@@ -31,6 +31,8 @@
 #include <deal.II/multigrid/mg_base.h>
 #include <deal.II/multigrid/mg_transfer_matrix_free.h>
 
+#include <deal.II/non_matching/mapping_info.h>
+
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -771,7 +773,7 @@ private:
   /**
    * MappingInfo object needed as Mapping argument by FEPointEvaluation.
    */
-  std::shared_ptr<NonMatching::MappingInfo<dim, dim>> mapping_info;
+  std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>> mapping_info;
 
   /**
    * Helper class for reading from and writing to global vectors and for
@@ -784,6 +786,12 @@ private:
    * Finite element of the coarse DoFHandler passed to reinit().
    */
   std::unique_ptr<FiniteElement<dim>> fe_coarse;
+
+  /**
+   * DoF indices of the fine cells, expressed in indices local to the MPI
+   * rank.
+   */
+  std::vector<unsigned int> level_dof_indices_fine;
 };
 
 

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 d15aa7a0a1de3bb2b1d6a0e8578a69356d86afdd..2b6115c83237f7a9b68e4728e87b0102a81997a3 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
@@ -3474,8 +3474,8 @@ namespace internal
 {
   namespace
   {
-    template <int dim>
-    std::shared_ptr<NonMatching::MappingInfo<dim>>
+    template <int dim, typename Number>
+    std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>>
     fill_mapping_info(const Utilities::MPI::RemotePointEvaluation<dim> &rpe)
     {
       const auto &cell_data = rpe.get_cell_data();
@@ -3486,10 +3486,10 @@ namespace internal
 
       for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
         {
-          typename Triangulation<dim>::active_cell_iterator cell = {
+          typename Triangulation<dim>::active_cell_iterator cell(
             &rpe.get_triangulation(),
             cell_data.cells[i].first,
-            cell_data.cells[i].second};
+            cell_data.cells[i].second);
 
           const ArrayView<const Point<dim>> unit_points(
             cell_data.reference_point_values.data() +
@@ -3497,19 +3497,14 @@ namespace internal
             cell_data.reference_point_ptrs[i + 1] -
               cell_data.reference_point_ptrs[i]);
 
-          std::vector<Point<dim>> unit_points_vector_cell;
-
-          unit_points_vector_cell.insert(unit_points_vector_cell.begin(),
-                                         unit_points.begin(),
-                                         unit_points.end());
-
           cell_iterators.emplace_back(cell);
-          unit_points_vector.emplace_back(unit_points_vector_cell);
+          unit_points_vector.emplace_back(unit_points.begin(),
+                                          unit_points.end());
         }
 
       auto mapping_info =
-        std::make_shared<NonMatching::MappingInfo<dim>>(rpe.get_mapping(),
-                                                        update_values);
+        std::make_shared<NonMatching::MappingInfo<dim, dim, Number>>(
+          rpe.get_mapping(), update_values);
       mapping_info->reinit_cells(cell_iterators, unit_points_vector);
 
       return mapping_info;
@@ -3541,12 +3536,11 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::
 
   // create partitioners and internal vectors
   {
-    IndexSet locally_relevant_dofs;
-    DoFTools::extract_locally_relevant_dofs(dof_handler_coarse,
-                                            locally_relevant_dofs);
+    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_relevant_dofs,
+                                      locally_active_dofs,
                                       dof_handler_coarse.get_communicator()));
 
     this->vec_coarse.reinit(this->partitioner_coarse);
@@ -3610,7 +3604,7 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::
   rpe.reinit(points, dof_handler_coarse.get_triangulation(), mapping_coarse);
 
   // set up MappingInfo for easier data access
-  mapping_info = internal::fill_mapping_info(rpe);
+  mapping_info = internal::fill_mapping_info<dim, Number>(rpe);
 
   // set up constraints
   const auto &cell_data = rpe.get_cell_data();
@@ -3621,11 +3615,11 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::
 
   for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
     {
-      typename DoFHandler<dim>::active_cell_iterator cell = {
+      typename DoFHandler<dim>::active_cell_iterator cell(
         &rpe.get_triangulation(),
         cell_data.cells[i].first,
         cell_data.cells[i].second,
-        &dof_handler_coarse};
+        &dof_handler_coarse);
 
       constraint_info.read_dof_indices(i,
                                        numbers::invalid_unsigned_int,
@@ -3697,23 +3691,24 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::
   if (rpe.is_map_unique() == false)
     {
       const auto evaluation_point_results_temp = evaluation_point_results;
-      evaluation_point_results.assign(rpe.get_point_ptrs().size() - 1, 0);
+      evaluation_point_results.clear();
+      evaluation_point_results.reserve(rpe.get_point_ptrs().size() - 1);
 
       const auto &ptr = rpe.get_point_ptrs();
 
       for (unsigned int i = 0; i < ptr.size() - 1; ++i)
         {
           const auto n_entries = ptr[i + 1] - ptr[i];
-          if (n_entries == 0)
-            continue;
 
           Number result = 0.0;
 
-          for (unsigned int j = 0; j < n_entries; ++j)
-            result += evaluation_point_results_temp[ptr[i] + j];
-          result /= n_entries;
-
-          evaluation_point_results[i] = result;
+          if (n_entries > 0)
+            {
+              for (unsigned int j = 0; j < n_entries; ++j)
+                result += evaluation_point_results_temp[ptr[i] + j];
+              result /= n_entries;
+            }
+          evaluation_point_results.push_back(result);
         }
     }
 
@@ -3817,7 +3812,10 @@ MGTwoLevelTransferNonNested<dim, LinearAlgebra::distributed::Vector<Number>>::
       &external_partitioner_fine)
 {
   this->internal_enable_inplace_operations_if_possible(
-    external_partitioner_coarse, external_partitioner_fine, constraint_info);
+    external_partitioner_coarse,
+    external_partitioner_fine,
+    constraint_info,
+    this->level_dof_indices_fine);
 }
 
 

diff --git a/tests/multigrid-global-coarsening/mg_transfer_util.h b/tests/multigrid-global-coarsening/mg_transfer_util.h
index 931981f9ae99254c3308e88f53639b1e0d44bb23..1031b7651fa0479b82fc525662e50d2074adbf73 100644 (file)
--- a/tests/multigrid-global-coarsening/mg_transfer_util.h
+++ b/tests/multigrid-global-coarsening/mg_transfer_util.h
@@ -208,29 +208,6 @@ test_non_nested_transfer(
     // via constraint matrix
     constraint_fine.distribute(dst);
 
-    // Visualize the result
-    {
-      DataOut<2> data_out;
-      data_out.attach_dof_handler(dof_handler_coarse);
-      data_out.add_data_vector(
-        src,
-        "coarse_solution",
-        DataOut_DoFData<dim, dim>::DataVectorType::type_dof_data);
-      data_out.build_patches();
-      std::ofstream output_coarse("coarse_sol.vtk");
-      data_out.write_vtk(output_coarse);
-      data_out.clear();
-      std::ofstream output_fine("prolonged_sol.vtk");
-      data_out.attach_dof_handler(dof_handler_fine);
-      data_out.add_data_vector(
-        dst,
-        "prolonged_solution",
-        DataOut_DoFData<dim, dim>::DataVectorType::type_dof_data);
-
-      data_out.build_patches();
-      data_out.write_vtk(output_fine);
-    }
-
     // print norms
     if (true)
       {

diff --git a/tests/multigrid-global-coarsening/non_nested_multigrid_01.output b/tests/multigrid-global-coarsening/non_nested_multigrid_01.mpirun=1.with_trilinos=true.output
similarity index 100%
rename from tests/multigrid-global-coarsening/non_nested_multigrid_01.output
rename to tests/multigrid-global-coarsening/non_nested_multigrid_01.mpirun=1.with_trilinos=true.output

diff --git a/tests/multigrid-global-coarsening/non_nested_multigrid_02.output b/tests/multigrid-global-coarsening/non_nested_multigrid_02.mpirun=1.with_trilinos=true.output
similarity index 100%
rename from tests/multigrid-global-coarsening/non_nested_multigrid_02.output
rename to tests/multigrid-global-coarsening/non_nested_multigrid_02.mpirun=1.with_trilinos=true.output

diff --git a/tests/multigrid-global-coarsening/non_nested_multigrid_03.output b/tests/multigrid-global-coarsening/non_nested_multigrid_03.mpirun=1.with_trilinos=true.output
similarity index 100%
rename from tests/multigrid-global-coarsening/non_nested_multigrid_03.output
rename to tests/multigrid-global-coarsening/non_nested_multigrid_03.mpirun=1.with_trilinos=true.output