From e34ecec095746046bbf8a8ac3bd1517f098c96fb Mon Sep 17 00:00:00 2001
From: Stefano Zampini <stefano.zampini@gmail.com>
Date: Thu, 10 Nov 2022 18:22:19 +0100
Subject: [PATCH] PETScWrappers:BlockSparseMatrix: create PETSc MATNEST

---
 .../deal.II/lac/petsc_block_sparse_matrix.h   | 22 +++++++++++-
 .../lac/petsc_parallel_block_sparse_matrix.cc | 35 +++++++++++++++++++
 tests/petsc/block_matrices_01.cc              |  4 +++
 tests/petsc/block_matrices_01.mpirun=1.output |  9 +++++
 4 files changed, 69 insertions(+), 1 deletion(-)

diff --git a/include/deal.II/lac/petsc_block_sparse_matrix.h b/include/deal.II/lac/petsc_block_sparse_matrix.h
index 913a694b87..8bb025586d 100644
--- a/include/deal.II/lac/petsc_block_sparse_matrix.h
+++ b/include/deal.II/lac/petsc_block_sparse_matrix.h
@@ -105,7 +105,7 @@ namespace PETScWrappers
       /**
        * Destructor.
        */
-      ~BlockSparseMatrix() override = default;
+      ~BlockSparseMatrix() override;
 
       /**
        * Pseudo copy operator only copying empty objects. The sizes of the
@@ -271,6 +271,26 @@ namespace PETScWrappers
        * protected.
        */
       using BlockMatrixBase<SparseMatrix>::clear;
+
+      /**
+       * Conversion operator to gain access to the underlying PETSc type. If you
+       * do this, you cut this class off some information it may need, so this
+       * conversion operator should only be used if you know what you do. In
+       * particular, it should only be used for read-only operations into the
+       * matrix.
+       */
+      operator Mat() const;
+
+      /**
+       * Return a reference to the underlying PETSc type. It can be used to
+       * modify the underlying data, so use it only when you know what you
+       * are doing.
+       */
+      Mat &
+      petsc_matrix();
+
+    private:
+      Mat petsc_nest_matrix = nullptr;
     };
 
 
diff --git a/source/lac/petsc_parallel_block_sparse_matrix.cc b/source/lac/petsc_parallel_block_sparse_matrix.cc
index 6b82eae960..d744d72f2b 100644
--- a/source/lac/petsc_parallel_block_sparse_matrix.cc
+++ b/source/lac/petsc_parallel_block_sparse_matrix.cc
@@ -14,6 +14,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/lac/petsc_block_sparse_matrix.h>
+#include <deal.II/lac/petsc_compatibility.h>
 
 #ifdef DEAL_II_WITH_PETSC
 
@@ -31,6 +32,11 @@ namespace PETScWrappers
       return *this;
     }
 
+    BlockSparseMatrix::~BlockSparseMatrix()
+    {
+      PetscErrorCode ierr = destroy_matrix(petsc_nest_matrix);
+      AssertThrow(ierr == 0, ExcPETScError(ierr));
+    }
 
 #  ifndef DOXYGEN
     void
@@ -111,6 +117,24 @@ namespace PETScWrappers
     BlockSparseMatrix::collect_sizes()
     {
       BaseClass::collect_sizes();
+
+      auto m = this->n_block_cols();
+      auto n = this->n_block_cols();
+
+      PetscErrorCode ierr = destroy_matrix(petsc_nest_matrix);
+      AssertThrow(ierr == 0, ExcPETScError(ierr));
+      std::vector<Mat> psub_objects(m * n);
+      for (unsigned int r = 0; r < m; r++)
+        for (unsigned int c = 0; c < n; c++)
+          psub_objects[r * n + c] = this->sub_objects[r][c]->petsc_matrix();
+      ierr = MatCreateNest(get_mpi_communicator(),
+                           m,
+                           NULL,
+                           n,
+                           NULL,
+                           psub_objects.data(),
+                           &petsc_nest_matrix);
+      AssertThrow(ierr == 0, ExcPETScError(ierr));
     }
 
     std::vector<IndexSet>
@@ -152,6 +176,17 @@ namespace PETScWrappers
       return block(0, 0).get_mpi_communicator();
     }
 
+    BlockSparseMatrix::operator Mat() const
+    {
+      return petsc_nest_matrix;
+    }
+
+    Mat &
+    BlockSparseMatrix::petsc_matrix()
+    {
+      return petsc_nest_matrix;
+    }
+
   } // namespace MPI
 } // namespace PETScWrappers
 
diff --git a/tests/petsc/block_matrices_01.cc b/tests/petsc/block_matrices_01.cc
index 51dc0875dd..193afc3b65 100644
--- a/tests/petsc/block_matrices_01.cc
+++ b/tests/petsc/block_matrices_01.cc
@@ -67,6 +67,10 @@ test()
   pbsm.reinit(rows, cols, bdsp, MPI_COMM_WORLD);
   deallog << "nonzeros BlockSparseMatrix: " << pbsm.n_nonzero_elements()
           << std::endl;
+
+  // Extract the PETSc MATNEST and use print from PETScWrappers::MatrixBase
+  PETScWrappers::MatrixBase tmp(pbsm.petsc_matrix());
+  tmp.print(deallog.get_file_stream());
 }
 
 
diff --git a/tests/petsc/block_matrices_01.mpirun=1.output b/tests/petsc/block_matrices_01.mpirun=1.output
index a9c5e15408..ca11cbdd02 100644
--- a/tests/petsc/block_matrices_01.mpirun=1.output
+++ b/tests/petsc/block_matrices_01.mpirun=1.output
@@ -1,3 +1,12 @@
 
 DEAL::nonzeros BlockSparsityPattern: 9
 DEAL::nonzeros BlockSparseMatrix: 9
+(0,0) 0.00000
+(1,0) 0.00000
+(1,1) 0.00000
+(2,2) 0.00000
+(3,2) 0.00000
+(3,3) 0.00000
+(4,2) 0.00000
+(4,3) 0.00000
+(4,4) 0.00000
-- 
2.39.5