Routine for summed RHS contribs

src/system/DistributedContext.cpp

@@ -61,6 +61,16 @@ namespace gmshfem::system
       return _distributedExtendedVector;
       #endif
     }
+    PetscErrorCode sumDistributedRHS(Vec localRHSdistributed, Vec* out) const {
+      PetscFunctionBeginUser;
+      Vec res;
+      PetscCall(VecDuplicate(_distributedVector, &res));
+      VecView(localRHSdistributed, PETSC_VIEWER_STDOUT_WORLD);
+      PetscCall(VecScatterBegin(_scatter, localRHSdistributed, res, ADD_VALUES, SCATTER_REVERSE));
+      PetscCall(VecScatterEnd(_scatter, localRHSdistributed, res, ADD_VALUES, SCATTER_REVERSE));
+      *out = res;
+      PetscFunctionReturn(PETSC_SUCCESS);
+    };
     ~DistributedContextImpl() {
       #ifdef HAVE_PETSC
       VecDestroy(&_distributedExtendedVector);
@@ -125,6 +135,16 @@ namespace gmshfem::system
 
     PetscFunctionReturn(PETSC_SUCCESS);
   }
+
+  template< class T_Scalar >
+  PetscErrorCode DistributedContext< T_Scalar >::sumDistributedRHS(Vec localRHSdistributed, Vec *out) const
+  {
+    PetscFunctionBeginUser;
+    if (!_impl)
+      throw common::Exception("Uninitialized DistributedContextImpl");
+    PetscCall(_impl->sumDistributedRHS(localRHSdistributed, out));
+    PetscFunctionReturn(PETSC_SUCCESS);
+  }
 #endif
 
   template <class T_Scalar>

src/system/DistributedContext.h

@@ -45,9 +45,12 @@ class DistributedContext {
     const Indices& localToGlobal() const;
     const Indices& localIDofOwned() const;
     const Indices& localIDofNonOwned() const;
+    // Vec is a distbuted vector. We scatter to an extended vector (duplicated overlaps) then set the local vals
     void readScatteredData(std::vector<T_Scalar>& values, Vec sol) const;
     #ifdef HAVE_PETSC
     PetscErrorCode getSubdomainsIS(IS* is) const;
+    // Map an extended vector into a global vector with sums. Output becomes a standalone vector
+    PetscErrorCode sumDistributedRHS(Vec localRHSdistributed, Vec* out) const;
     #endif
     DistributedContext(Indices&& localToGlobal, Indices&& owned, Indices&& nonOwned);
     ~DistributedContext();

src/system/VectorFactory.h

@@ -59,6 +59,7 @@ namespace gmshfem::system
     Vec getPetsc() const;
     #ifdef HAVE_MPI
     Vec getPetscDistributed(const DistributedContext<T_Scalar>& distributedContext) const;
+    PetscErrorCode getPetscDistributedSummedRHS(const DistributedContext<T_Scalar>& distributedContext, Vec* out) const;
     // Return PetscScalar* but not declared yet
     void* generateRawLocalData(const DistributedContext<T_Scalar>& distributedContext) const;