From e9d19eb1574f38d4716727917d273cd25adf43a5 Mon Sep 17 00:00:00 2001
From: Boris Martin <boris.martin@uliege.be>
Date: Wed, 21 May 2025 14:30:51 +0200
Subject: [PATCH] Routine for summed RHS contribs

---
 src/system/DistributedContext.cpp | 20 ++++++++++++++++++++
 src/system/DistributedContext.h   |  3 +++
 src/system/VectorFactory.h        |  1 +
 3 files changed, 24 insertions(+)

diff --git a/src/system/DistributedContext.cpp b/src/system/DistributedContext.cpp
index 99c37e11..988a8f3a 100644
--- a/src/system/DistributedContext.cpp
+++ b/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>
diff --git a/src/system/DistributedContext.h b/src/system/DistributedContext.h
index b037e163..eff68c99 100644
--- a/src/system/DistributedContext.h
+++ b/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();
diff --git a/src/system/VectorFactory.h b/src/system/VectorFactory.h
index e27d5013..77bc2dbe 100644
--- a/src/system/VectorFactory.h
+++ b/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; 
     
-- 
GitLab