added getGlobalRHS()

src/problem/Formulation.cpp

@@ -1432,6 +1432,41 @@ namespace gmshddm
       return _rhs;
     }
 
+    template< class T_Scalar >
+    gmshfem::algebra::Vector< T_Scalar > Formulation< T_Scalar >::getGlobalRHS()
+    {
+      Vec GlobalRHS;
+      VecScatter ctx;
+      Vec DistributedRHS;
+      _buildPetscRHS(&DistributedRHS);
+
+      // Allocate locally the full vector
+      PetscInt n, N;
+      VecGetLocalSize(DistributedRHS, &n);
+      VecGetSize(DistributedRHS, &N);
+      //VecCreateSeq(PETSC_COMM_SELF, N, &GlobalRHS);
+
+      // Scatter: RHS (distributed) -> GlobalRHS (full copy on each process)
+      VecScatterCreateToAll(DistributedRHS, &ctx, &GlobalRHS);
+      VecScatterBegin(ctx, DistributedRHS, GlobalRHS, INSERT_VALUES, SCATTER_FORWARD);
+      VecScatterEnd(ctx, DistributedRHS, GlobalRHS, INSERT_VALUES, SCATTER_FORWARD);
+
+      // At this point, GlobalRHS on each process contains the entire vector
+      const PetscScalar *petscArray;
+      gmshfem::algebra::Vector<T_Scalar> result(N);
+      PetscInt petscArraySize;
+      VecGetArrayRead(GlobalRHS, &petscArray);
+      for (unsigned i = 0; i < N; ++i) {
+        result[i] = petscArray[i];
+      }
+      VecRestoreArrayRead(GlobalRHS, &petscArray);
+
+      // Clean up
+      VecScatterDestroy(&ctx);
+      VecDestroy(&GlobalRHS);
+      return result;
+    }
+
     template< class T_Scalar >
     unsigned int Formulation< T_Scalar >::numberOfIterations() const
     {
@@ -1637,13 +1672,15 @@ namespace gmshddm
     {
       gmshfem::common::Timer time;
       time.tick();
-      if(mpi::getMPIRank() == 0) {
-        gmshfem::msg::info << "Starting to build a block of RHS." << gmshfem::msg::endl;
-      }
 
       auto locSize = _rhs.size();
       auto nRHS = _activeShots.size();
 
+      if(mpi::getMPIRank() == 0) {
+        gmshfem::msg::info << "Starting to build a block of " << nRHS << "RHS." << gmshfem::msg::endl;
+      }
+
+
       MatCreateDense(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, PETSC_DETERMINE, nRHS, NULL, RHS);
       MatSetUp(*RHS);
 

src/problem/Formulation.h

@@ -138,6 +138,7 @@ namespace gmshddm
 
       void setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution);
       const gmshfem::algebra::Vector< T_Scalar > &getRHS() const;
+      gmshfem::algebra::Vector< T_Scalar > getGlobalRHS(); //Built on the fly, so not a reference. Modifies internally the Petsc representation of the local RHS.
 
 
       unsigned int numberOfIterations() const;