No commit message

No commit message

Solver/dofManager.h

@@ -169,6 +169,17 @@ class dofManager : public dofManagerBase{
     }
     return false;
   }
+  
+  inline bool isAnUnknown(Dof key) const
+  {
+    if(ghostValue.find(key) == ghostValue.end())
+    {
+      if(unknown.find(key) != unknown.end())
+        return true;
+    }
+    return false;
+  }
+
   inline bool isFixed(long int ent, int type) const
   {
     return isFixed(Dof(ent, type));
@@ -210,6 +221,21 @@ class dofManager : public dofManagerBase{
     Vals.resize(originalSize + ndofs);
     for (int i = 0; i < ndofs; ++i) getDofValue(keys[i], Vals[originalSize+i]);
   }
+  
+  inline bool getAnUnknown(Dof key,  dataVec &val) const
+  {
+    if(ghostValue.find(key) == ghostValue.end())
+    {
+      std::map<Dof, int>::const_iterator it = unknown.find(key);
+      if (it != unknown.end())
+      {
+        _current->getFromSolution(it->second, val);
+        return true;
+      }
+    }
+    return false;
+  }
+  
   virtual inline void getDofValue(Dof key,  dataVec &val) const
   {
     {

Solver/linearSystem.h

@@ -20,6 +20,7 @@ class linearSystemBase {
   virtual void clear() = 0;
   virtual void zeroMatrix() = 0;
   virtual void zeroRightHandSide() = 0;
+  virtual void zeroSolution() = 0;
   virtual int systemSolve() = 0;
   void setParameter (std::string key, std::string value);
   virtual void insertInSparsityPattern(int _row, int _col){};
@@ -36,6 +37,7 @@ class linearSystem : public linearSystemBase {
   virtual void addToRightHandSide(int _row, const scalar &val) = 0;
   virtual void getFromRightHandSide(int _row, scalar &val) const = 0;
   virtual void getFromSolution(int _row, scalar &val) const = 0;
+  virtual void addToSolution(int _row, const scalar &val) = 0;
 };
 
 #endif

Solver/linearSystemCSR.h

@@ -125,6 +125,10 @@ class linearSystemCSR : public linearSystem<scalar> {
   {
     if(val != 0.0) (*_b)[row] += val;
   }
+  virtual void addToSolution(int row, const scalar &val)
+  {
+    if(val != 0.0) (*_x)[row] += val;
+  }
   virtual void getFromRightHandSide(int row, scalar &val) const
   {
     val = (*_b)[row];
@@ -145,6 +149,12 @@ class linearSystemCSR : public linearSystem<scalar> {
     if (!_b) return;
     for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
   }
+  virtual void zeroSolution()
+  {
+    if (!_x) return;
+    for(unsigned int i = 0; i < _x->size(); i++) (*_x)[i] = 0.;
+  }
+
   virtual double normInfRightHandSide() const
   {
     double nor = 0.;

Solver/linearSystemFull.h

@@ -54,6 +54,10 @@ class linearSystemFull : public linearSystem<scalar> {
   {
     if(val != 0.0) (*_b)(row) += val;
   }
+  virtual void addToSolution(int row, const scalar &val)
+  {
+    if(val != 0.0) (*_x)(row) += val;
+  }
   virtual void getFromRightHandSide(int row, scalar &val) const
   {
     val = (*_b)(row);
@@ -70,6 +74,10 @@ class linearSystemFull : public linearSystem<scalar> {
   {
     _b->setAll(0.);
   }
+  virtual void zeroSolution()
+  {
+    _x->setAll(0.);
+  }
   virtual double normInfRightHandSide() const{
     double nor = 0.;
     double temp;

Solver/linearSystemGMM.h

@@ -67,6 +67,12 @@ class linearSystemGmm : public linearSystem<scalar> {
   {
     val = (*_b)[row];
   }
+
+  virtual void addToSolution(int row, const scalar &val)
+  {
+    if(val != 0.0) (*_x)[row] += val;
+  }
+
   virtual void getFromSolution(int row, scalar &val) const
   {
     val = (*_x)[row];
@@ -79,6 +85,11 @@ class linearSystemGmm : public linearSystem<scalar> {
   {
     for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
   }
+  virtual void zeroSolution()
+  {
+    for(unsigned int i = 0; i < _x->size(); i++) (*_x)[i] = 0.;
+  }
+
   virtual double normInfRightHandSide() const {
     double nor = 0.;
     double temp;
@@ -120,9 +131,11 @@ class linearSystemGmm : public linearSystem<scalar> {
   virtual void getFromMatrix(int row, int col, scalar &val) const {}
   virtual void addToRightHandSide(int row, const scalar &val) {}
   virtual void getFromRightHandSide(int row, scalar &val) const {}
+  virtual void addToSolution(int row, const scalar &val) {}
   virtual void getFromSolution(int row, scalar &val) const {}
   virtual void zeroMatrix() {}
   virtual void zeroRightHandSide() {}
+  virtual void zeroSolution() {}
   virtual int systemSolve() { return 0; }
   virtual double normInfRightHandSide() const { return 0.; }
   void setPrec(double p){}

Solver/linearSystemPETSc.cpp

@@ -86,6 +86,16 @@ void linearSystemPETScBlockDouble::getFromRightHandSide(int row, fullMatrix<doub
   }
 }
 
+void linearSystemPETScBlockDouble::addToSolution(int row, const fullMatrix<double> &val)
+{
+  for (int ii = 0; ii < _blockSize; ii++) {
+    PetscInt i = row * _blockSize + ii;
+    PetscScalar v = val(ii, 0);
+    VecSetValues(_x, 1, &i, &v, ADD_VALUES);
+  }
+}
+
+
 void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<double> &val) const
 {
   for (int i = 0; i < _blockSize; i++) {
@@ -236,6 +246,16 @@ void linearSystemPETScBlockDouble::zeroRightHandSide()
   }
 }
 
+void linearSystemPETScBlockDouble::zeroSolution()
+{
+  if (_isAllocated) {
+    VecAssemblyBegin(_x);
+    VecAssemblyEnd(_x);
+    VecZeroEntries(_x);
+  }
+}
+
+
 linearSystemPETScBlockDouble::linearSystemPETScBlockDouble()
 {
   _entriesPreAllocated = false;

Solver/linearSystemPETSc.h

@@ -68,9 +68,11 @@ class linearSystemPETSc : public linearSystem<scalar> {
   virtual void getFromRightHandSide(int row, scalar &val) const;
   virtual double normInfRightHandSide() const;
   virtual void addToMatrix(int row, int col, const scalar &val);
+  virtual void addToSolution(int row, const scalar &val);
   virtual void getFromSolution(int row, scalar &val) const;
   virtual void zeroMatrix();
   virtual void zeroRightHandSide();
+  virtual void zeroSolution();
   virtual int systemSolve();
   Mat &getMatrix(){ return _a; }
   std::vector<scalar> getData();
@@ -89,6 +91,7 @@ class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
   void _kspCreate();
   virtual void addToMatrix(int row, int col, const fullMatrix<double> &val);
   virtual void addToRightHandSide(int row, const fullMatrix<double> &val);
+  virtual void addToSolution(int row,  const fullMatrix<double> &val);
   virtual void getFromMatrix(int row, int col, fullMatrix<double> &val ) const;
   virtual void getFromRightHandSide(int row, fullMatrix<double> &val) const;
   virtual void getFromSolution(int row, fullMatrix<double> &val) const;
@@ -99,6 +102,7 @@ class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
   void clear();
   void zeroMatrix();
   void zeroRightHandSide();
+  void zeroSolution();
   double normInfRightHandSide() const;
   void insertInSparsityPattern (int i, int j);
   linearSystemPETScBlockDouble();
@@ -119,10 +123,12 @@ class linearSystemPETSc : public linearSystem<scalar> {
   virtual void addToMatrix(int row, int col, const scalar &val) {}
   virtual void getFromMatrix(int row, int col, scalar &val) const {}
   virtual void addToRightHandSide(int row, const scalar &val) {}
+  virtual void addToSolution(int row, const scalar &val) {}
   virtual void getFromRightHandSide(int row, scalar &val) const {}
   virtual void getFromSolution(int row, scalar &val) const {}
   virtual void zeroMatrix() {}
   virtual void zeroRightHandSide() {}
+  virtual void zeroSolution() {}
   virtual int systemSolve() { return 0; }
   virtual double normInfRightHandSide() const{return 0;}
 };

Solver/linearSystemPETSc.hpp

@@ -213,6 +213,14 @@ void linearSystemPETSc<scalar>::addToMatrix(int row, int col, const scalar &val)
   _try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
 }
 
+template <class scalar>
+void linearSystemPETSc<scalar>::addToSolution(int row, const scalar &val)
+{
+  PetscInt i = row;
+  PetscScalar s = val;
+  _try(VecSetValues(_x, 1, &i, &s, ADD_VALUES));
+}
+
 template <class scalar>
 void linearSystemPETSc<scalar>::getFromSolution(int row, scalar &val) const
 {
@@ -247,6 +255,17 @@ void linearSystemPETSc<scalar>::zeroRightHandSide()
   }
 }
 
+template <class scalar>
+void linearSystemPETSc<scalar>::zeroSolution()
+{
+  if (_isAllocated) {
+    _try(VecAssemblyBegin(_x));
+    _try(VecAssemblyEnd(_x));
+    _try(VecZeroEntries(_x));
+  }
+}
+
+
 template <class scalar>
 int linearSystemPETSc<scalar>::systemSolve()
 {