Final implementation + update documentation and test case

Co-authored-by: Doug Shi-Dong <doug.shidong@gmail.com>

diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h
index 341e3429c532e406b28e4a77c578f0e33d5430eb..fafa5d3d7b40b6de66a396687c941fc25de60649 100644 (file)
--- a/include/deal.II/numerics/vector_tools.h
+++ b/include/deal.II/numerics/vector_tools.h
@@ -641,7 +641,7 @@ namespace VectorTools
    * Interpolate different finite element spaces. The interpolation of vector
    * @p data_1 (which is assumed to be ghosted, see @ref GlossGhostedVector)
    * is executed from the FE space represented by @p dof_1
-   * to the vector @p data_2 (which must be non-ghosted) on FE space @p dof_2.
+   * to the vector @p data_2 on FE space @p dof_2.
    * The interpolation on each cell is represented by the matrix @p transfer.
    * Curved boundaries are neglected so far.
    *

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 65316239e94e350b26e1e84a4751510c4feb7df8..c95d4b7b7ea8bd02ebbafc3115b00a5fa0e7b970 100644 (file)
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -602,6 +602,9 @@ namespace VectorTools
     Vector<number> cell_data_1(dof_1.get_fe().dofs_per_cell);
     Vector<number> cell_data_2(dof_2.get_fe().dofs_per_cell);
 
+    // Reset output vector.
+    data_2 = static_cast<number>(0);
+
     // Store how many cells share each dof (unghosted).
     OutVector touch_count;
     touch_count.reinit(data_2);
@@ -660,11 +663,12 @@ namespace VectorTools
         ::dealii::internal::ElementAccess<OutVector>::set(value, i, data_2);
       }
 
-    // Compress data_2 to set the proper values on all the processors.
+    // Compress data_2 to set the proper values on all the parallel processes.
     data_2.compress(VectorOperation::insert);
   }
 
 
+
   template <int dim,
             int spacedim,
             typename VectorType,

diff --git a/tests/vector_tools/interpolate.cc b/tests/vector_tools/interpolate.cc
index c52b5261dcf4db514b4cc24a1b7a2844b1fbba73..f0f86ade787032ab75cc1405a54ddc677e40a643 100644 (file)
--- a/tests/vector_tools/interpolate.cc
+++ b/tests/vector_tools/interpolate.cc
@@ -1,3 +1,28 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2019 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.
+//
+// ---------------------------------------------------------------------
+
+// The test initializes a non-zero solution over a hypercube using the
+// finite element space FESystem(FE_Q(2), dim). The solution is then
+// interpolated to the finite element space FESystem(FE_Q(3), dim) using
+// the parallel VectorTools::interpolate() function.
+//
+// Since the solution is continuous across the elements, the solution
+// values at the new support points of FESystem(FE_Q(3), dim) should be
+// equal to the old solution defined in FESystem(FE_Q(2), dim) at those
+// same support points.
+
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/distributed/tria.h>
@@ -16,134 +41,11 @@
 
 template <int dim>
 void
-run_1d()
-{
-  using Vector = dealii::Vector<double>;
-
-  // Create grid
-  dealii::Triangulation<dim> triangulation(
-    typename dealii::Triangulation<dim>::MeshSmoothing(
-      dealii::Triangulation<dim>::smoothing_on_refinement |
-      dealii::Triangulation<dim>::smoothing_on_coarsening));
-  const unsigned int n_1d_cells = 2;
-  dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells);
-
-  // Create low order system
-  dealii::DoFHandler<dim> dof_handler_coarse(triangulation);
-
-  const unsigned int          degree_coarse = 2;
-  const dealii::FE_Q<dim>     fe_q_coarse(degree_coarse);
-  const dealii::FESystem<dim> fe_system_coarse(fe_q_coarse, dim);
-  dof_handler_coarse.initialize(triangulation, fe_system_coarse);
-  dof_handler_coarse.distribute_dofs(fe_system_coarse);
-
-  Vector solution_coarse;
-  solution_coarse.reinit(dof_handler_coarse.n_dofs());
-  // Initialize dummy coarse solution
-  // solution_coarse.add(1.0);
-  for (unsigned int idof = 0; idof < dof_handler_coarse.n_dofs(); ++idof)
-    {
-      solution_coarse[idof] = idof;
-    }
-
-  const unsigned int degree_fine = 3;
-
-  dealii::DoFHandler<dim> dof_handler_fine(triangulation);
-
-  const dealii::FE_Q<dim>     fe_q_fine(degree_fine);
-  const dealii::FESystem<dim> fe_system_fine(fe_q_fine, dim);
-  dof_handler_fine.initialize(triangulation, fe_system_fine);
-  dof_handler_fine.distribute_dofs(fe_system_fine);
-
-  Vector solution_fine;
-  solution_fine.reinit(dof_handler_fine.n_dofs());
-  solution_fine.add(1.0); // Initialize to non-zero
-
-  // Interpolate the solution
-  dealii::FullMatrix<double> interpolation_matrix(
-    fe_system_fine.dofs_per_cell, fe_system_coarse.dofs_per_cell);
-  fe_system_fine.get_interpolation_matrix(fe_system_coarse,
-                                          interpolation_matrix);
-  dealii::VectorTools::interpolate(dof_handler_coarse,
-                                   dof_handler_fine,
-                                   interpolation_matrix,
-                                   solution_coarse,
-                                   solution_fine);
-
-
-  const unsigned int n_dofs_coarse = fe_system_coarse.dofs_per_cell;
-  const unsigned int n_dofs_fine   = fe_system_fine.dofs_per_cell;
-  const std::vector<dealii::Point<dim>> &points =
-    fe_system_fine.get_unit_support_points();
-
-  // Check that interpolated solution matches at the support points
-  std::vector<dealii::types::global_dof_index> dof_indices_coarse(
-    n_dofs_coarse);
-  std::vector<dealii::types::global_dof_index> dof_indices_fine(n_dofs_fine);
-
-  auto cell_coarse = dof_handler_coarse.begin_active();
-  auto cell_fine   = dof_handler_fine.begin_active();
-  auto endcell     = dof_handler_fine.end();
-
-  bool matching = true;
-  for (; cell_fine != dof_handler_fine.end(); ++cell_fine, ++cell_coarse)
-    {
-      if (!cell_fine->is_locally_owned())
-        continue;
-      cell_fine->get_dof_indices(dof_indices_fine);
-      cell_coarse->get_dof_indices(dof_indices_coarse);
-      for (unsigned int iquad = 0; iquad < points.size(); ++iquad)
-        {
-          dealii::Tensor<1, dim, double> local_solution_coarse;
-          dealii::Tensor<1, dim, double> local_solution_fine;
-          for (unsigned int idof = 0; idof < n_dofs_coarse; ++idof)
-            {
-              const unsigned int axis =
-                fe_system_coarse.system_to_component_index(idof).first;
-              local_solution_coarse[axis] +=
-                solution_coarse[dof_indices_coarse[idof]] *
-                fe_system_coarse.shape_value(idof, points[iquad]);
-            }
-          for (unsigned int idof = 0; idof < n_dofs_fine; ++idof)
-            {
-              const unsigned int axis =
-                fe_system_fine.system_to_component_index(idof).first;
-              local_solution_fine[axis] +=
-                solution_fine[dof_indices_fine[idof]] *
-                fe_system_fine.shape_value(idof, points[iquad]);
-            }
-          dealii::Tensor<1, dim, double> diff = local_solution_fine;
-          diff -= local_solution_coarse;
-          const double diff_norm = diff.norm();
-
-          std::cout << "Fine solution " << local_solution_fine
-                    << " Coarse solution " << local_solution_coarse
-                    << " Difference " << diff_norm << std::endl;
-          if (diff_norm > 1e-12)
-            matching = false;
-        }
-    }
-  if (matching)
-    dealii::deallog << "OK" << std::endl;
-  if (!matching)
-    dealii::deallog << "Non matching" << std::endl;
-}
-
-template <int dim>
-void
-run()
+run(const dealii::Triangulation<dim> &triangulation)
 {
   using Vector = dealii::LinearAlgebra::distributed::Vector<double>;
-  // Create grid
-  dealii::parallel::distributed::Triangulation<dim> triangulation(
-    MPI_COMM_WORLD,
-    typename dealii::Triangulation<dim>::MeshSmoothing(
-      dealii::Triangulation<dim>::smoothing_on_refinement |
-      dealii::Triangulation<dim>::smoothing_on_coarsening));
-  const unsigned int n_1d_cells = 2;
-  dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells);
 
-  // Create low order system
+  // Create low order system.
   dealii::DoFHandler<dim> dof_handler_coarse(triangulation);
 
   const unsigned int          degree_coarse = 2;
@@ -160,10 +62,19 @@ run()
                                                   ghost_dofs_coarse);
   dealii::IndexSet locally_relevant_dofs_coarse = ghost_dofs_coarse;
   ghost_dofs_coarse.subtract_set(locally_owned_dofs_coarse);
-  solution_coarse.reinit(locally_owned_dofs_coarse,
-                         ghost_dofs_coarse,
-                         MPI_COMM_WORLD);
-  // Initialize dummy coarse solution
+
+  if (dim == 1)
+    {
+      solution_coarse.reinit(locally_owned_dofs_coarse, MPI_COMM_WORLD);
+    }
+  else
+    {
+      solution_coarse.reinit(locally_owned_dofs_coarse,
+                             ghost_dofs_coarse,
+                             MPI_COMM_WORLD);
+    }
+
+  // Initialize dummy coarse solution.
   for (unsigned int idof = 0; idof < dof_handler_coarse.n_dofs(); ++idof)
     {
       if (locally_owned_dofs_coarse.is_element(idof))
@@ -190,13 +101,19 @@ run()
                                                   ghost_dofs_fine);
   dealii::IndexSet locally_relevant_dofs_fine = ghost_dofs_fine;
   ghost_dofs_fine.subtract_set(locally_owned_dofs_fine);
-  solution_fine.reinit(locally_owned_dofs_fine,
-                       ghost_dofs_fine,
-                       MPI_COMM_WORLD);
-  solution_fine.add(1.0); // Initialize to non-zero
-  solution_fine.update_ghost_values();
 
-  // Interpolate the solution
+  if (dim == 1)
+    {
+      solution_fine.reinit(locally_owned_dofs_fine, MPI_COMM_WORLD);
+    }
+  else
+    {
+      solution_fine.reinit(locally_owned_dofs_fine,
+                           ghost_dofs_fine,
+                           MPI_COMM_WORLD);
+    }
+
+  // Interpolate the solution.
   dealii::FullMatrix<double> interpolation_matrix(
     fe_system_fine.dofs_per_cell, fe_system_coarse.dofs_per_cell);
   fe_system_fine.get_interpolation_matrix(fe_system_coarse,
@@ -206,14 +123,14 @@ run()
                                    interpolation_matrix,
                                    solution_coarse,
                                    solution_fine);
-
+  solution_fine.update_ghost_values();
 
   const unsigned int n_dofs_coarse = fe_system_coarse.dofs_per_cell;
   const unsigned int n_dofs_fine   = fe_system_fine.dofs_per_cell;
   const std::vector<dealii::Point<dim>> &points =
     fe_system_fine.get_unit_support_points();
 
-  // Check that interpolated solution matches at the support points
+  // Check that interpolated solution matches at quadrature points.
   std::vector<dealii::types::global_dof_index> dof_indices_coarse(
     n_dofs_coarse);
   std::vector<dealii::types::global_dof_index> dof_indices_fine(n_dofs_fine);
@@ -223,16 +140,19 @@ run()
   auto endcell     = dof_handler_fine.end();
 
   bool matching = true;
-  for (; cell_fine != dof_handler_fine.end(); ++cell_fine, ++cell_coarse)
+  for (; cell_fine != endcell; ++cell_fine, ++cell_coarse)
     {
       if (!cell_fine->is_locally_owned())
         continue;
+
       cell_fine->get_dof_indices(dof_indices_fine);
       cell_coarse->get_dof_indices(dof_indices_coarse);
+
       for (unsigned int iquad = 0; iquad < points.size(); ++iquad)
         {
           dealii::Tensor<1, dim, double> local_solution_coarse;
           dealii::Tensor<1, dim, double> local_solution_fine;
+
           for (unsigned int idof = 0; idof < n_dofs_coarse; ++idof)
             {
               const unsigned int axis =
@@ -241,6 +161,7 @@ run()
                 solution_coarse[dof_indices_coarse[idof]] *
                 fe_system_coarse.shape_value(idof, points[iquad]);
             }
+
           for (unsigned int idof = 0; idof < n_dofs_fine; ++idof)
             {
               const unsigned int axis =
@@ -249,21 +170,31 @@ run()
                 solution_fine[dof_indices_fine[idof]] *
                 fe_system_fine.shape_value(idof, points[iquad]);
             }
+
           dealii::Tensor<1, dim, double> diff = local_solution_fine;
           diff -= local_solution_coarse;
+
           const double diff_norm = diff.norm();
 
-          std::cout << "Fine solution " << local_solution_fine
-                    << " Coarse solution " << local_solution_coarse
-                    << " Difference " << diff_norm << std::endl;
+          std::cout << "Coarse solution = " << local_solution_coarse
+                    << ", fine solution = " << local_solution_fine
+                    << ", difference = " << diff_norm << std::endl;
+
           if (diff_norm > 1e-12)
-            matching = false;
+            {
+              matching = false;
+            }
         }
     }
+
   if (matching)
-    dealii::deallog << "OK" << std::endl;
-  if (!matching)
-    dealii::deallog << "Non matching" << std::endl;
+    {
+      dealii::deallog << "OK" << std::endl;
+    }
+  else
+    {
+      dealii::deallog << "Non matching" << std::endl;
+    }
 }
 
 int
@@ -272,11 +203,34 @@ main(int argc, char *argv[])
   dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
   initlog();
 
-  // dealii::deallog.push("1d");
-  dealii::deallog << "1d" << std::endl;
-  run_1d<1>();
-  dealii::deallog << "2d" << std::endl;
-  run<2>();
-  dealii::deallog << "3d" << std::endl;
-  run<3>();
+  const unsigned int n_1d_cells = 2;
+
+  {
+    // 1D
+    dealii::Triangulation<1> triangulation;
+    dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells);
+
+    dealii::deallog << "1D" << std::endl;
+    run<1>(triangulation);
+  }
+
+  // 2D.
+  {
+    dealii::parallel::distributed::Triangulation<2> triangulation(
+      MPI_COMM_WORLD);
+    dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells);
+
+    dealii::deallog << "2D" << std::endl;
+    run<2>(triangulation);
+  }
+
+  // 3D.
+  {
+    dealii::parallel::distributed::Triangulation<3> triangulation(
+      MPI_COMM_WORLD);
+    dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells);
+
+    dealii::deallog << "3D" << std::endl;
+    run<3>(triangulation);
+  }
 }

diff --git a/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output b/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output
index 4b883267c3ce38cfa0d04c1ac93f446a59cdacd0..661e4988ab67a8db3c3d6ba3b78bb3d14a10ceab 100644 (file)
--- a/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output
+++ b/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output
@@ -1,7 +1,7 @@
 
-DEAL::1d
+DEAL::1D
 DEAL::OK
-DEAL::2d
+DEAL::2D
 DEAL::OK
-DEAL::3d
+DEAL::3D
 DEAL::OK

diff --git a/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output b/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output
index 4b883267c3ce38cfa0d04c1ac93f446a59cdacd0..661e4988ab67a8db3c3d6ba3b78bb3d14a10ceab 100644 (file)
--- a/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output
+++ b/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output
@@ -1,7 +1,7 @@
 
-DEAL::1d
+DEAL::1D
 DEAL::OK
-DEAL::2d
+DEAL::2D
 DEAL::OK
-DEAL::3d
+DEAL::3D
 DEAL::OK