diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 0a95a318fc901a62163be271b5d0b1382476d5a7..34d6d3fa04aa28b5dd1c79517790d970fe5414b4 100644
--- a/Solver/dofManager.h
+++ b/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
   {
     {
diff --git a/Solver/linearSystem.h b/Solver/linearSystem.h
index 6a500299943883fa5b6d1f655f7e95a7f7c5a8f2..89285413f338ef94a0e9c212b3cac462253ac543 100644
--- a/Solver/linearSystem.h
+++ b/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
diff --git a/Solver/linearSystemCSR.h b/Solver/linearSystemCSR.h
index 17ae925813ab9a9e278e41d4926150eef5432cfd..ac51f15ba674b6540dc9bca0074a80f6da414289 100644
--- a/Solver/linearSystemCSR.h
+++ b/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.;
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index 2703f41acfec06340cd7e84c37dbccc0b9a2180a..4e22516cf527600c90f2a10a979ae01f5d00457f 100644
--- a/Solver/linearSystemFull.h
+++ b/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;
diff --git a/Solver/linearSystemGMM.h b/Solver/linearSystemGMM.h
index 0c0daa0810427be68d835a8cd9f26dac4c56fa01..96072966fb607e539de21571a283eeff21d4e61a 100644
--- a/Solver/linearSystemGMM.h
+++ b/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){}
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 2ea311cdf982a7f3987a3c2be52451373d95abd7..725abf46e3aac12b30726e507a22588bbc045267 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/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;
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 045c7ec2bdf617d7427c475cb2340bfcb10044a3..24dcfe64df2ad699d656a501ba97a8b59775ed29 100644
--- a/Solver/linearSystemPETSc.h
+++ b/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;}
 };
diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp
index 4cfd224d9665dea5c524dc8fa0180c3499e7932a..e017cf133a9489d5cc431e575e37fc12fa8eb611 100644
--- a/Solver/linearSystemPETSc.hpp
+++ b/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()
 {