Add solveAgainLaplacian() for multiple RHS solve with same matrix

src/maths/Laplacian.hpp

@@ -64,6 +64,7 @@ virtual void addRHSValue(int row, double val) = 0;
 virtual double getRHSValue(const int row) = 0;
 virtual void zeroRHS() = 0;
 virtual double solveLaplacian() = 0;
+virtual double solveAgainLaplacian() = 0;
 virtual std::vector<double> getSolution() = 0;
 virtual std::vector<double> error() = 0;
 virtual void dumpMatrix() = 0;

src/maths/solvers/LaplacianPetsc.cpp

@@ -37,8 +37,6 @@ void LaplacianPetsc::zeroMatrix()
 
 double LaplacianPetsc::solveLaplacian()
 {
-  PetscInt nbIteration;
-  PetscReal norm;
   PetscBool assembled;
 
   if (!MatAssembled(_LaplacianMatrix,&assembled)) {
@@ -46,8 +44,6 @@ double LaplacianPetsc::solveLaplacian()
     MatAssemblyEnd(_LaplacianMatrix, MAT_FINAL_ASSEMBLY);
   }
 
-  VecAssemblyBegin(_RHS);
-  VecAssemblyEnd(_RHS);
   VecDuplicate(_RHS, &_dum);
 
   KSPSetOperators(_ksp, _LaplacianMatrix, _LaplacianMatrix);
@@ -62,6 +58,18 @@ double LaplacianPetsc::solveLaplacian()
       MatSetNullSpace(_LaplacianMatrix, _nullsp);
     }
   }
+
+  return LaplacianPetsc::solveAgainLaplacian();
+}
+
+double LaplacianPetsc::solveAgainLaplacian()
+{
+  PetscInt nbIteration;
+  PetscReal norm;
+
+  VecAssemblyBegin(_RHS);
+  VecAssemblyEnd(_RHS);
+
   KSPSolve(_ksp, _RHS, _dum);
   KSPGetIterationNumber(_ksp, &nbIteration);
   KSPGetResidualNorm(_ksp, &norm);

src/maths/solvers/LaplacianPetsc.hpp

@@ -111,6 +111,13 @@ void dumpRHS();
  */
 double solveLaplacian();
 
+/**
+ * @brief Solve again the LaplacianPetsc (same system, different RHS)
+ *
+ * @brief[out] return residual norm
+ */
+double solveAgainLaplacian();
+
 /**
  * @brief get Solution converted in std::vector
  */