Make the PETSc to deal.II vector copy generic.

We may as well allow copying instances of the base class: the operations are the
same.

diff --git a/include/deal.II/lac/petsc_parallel_vector.h b/include/deal.II/lac/petsc_parallel_vector.h
index 804bb079ca9ec362281d9dfad24a81b70e6fed77..2f07585ff34ad6c47cc1bcfb7f996287e6923f83 100644 (file)
--- a/include/deal.II/lac/petsc_parallel_vector.h
+++ b/include/deal.II/lac/petsc_parallel_vector.h
@@ -36,6 +36,9 @@ DEAL_II_NAMESPACE_OPEN
  */
 namespace PETScWrappers
 {
+  // forward declaration for copy constructor from sequential vector
+  class Vector;
+
   /**
    * Namespace for PETSc classes that work in parallel over MPI, such as
    * distributed vectors and matrices.

diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h
index 99a89bc60133b4bd682f1ae9a76e47ac676750b8..b34b9d5ceff32f31d2cf8711d80dbed809c21b52 100644 (file)
--- a/include/deal.II/lac/vector.h
+++ b/include/deal.II/lac/vector.h
@@ -45,11 +45,7 @@ DEAL_II_NAMESPACE_OPEN
 #ifdef DEAL_II_WITH_PETSC
 namespace PETScWrappers
 {
-  class Vector;
-  namespace MPI
-  {
-    class Vector;
-  }
+  class VectorBase;
 }
 #endif
 
@@ -208,23 +204,17 @@ public:
 
 #ifdef DEAL_II_WITH_PETSC
   /**
-   * Another copy constructor: copy the values from a sequential PETSc wrapper
-   * vector class. This copy constructor is only available if PETSc was
-   * detected during configuration time.
-   */
-  explicit Vector (const PETScWrappers::Vector &v);
-
-  /**
-   * Another copy constructor: copy the values from a parallel PETSc wrapper
-   * vector class. This copy constructor is only available if PETSc was
-   * detected during configuration time.
+   * Another copy constructor: copy the values from a PETSc vector class. This
+   * copy constructor is only available if PETSc was detected during
+   * configuration time.
    *
    * Note that due to the communication model used in MPI, this operation can
-   * only succeed if all processes do it at the same time. I.e., it is not
-   * possible for only one process to obtain a copy of a parallel vector while
-   * the other jobs do something else.
+   * only succeed if all processes do it at the same time when <code>v</code>
+   * is a distributed vector: It is not possible for only one process to
+   * obtain a copy of a parallel vector while the other jobs do something
+   * else.
    */
-  explicit Vector (const PETScWrappers::MPI::Vector &v);
+  explicit Vector (const PETScWrappers::VectorBase &v);
 #endif
 
 #ifdef DEAL_II_WITH_TRILINOS
@@ -382,25 +372,18 @@ public:
 
 #ifdef DEAL_II_WITH_PETSC
   /**
-   * Another copy operator: copy the values from a sequential PETSc wrapper
-   * vector class. This operator is only available if PETSc was detected
-   * during configuration time.
-   */
-  Vector<Number> &
-  operator= (const PETScWrappers::Vector &v);
-
-  /**
-   * Another copy operator: copy the values from a parallel PETSc wrapper
-   * vector class. This operator is only available if PETSc was detected
-   * during configuration time.
+   * Another copy operator: copy the values from a PETSc wrapper vector
+   * class. This operator is only available if PETSc was detected during
+   * configuration time.
    *
    * Note that due to the communication model used in MPI, this operation can
-   * only succeed if all processes do it at the same time. I.e., it is not
-   * possible for only one process to obtain a copy of a parallel vector while
-   * the other jobs do something else.
+   * only succeed if all processes do it at the same time when <code>v</code>
+   * is a distributed vector: It is not possible for only one process to
+   * obtain a copy of a parallel vector while the other jobs do something
+   * else.
    */
   Vector<Number> &
-  operator= (const PETScWrappers::MPI::Vector &v);
+  operator= (const PETScWrappers::VectorBase &v);
 #endif
 
 

diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index 38554edc26d17d6261eb69ce3f67ed56d3873be1..df73ec9ed9a70c1457bbb5ed098d1bdd7721c2ec 100644 (file)
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -25,8 +25,7 @@
 #include <deal.II/lac/vector_operations_internal.h>
 
 #ifdef DEAL_II_WITH_PETSC
-#  include <deal.II/lac/petsc_vector.h>
-#  include <deal.II/lac/petsc_parallel_vector.h>
+#  include <deal.II/lac/petsc_vector_base.h>
 #endif
 
 #ifdef DEAL_II_WITH_TRILINOS
@@ -34,11 +33,12 @@
 #endif
 
 
+#include <algorithm>
 #include <cmath>
 #include <cstring>
-#include <algorithm>
-#include <iostream>
 #include <iomanip>
+#include <iostream>
+#include <memory>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -95,38 +95,56 @@ Vector<Number>::Vector (const Vector<OtherNumber> &v)
 
 
 #ifdef DEAL_II_WITH_PETSC
-
-template <typename Number>
-Vector<Number>::Vector (const PETScWrappers::Vector &v)
-  :
-  Subscriptor(),
-  vec_size(v.size()),
-  max_vec_size(v.size()),
-  val(nullptr)
+namespace internal
 {
-  if (vec_size != 0)
+  template <typename Number>
+  void
+  copy_petsc_vector(const PETScWrappers::VectorBase &v,
+                    ::dealii::Vector<Number> &out)
+  {
+    // Create a sequential PETSc vector and then copy over the entries into
+    // the deal.II vector.
+    //
+    // Wrap the sequential vector with a custom deleter.
+    std::unique_ptr<Vec, std::function<void(Vec *)>> sequential_vector
+                                                  (new Vec(), [](Vec *ptr)
     {
-      allocate();
-
-      // get a representation of the vector
-      // and copy it
-      PetscScalar *start_ptr;
-      PetscErrorCode ierr = VecGetArray (static_cast<const Vec &>(v), &start_ptr);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-      internal::VectorOperations::copy (start_ptr, start_ptr+vec_size, begin());
-
-      // restore the representation of the
-      // vector
-      ierr = VecRestoreArray (static_cast<const Vec &>(v), &start_ptr);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
-    }
+#if DEAL_II_PETSC_VERSION_LT(3,2,0)
+      const PetscErrorCode ierr = VecDestroy (*ptr);
+#else
+      const PetscErrorCode ierr = VecDestroy (ptr);
+#endif
+      (void)ierr;
+      AssertNothrow (ierr == 0, ExcPETScError(ierr));
+    });
+
+    PetscErrorCode ierr = VecCreateSeq (v.get_mpi_communicator(),
+                                        v.size(),
+                                        sequential_vector.get());
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+    ierr = VecCopy (v, *sequential_vector);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    PetscScalar *start_ptr;
+    ierr = VecGetArray(*sequential_vector, &start_ptr);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    const auto v_size = v.size();
+    if (out.size() != v_size)
+      out.reinit (v_size, true);
+
+    internal::VectorOperations::copy (start_ptr,
+                                      start_ptr + out.size(),
+                                      out.begin());
+    ierr = VecRestoreArray (*sequential_vector, &start_ptr);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+  }
 }
 
 
 
 template <typename Number>
-Vector<Number>::Vector (const PETScWrappers::MPI::Vector &v)
+Vector<Number>::Vector (const PETScWrappers::VectorBase &v)
   :
   Subscriptor(),
   vec_size(0),
@@ -135,14 +153,9 @@ Vector<Number>::Vector (const PETScWrappers::MPI::Vector &v)
 {
   if (v.size() != 0)
     {
-      // do this in a two-stage process:
-      // first convert to a sequential petsc
-      // vector, then copy that
-      PETScWrappers::Vector seq (v);
-      *this = seq;
+      internal::copy_petsc_vector(v, *this);
     }
 }
-
 #endif
 
 
@@ -862,47 +875,13 @@ Vector<Number>::operator= (const BlockVector<Number> &v)
 
 
 #ifdef DEAL_II_WITH_PETSC
-
 template <typename Number>
 Vector<Number> &
-Vector<Number>::operator= (const PETScWrappers::Vector &v)
+Vector<Number>::operator= (const PETScWrappers::VectorBase &v)
 {
-  if (v.size() != vec_size)
-    reinit (v.size(), true);
-  if (vec_size != 0)
-    {
-      // get a representation of the vector
-      // and copy it
-      PetscScalar *start_ptr;
-      PetscErrorCode ierr = VecGetArray (static_cast<const Vec &>(v), &start_ptr);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-      internal::VectorOperations::copy (start_ptr, start_ptr+vec_size, begin());
-
-      // restore the representation of the
-      // vector
-      ierr = VecRestoreArray (static_cast<const Vec &>(v), &start_ptr);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
-    }
-
+  internal::copy_petsc_vector(v, *this);
   return *this;
 }
-
-
-
-template <typename Number>
-Vector<Number> &
-Vector<Number>::operator= (const PETScWrappers::MPI::Vector &v)
-{
-  // do this in a two-stage process:
-  // first convert to a sequential petsc
-  // vector, then copy that
-  PETScWrappers::Vector seq (v);
-  *this = seq;
-
-  return *this;
-}
-
 #endif