Merge pull request #14735 from luca-heltai/stefanozampini/issue-14656

partially address #14656

cmake/configure/configure_20_petsc.cmake

@@ -44,7 +44,7 @@ macro(feature_petsc_find_external var)
     endif()
 
     #
-    # Petsc has to be configured with the same MPI configuration as
+    # PETSc has to be configured with the same MPI configuration as
     # deal.II.
     #
     # petscconf.h should export PETSC_HAVE_MPIUNI 1 in case mpi support is
@@ -68,7 +68,7 @@ macro(feature_petsc_find_external var)
     endif()
 
     #
-    # Petsc has to be configured with the same number of bits for indices as
+    # PETSc has to be configured with the same number of bits for indices as
     # deal.II.
     #
     # petscconf.h should export PETSC_WITH_64BIT_INDICES 1 in case 64bits

examples/step-44/doc/results.dox

@@ -278,7 +278,7 @@ There are a number of obvious extensions for this work:
 - The program has been developed for solving problems on single-node
   multicore machines. With a little effort, the program could be
   extended to a large-scale computing environment through the use of
-  Petsc or Trilinos, using a similar technique to that demonstrated in
+  PETSc or Trilinos, using a similar technique to that demonstrated in
   step-40. This would mostly involve changes to the setup, assembly,
   <code>PointHistory</code> and linear solver routines.
 - As this program assumes quasi-static equilibrium, extensions to

include/deal.II/lac/petsc_block_sparse_matrix.h

@@ -103,7 +103,11 @@ namespace PETScWrappers
       BlockSparseMatrix() = default;
 
       /**
-       * Create a BlockSparseMatrix with a PETSc Mat
+       * Create a BlockSparseMatrix with a PETSc Mat that describes the entire
+       * block matrix.
+       * It infers the number of blocks from the Mat if it is of type MATNEST,
+       * otherwise the block operator will only have a single block.
+       * Internally, we always store a MATNEST matrix.
        */
       explicit BlockSparseMatrix(const Mat &);
 
@@ -173,6 +177,14 @@ namespace PETScWrappers
              const MPI_Comm &                   com);
 
 
+      /**
+       * This method associates the PETSc Mat to the instance of the class.
+       * Infers the number of blocks from A if it is of type MATNEST, otherwise
+       * the block operator will only have a single block.
+       */
+      void
+      reinit(Mat A);
+
 
       /**
        * Matrix-vector multiplication: let $dst = M*src$ with $M$ being this
@@ -292,29 +304,18 @@ namespace PETScWrappers
        * are doing.
        *
        * The PETSc Mat object returned here describes the *entire* matrix,
-       * not just one of its blocks. Internally, this is done by creating
-       * a "nested" matrix using Petsc's MatCreateNest object whose individual
+       * not just one of its blocks. Internally, this is done using
+       * a "nested" matrix using PETSc's MATNEST object whose individual
        * blocks are the blocks of this matrix.
        */
       Mat &
       petsc_matrix();
 
-      /**
-       * This method assigns the PETSc Mat to the instance of the class.
-       *
-       * Note that the matrix is not copied: instead, the instance of this class
-       * is initialized to use the given matrix. This is useful if you want to
-       * interpret a PETSc Mat object as a deal.II BlockMatrix, and you already
-       * have a BlockMatrix object that you want to use for this purpose.
-       */
-      void
-      assign_petsc_matrix(Mat A);
-
     private:
       /**
        * A PETSc Mat object that describes the entire block matrix.
        * Internally, this is done by creating
-       * a "nested" matrix using Petsc's MatCreateNest object whose individual
+       * a "nested" matrix using PETSc's MATNEST object whose individual
        * blocks are the blocks of this matrix.
        */
       Mat petsc_nest_matrix = nullptr;
@@ -329,7 +330,7 @@ namespace PETScWrappers
     inline BlockSparseMatrix::BlockSparseMatrix(const Mat &A)
       : BlockSparseMatrix()
     {
-      this->assign_petsc_matrix(A);
+      this->reinit(A);
     }
 
     inline BlockSparseMatrix &

include/deal.II/lac/petsc_block_vector.h

@@ -132,10 +132,12 @@ namespace PETScWrappers
 
       /**
        * Create a BlockVector with a PETSc Vec
+       * It infers the number of blocks from the Vec if it is of type VECNEST,
+       * otherwise the block vector will only have a single block.
+       * Internally, we always store a VECNEST vector.
        */
       explicit BlockVector(Vec v);
 
-
       /**
        * Destructor. Clears memory
        */
@@ -155,15 +157,12 @@ namespace PETScWrappers
       operator=(const BlockVector &V);
 
       /**
-       * This method assigns the given PETSc Vec to the instance of the class.
-       *
-       * Note that the vector is not copied: instead, the instance of this class
-       * is initialized to use the given vector. This is useful if you want to
-       * interpret a PETSc vector as a deal.II vector, and you already have a
-       * BlockVector that you want to use for this purpose.
+       * This method associates the PETSc Vec to the instance of the class.
+       * Infers the number of blocks from v if it is of type VECNEST, otherwise
+       * the block vector will only have a single block.
        */
       void
-      assign_petsc_vector(Vec v);
+      reinit(Vec v);
 
       /**
        * Reinitialize the BlockVector to contain @p n_blocks of size @p
@@ -399,7 +398,7 @@ namespace PETScWrappers
       : BlockVectorBase<Vector>()
       , petsc_nest_vector(nullptr)
     {
-      this->assign_petsc_vector(v);
+      this->reinit(v);
     }
 
     inline BlockVector &

include/deal.II/lac/petsc_matrix_base.h

@@ -345,11 +345,13 @@ namespace PETScWrappers
     virtual ~MatrixBase() override;
 
     /**
-     * This method assigns the PETSc Mat to the instance of the class.
-     *
+     * This method associates the PETSc Mat to the instance of the class.
+     * This is particularly useful when performing PETSc to Deal.II operations
+     * since it allows to reuse the Deal.II MatrixBase and the PETSc Mat
+     * without incurring in memory copies.
      */
     void
-    assign_petsc_matrix(Mat A);
+    reinit(Mat A);
 
     /**
      * This operator assigns a scalar to a matrix. Since this does usually not

include/deal.II/lac/petsc_sparse_matrix.h

@@ -77,7 +77,9 @@ namespace PETScWrappers
 
     /**
      * Initialize a Matrix from a PETSc Mat object. Note that we do not copy
-     * the matrix.
+     * the matrix. The Mat object is referenced by the newly created instance
+     * of the class using PetscObjectReference. This is in line with the PETSc
+     * approach to object ownership, which mimicks std::shared_ptr.
      */
     explicit SparseMatrix(const Mat &);
 
@@ -396,7 +398,9 @@ namespace PETScWrappers
 
       /**
        * Initialize a SparseMatrix from a PETSc Mat object. Note that we do not
-       * copy the matrix.
+       * copy the matrix. The Mat object is referenced by the newly created
+       * instance of the class using PetscObjectReference. This is in line with
+       * the PETSc approach to object ownership, which mimicks std::shared_ptr.
        */
       explicit SparseMatrix(const Mat &);
 

include/deal.II/lac/petsc_vector.h

@@ -289,6 +289,8 @@ namespace PETScWrappers
       Vector &
       operator=(const dealii::Vector<number> &v);
 
+      using VectorBase::reinit;
+
       /**
        * Change the dimension of the vector to @p N. It is unspecified how
        * resizing the vector affects the memory allocation of this object;
@@ -476,7 +478,7 @@ namespace PETScWrappers
       // in this case a) again takes up a whole lot of memory on the heap,
       // and b) is totally dumb since its content would simply be the
       // sequence 0,1,2,3,...,n. the best of all worlds would probably be a
-      // function in Petsc that would take a pointer to an array of
+      // function in PETSc that would take a pointer to an array of
       // PetscScalar values and simply copy n elements verbatim into the
       // vector...
       for (size_type i = 0; i < v.size(); ++i)

include/deal.II/lac/petsc_vector_base.h

@@ -327,13 +327,13 @@ namespace PETScWrappers
     operator=(const PetscScalar s);
 
     /**
-     * This method assigns the PETSc Vec to the instance of the class.
-     * This is particularly useful when performing PETSc to deal.II operations
-     * since it allows one to reuse the VectorBase and the PETSc Vec
+     * This method associates the PETSc Vec to the instance of the class.
+     * This is particularly useful when performing PETSc to Deal.II operations
+     * since it allows to reuse the Deal.II VectorBase and the PETSc Vec
      * without incurring in memory copies.
      */
     void
-    assign_petsc_vector(Vec v);
+    reinit(Vec v);
 
     /**
      * Test for equality. This function assumes that the present vector and

include/deal.II/matrix_free/vector_access_internal.h

@@ -37,7 +37,7 @@ namespace internal
   // below we use type-traits from matrix-free/type_traits.h
 
   // access to generic const vectors that have operator ().
-  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
+  // FIXME: this is wrong for Trilinos/PETSc MPI vectors
   // where we should first do Partitioner::local_to_global()
   template <
     typename VectorType,
@@ -51,7 +51,7 @@ namespace internal
 
 
   // access to generic non-const vectors that have operator ().
-  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
+  // FIXME: this is wrong for Trilinos/PETSc MPI vectors
   // where we should first do Partitioner::local_to_global()
   template <
     typename VectorType,

source/lac/petsc_matrix_base.cc

@@ -90,7 +90,7 @@ namespace PETScWrappers
   }
 
   void
-  MatrixBase::assign_petsc_matrix(Mat A)
+  MatrixBase::reinit(Mat A)
   {
     AssertThrow(last_action == ::dealii::VectorOperation::unknown,
                 ExcMessage("Cannot assign a new Mat."));

source/lac/petsc_parallel_block_sparse_matrix.cc

@@ -210,7 +210,7 @@ namespace PETScWrappers
     }
 
     void
-    BlockSparseMatrix::assign_petsc_matrix(Mat A)
+    BlockSparseMatrix::reinit(Mat A)
     {
       clear();
 

source/lac/petsc_parallel_block_vector.cc

@@ -46,7 +46,7 @@ namespace PETScWrappers
     }
 
     void
-    BlockVector::assign_petsc_vector(Vec v)
+    BlockVector::reinit(Vec v)
     {
       PetscBool isnest;
 
@@ -81,7 +81,7 @@ namespace PETScWrappers
       this->components.resize(nb);
       for (unsigned int i = 0; i < nb; ++i)
         {
-          this->components[i].assign_petsc_vector(sv[i]);
+          this->components[i].reinit(sv[i]);
         }
 
       this->collect_sizes();

source/lac/petsc_vector_base.cc

@@ -181,7 +181,7 @@ namespace PETScWrappers
 
 
   void
-  VectorBase::assign_petsc_vector(Vec v)
+  VectorBase::reinit(Vec v)
   {
     /* TODO GHOSTED */
     AssertThrow(last_action == ::dealii::VectorOperation::unknown,