From 4fd8cc6461771f1216837b8f3ffd3b1feceaf992 Mon Sep 17 00:00:00 2001
From: Stefano Zampini <stefano.zampini@gmail.com>
Date: Tue, 18 Apr 2023 11:19:08 +0300
Subject: [PATCH] PETScWrappers::BlockVector add constructors from arrays of
 Vecs

---
 include/deal.II/lac/petsc_block_vector.h | 55 ++++++++++++++++++++----
 tests/petsc/copy_to_dealvec_block.cc     | 15 +++++++
 2 files changed, 61 insertions(+), 9 deletions(-)

diff --git a/include/deal.II/lac/petsc_block_vector.h b/include/deal.II/lac/petsc_block_vector.h
index 17e27e3b87..cc751022a3 100644
--- a/include/deal.II/lac/petsc_block_vector.h
+++ b/include/deal.II/lac/petsc_block_vector.h
@@ -138,6 +138,12 @@ namespace PETScWrappers
        */
       explicit BlockVector(Vec v);
 
+      /**
+       * Create a BlockVector with an array of PETSc vectors.
+       */
+      template <size_t num_blocks>
+      BlockVector(const std::array<Vec, num_blocks> &);
+
       /**
        * Destructor. Clears memory
        */
@@ -361,7 +367,8 @@ namespace PETScWrappers
     /*--------------------- Inline functions --------------------------------*/
 
     inline BlockVector::BlockVector()
-      : petsc_nest_vector(nullptr)
+      : BlockVectorBase<Vector>()
+      , petsc_nest_vector(nullptr)
     {}
 
 
@@ -370,7 +377,7 @@ namespace PETScWrappers
                                     const MPI_Comm &   communicator,
                                     const size_type    block_size,
                                     const size_type    locally_owned_size)
-      : petsc_nest_vector(nullptr)
+      : BlockVector()
     {
       reinit(n_blocks, communicator, block_size, locally_owned_size);
     }
@@ -381,15 +388,15 @@ namespace PETScWrappers
       const std::vector<size_type> &block_sizes,
       const MPI_Comm &              communicator,
       const std::vector<size_type> &local_elements)
-      : petsc_nest_vector(nullptr)
+      : BlockVector()
     {
       reinit(block_sizes, communicator, local_elements, false);
     }
 
 
+
     inline BlockVector::BlockVector(const BlockVector &v)
-      : BlockVectorBase<Vector>()
-      , petsc_nest_vector(nullptr)
+      : BlockVector()
     {
       this->block_indices = v.block_indices;
 
@@ -400,30 +407,51 @@ namespace PETScWrappers
       this->collect_sizes();
     }
 
+
+
     inline BlockVector::BlockVector(
       const std::vector<IndexSet> &parallel_partitioning,
       const MPI_Comm &             communicator)
-      : petsc_nest_vector(nullptr)
+      : BlockVector()
     {
       reinit(parallel_partitioning, communicator);
     }
 
+
+
     inline BlockVector::BlockVector(
       const std::vector<IndexSet> &parallel_partitioning,
       const std::vector<IndexSet> &ghost_indices,
       const MPI_Comm &             communicator)
-      : petsc_nest_vector(nullptr)
+      : BlockVector()
     {
       reinit(parallel_partitioning, ghost_indices, communicator);
     }
 
+
+
     inline BlockVector::BlockVector(Vec v)
-      : BlockVectorBase<Vector>()
-      , petsc_nest_vector(nullptr)
+      : BlockVector()
     {
       this->reinit(v);
     }
 
+
+
+    template <size_t num_blocks>
+    inline BlockVector::BlockVector(const std::array<Vec, num_blocks> &arrayV)
+      : BlockVector()
+    {
+      this->block_indices.reinit(num_blocks, 0);
+
+      this->components.resize(num_blocks);
+      for (auto i = 0; i < num_blocks; ++i)
+        this->components[i].reinit(arrayV[i]);
+      this->collect_sizes();
+    }
+
+
+
     inline BlockVector &
     BlockVector::operator=(const value_type s)
     {
@@ -431,6 +459,8 @@ namespace PETScWrappers
       return *this;
     }
 
+
+
     inline BlockVector &
     BlockVector::operator=(const BlockVector &v)
     {
@@ -500,6 +530,8 @@ namespace PETScWrappers
       this->collect_sizes();
     }
 
+
+
     inline void
     BlockVector::reinit(const std::vector<IndexSet> &parallel_partitioning,
                         const MPI_Comm &             communicator)
@@ -516,6 +548,8 @@ namespace PETScWrappers
       this->collect_sizes();
     }
 
+
+
     inline void
     BlockVector::reinit(const std::vector<IndexSet> &parallel_partitioning,
                         const std::vector<IndexSet> &ghost_entries,
@@ -550,6 +584,8 @@ namespace PETScWrappers
       return comm;
     }
 
+
+
     inline bool
     BlockVector::has_ghost_elements() const
     {
@@ -562,6 +598,7 @@ namespace PETScWrappers
     }
 
 
+
     inline void
     BlockVector::swap(BlockVector &v)
     {
diff --git a/tests/petsc/copy_to_dealvec_block.cc b/tests/petsc/copy_to_dealvec_block.cc
index c6ede99bb4..3654ff4561 100644
--- a/tests/petsc/copy_to_dealvec_block.cc
+++ b/tests/petsc/copy_to_dealvec_block.cc
@@ -148,6 +148,7 @@ test()
   PETScWrappers::MPI::BlockVector vb2(v.petsc_vector());
   Assert(vb2.n_blocks() == v.n_blocks(), ExcInternalError());
   Assert(vb2.size() == v.size(), ExcInternalError());
+  Assert(vb2.petsc_vector() == v.petsc_vector(), ExcInternalError());
   for (unsigned int bl = 0; bl < 2; ++bl)
     {
       Assert(vb2.block(bl).size() == v.block(bl).size(), ExcInternalError());
@@ -155,6 +156,20 @@ test()
              ExcInternalError());
     }
 
+  // Create new block vector from an array of PETSc vectors
+  std::array<Vec, 2> arrayVecs = {
+    {vb.block(0).petsc_vector(), vb.block(1).petsc_vector()}};
+  PETScWrappers::MPI::BlockVector vb3(arrayVecs);
+  Assert(vb3.n_blocks() == vb.n_blocks(), ExcInternalError());
+  Assert(vb3.size() == vb.size(), ExcInternalError());
+  for (unsigned int bl = 0; bl < 2; ++bl)
+    {
+      Assert(vb3.block(bl).size() == vb.block(bl).size(), ExcInternalError());
+      Assert(vb3.block(bl).petsc_vector() == vb.block(bl).petsc_vector(),
+             ExcInternalError());
+    }
+
+
   // Test swap
   auto old_v_vb2 = vb2.petsc_vector();
   auto old_v_vb  = vb.petsc_vector();
-- 
2.39.5