Implement optimized copy for generic vectors

diff --git a/include/deal.II/multigrid/mg_transfer.h b/include/deal.II/multigrid/mg_transfer.h
index 2dc701135f89acc5a5f7a46d9bee7544427fa624..b5316b6e68c33c002f2e48f0e3bc432b09e257af 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer.h
+++ b/include/deal.II/multigrid/mg_transfer.h
@@ -233,6 +233,14 @@ protected:
   std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > >
   copy_indices_level_mine;
 
+  /**
+   * Stores whether the copy operation from the global to the level vector is
+   * actually a plain copy to the finest level. This means that the grid has
+   * no adaptive refinement and the numbering on the finest multigrid level is
+   * the same as in the global case.
+   */
+  bool perform_plain_copy;
+
   /**
    * The vector that stores what has been given to the
    * set_component_to_block_map() function.

diff --git a/include/deal.II/multigrid/mg_transfer.templates.h b/include/deal.II/multigrid/mg_transfer.templates.h
index a3bd668bdf4aaa2d6514f1bd64ace334d0412cfc..59bb8f241382764b042e8b54a9fba43ce3d39c79 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer.templates.h
+++ b/include/deal.II/multigrid/mg_transfer.templates.h
@@ -161,6 +161,34 @@ namespace
 /* ------------------ MGLevelGlobalTransfer<VectorType> ----------------- */
 
 
+namespace internal
+{
+  // generic copy function of two different vectors -> need to access each
+  // individual entry
+  template <typename T, typename V>
+  void copy_vector (const std::vector<std::pair<types::global_dof_index,types::global_dof_index> > &copy_indices,
+                    const T &src,
+                    V &dst)
+  {
+    // we should have i->second == i->first, therefore we can use the same
+    // function for both copying to mg as well as copying from mg
+    for (std::vector<std::pair<types::global_dof_index, types::global_dof_index> >::
+         const_iterator i = copy_indices.begin(); i != copy_indices.end(); ++i)
+      dst(i->first) = src(i->first);
+    dst.compress(VectorOperation::insert);
+  }
+
+  // specialized copy function for the same vector
+  template <typename T>
+  void copy_vector (const std::vector<std::pair<types::global_dof_index,types::global_dof_index> > &,
+                    const T &src,
+                    T &dst)
+  {
+    dst = src;
+  }
+}
+
+
 template <typename VectorType>
 template <int dim, class InVector, int spacedim>
 void
@@ -175,6 +203,24 @@ MGLevelGlobalTransfer<VectorType>::copy_to_mg
   std::cout << "copy_to_mg src " << src.l2_norm() << std::endl;
   MPI_Barrier(MPI_COMM_WORLD);
 #endif
+
+  if (perform_plain_copy)
+    {
+      // if the finest multigrid level covers the whole domain (i.e., no
+      // adaptive refinement) and the numbering of the finest level DoFs and
+      // the global DoFs are the same, we can do a plain copy
+      AssertDimension(dst[dst.max_level()].size(), src.size());
+      internal::copy_vector(copy_indices[dst.max_level()], src, dst[dst.max_level()]);
+
+      // do the initial restriction
+      for (unsigned int level=mg_dof_handler.get_triangulation().n_global_levels()-1; level != 0; )
+        {
+          --level;
+          this->restrict_and_add (level+1, dst[level], dst[level+1]);
+        }
+      return;
+    }
+
   for (unsigned int level=mg_dof_handler.get_triangulation().n_global_levels(); level != 0;)
     {
       --level;
@@ -225,6 +271,13 @@ MGLevelGlobalTransfer<VectorType>::copy_from_mg
  OutVector                       &dst,
  const MGLevelObject<VectorType> &src) const
 {
+  if (perform_plain_copy)
+    {
+      AssertDimension(dst.size(), src[src.max_level()].size());
+      internal::copy_vector(copy_indices[src.max_level()], src[src.max_level()], dst);
+      return;
+    }
+
   // For non-DG: degrees of freedom in the refinement face may need special
   // attention, since they belong to the coarse level, but have fine level
   // basis functions
@@ -331,6 +384,8 @@ MGLevelGlobalTransfer<parallel::distributed::Vector<Number> >::copy_to_mg
       VectorView<Number>  dst_view (src.local_size(), dst[dst.max_level()].begin());
       VectorView<Number2> src_view (src.local_size(), src.begin());
       static_cast<Vector<Number> &>(dst_view) = static_cast<Vector<Number2> &>(src_view);
+
+      // do the initial restriction
       for (unsigned int level=mg_dof_handler.get_triangulation().n_global_levels()-1; level != 0; )
         {
           --level;
@@ -387,10 +442,6 @@ MGLevelGlobalTransfer<parallel::distributed::Vector<Number> >::copy_from_mg
  parallel::distributed::Vector<Number2>                      &dst,
  const MGLevelObject<parallel::distributed::Vector<Number> > &src) const
 {
-  // For non-DG: degrees of freedom in the refinement face may need special
-  // attention, since they belong to the coarse level, but have fine level
-  // basis functions
-
   if (perform_plain_copy)
     {
       // In this case, we can simply copy the local range (in parallel by
@@ -404,6 +455,10 @@ MGLevelGlobalTransfer<parallel::distributed::Vector<Number> >::copy_from_mg
       return;
     }
 
+  // For non-DG: degrees of freedom in the refinement face may need special
+  // attention, since they belong to the coarse level, but have fine level
+  // basis functions
+
   dst = 0;
   for (unsigned int level=0; level<mg_dof_handler.get_triangulation().n_global_levels(); ++level)
     {

diff --git a/source/multigrid/mg_level_global_transfer.cc b/source/multigrid/mg_level_global_transfer.cc
index 8e5b6b55b5754d00a9fcf5b48598e20f61bd5a94..8fec5c832c959607daf2c1d05e9ee2a6fd77d594 100644 (file)
--- a/source/multigrid/mg_level_global_transfer.cc
+++ b/source/multigrid/mg_level_global_transfer.cc
@@ -271,6 +271,38 @@ MGLevelGlobalTransfer<VectorType>::fill_and_communicate_copy_indices
 {
   fill_copy_indices(mg_dof, mg_constrained_dofs, copy_indices,
                     copy_indices_global_mine, copy_indices_level_mine);
+
+  // check if we can run a plain copy operation between the global DoFs and
+  // the finest level.
+  perform_plain_copy =
+    (copy_indices.back().size() == mg_dof.locally_owned_dofs().n_elements())
+    &&
+    (mg_dof.locally_owned_dofs().n_elements() ==
+     mg_dof.locally_owned_mg_dofs(mg_dof.get_triangulation().n_global_levels()-1).n_elements());
+  if (perform_plain_copy)
+    {
+      AssertDimension(copy_indices_global_mine.back().size(), 0);
+      AssertDimension(copy_indices_level_mine.back().size(), 0);
+
+      // check whether there is a renumbering of degrees of freedom on
+      // either the finest level or the global dofs, which means that we
+      // cannot apply a plain copy
+      for (unsigned int i=0; i<copy_indices.back().size(); ++i)
+        if (copy_indices.back()[i].first != copy_indices.back()[i].second)
+          {
+            perform_plain_copy = false;
+            break;
+          }
+    }
+  const parallel::Triangulation<dim, spacedim> *ptria =
+    dynamic_cast<const parallel::Triangulation<dim, spacedim> *>
+    (&mg_dof.get_tria());
+  const MPI_Comm mpi_communicator = ptria != 0 ? ptria->get_communicator() :
+                                    MPI_COMM_SELF;
+  perform_plain_copy =
+    Utilities::MPI::min(static_cast<int>(perform_plain_copy),
+                        mpi_communicator);
+
 }