diff --git a/Common/gmshpy.i b/Common/gmshpy.i
index ce22c6191880057819bc1ca51433653487de48ca..be5248da9a00e40c2bcfd94a70aceda1ae913a7d 100644
--- a/Common/gmshpy.i
+++ b/Common/gmshpy.i
@@ -11,7 +11,6 @@
   #include "linearSystem.h"
   #include "linearSystemFull.h"
   #include "linearSystemPETSc.h"
-  #include "linearSystemPETSc.cpp" // needed for the specialization
   #include "linearSystemCSR.h"
   #include "GEntity.h"
   #include "GVertex.h"
@@ -79,15 +78,25 @@ namespace std {
    %template(GFaceList) list<GFace*, std::allocator<GFace*> >;
 }
 
+
+%include "GmshConfig.h"
 %include "Context.h"
 %include "fullMatrix.h"
+%template(fullMatrixDouble) fullMatrix<double>;
+%template(fullVectorDouble) fullVector<double>;
 %include "dofManager.h"
+%template(dofManagerDouble) dofManager<double>;
 %include "GModel.h"
 %include "function.h"
 %include "linearSystem.h"
+%template(linearSystemDouble) linearSystem<double>;
+%template(linearSystemFullMatrixDouble) linearSystem<fullMatrix<double> >;
 %include "linearSystemFull.h"
+%template(linearSystemFullDouble) linearSystemFull<double> ;
 %include "linearSystemPETSc.h"
+%template(linearSystemPETScDouble) linearSystemPETSc<double>;
 %include "linearSystemCSR.h"
+%template(linearSystemTAUCSDouble) linearSystemCSRTaucs<double>;
 %include "GEntity.h"
 %include "GVertex.h"
 %include "GEdge.h"
@@ -104,7 +113,6 @@ namespace std {
 %include "polynomialBasis.h"
 %include "Gauss.h"
 %include "meshPartitionOptions.h"
-%include "linearSystemCSR.h"
 %include "elasticitySolver.h"
 %include "PView.h"
 %include "PViewData.h"
@@ -115,15 +123,5 @@ namespace std {
 %include "SPoint3.h"
 %include "SPoint2.h"
 %include "GPoint.h"  
-
-%template(dofManagerDouble) dofManager<double>;
-%template(linearSystemDouble) linearSystem<double>;
-%template(linearSystemFullMatrixDouble) linearSystem<fullMatrix<double> >;
-%template(linearSystemFullDouble) linearSystemFull<double> ;
-%template(linearSystemPETScDouble) linearSystemPETSc<double>;
-%template(linearSystemTAUCSDouble) linearSystemCSRTaucs<double>;
-%template(fullMatrixDouble) fullMatrix<double>;
-%template(fullVectorDouble) fullVector<double>;
-%template(linearSystemPETScBlockDouble) linearSystemPETSc<fullMatrix<double> >;
 %include "functionPython.h"
 
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 57b212ff7e0f2c090d64fca57bb812a1287762ea..249bef5b8907d75ef55d05f8adb86515cb1be3e9 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -86,7 +86,6 @@ class elasticitySolver
 
   void addDirichletBCLua (int dim, int entityId, int component, std::string luaFunctionName, lua_State *L);
   void addNeumannBCLua (int dim, int entityId, std::string luaFunctionName, lua_State *L);
-  void addNeumannBCFct (int dim, int entityId, const function* luaFunction, lua_State *L);
 
   #endif
 
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 2815bd10c6a6129ce509e8862bb2027ef295c511..5cc15f0b68df6a563e4d8c19c69703d9dfa4517d 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -6,8 +6,15 @@
 #include <stdlib.h>
 #include "GmshMessage.h"
 
-template <>
-void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
+void linearSystemPETScBlockDouble::_kspCreate() {
+  KSPCreate(PETSC_COMM_WORLD, &_ksp);
+  if (this->_parameters.count("petscPrefix"))
+    KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str());
+  KSPSetFromOptions(_ksp);
+  _kspAllocated = true;
+}
+
+void linearSystemPETScBlockDouble::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
 {
   if (!_entriesPreAllocated)
     preAllocateEntries();
@@ -19,7 +26,7 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col,
       modval(jj, ii) = buff;
     }
   PetscInt i = row, j = col;
-  _try(MatSetValuesBlocked(_a, 1, &i, 1, &j, &modval(0,0), ADD_VALUES));
+  MatSetValuesBlocked(_a, 1, &i, 1, &j, &modval(0,0), ADD_VALUES);
   //transpose back so that the original matrix is not modified
   for (int ii = 0; ii < val.size1(); ii++)
     for (int jj = 0; jj < ii; jj++) {
@@ -29,99 +36,175 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col,
     }
 }
 
-template<>
-void linearSystemPETSc<fullMatrix<PetscScalar> >::addToRightHandSide(int row, const fullMatrix<PetscScalar> &val)
+void linearSystemPETScBlockDouble::addToRightHandSide(int row, const fullMatrix<PetscScalar> &val)
 {
   for (int ii = 0; ii < _blockSize; ii++) {
     PetscInt i = row * _blockSize + ii;
     PetscScalar v = val(ii, 0);
-    _try(VecSetValues(_b, 1, &i, &v, ADD_VALUES));
+    VecSetValues(_b, 1, &i, &v, ADD_VALUES);
   }
 }
 
-template<>
-void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromMatrix(int row, int col, fullMatrix<PetscScalar> &val ) const
+void linearSystemPETScBlockDouble::getFromMatrix(int row, int col, fullMatrix<PetscScalar> &val ) const
 {
   Msg::Error("getFromMatrix not implemented for PETSc");
 }
 
-template<>
-void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromRightHandSide(int row, fullMatrix<PetscScalar> &val) const
+void linearSystemPETScBlockDouble::getFromRightHandSide(int row, fullMatrix<PetscScalar> &val) const
 {
-#if defined(PETSC_USE_COMPLEX)
-  PetscScalar *tmp;
-  _try(VecGetArray(_b, &tmp));
-  for (int i = 0; i < _blockSize; i++) {
-    PetscScalar s = tmp[row * _blockSize + i];
-    val(i,0) = s.real();
-  }
-  _try(VecRestoreArray(_b, &tmp));
-#else
   for (int i = 0; i < _blockSize; i++) {
     int ii = row*_blockSize +i;
     VecGetValues ( _b, 1, &ii, &val(i,0));
   }
-#endif
 }
 
-template<>
-void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullMatrix<PetscScalar> &val) const {
-#if defined(PETSC_USE_COMPLEX)
-  PetscScalar *tmp;
-  _try(VecGetArray(_x, &tmp));
-  for (int i = 0; i < _blockSize; i++) {
-    PetscScalar s = tmp[row * _blockSize + i];
-    val(i,0) = s.real();
-  }
-  _try(VecRestoreArray(_x, &tmp));
-#else
+void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<PetscScalar> &val) const 
+{
   for (int i = 0; i < _blockSize; i++) {
     int ii = row*_blockSize +i;
     VecGetValues ( _x, 1, &ii, &val(i,0));
   }
-#endif
 }
 
-template<>
-void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows) 
+void linearSystemPETScBlockDouble::allocate(int nbRows) 
 {
   if (this->_parameters.count("petscOptions"))
-    _try(PetscOptionsInsertString(this->_parameters["petscOptions"].c_str()));
+    PetscOptionsInsertString(this->_parameters["petscOptions"].c_str());
   _blockSize = strtol (_parameters["blockSize"].c_str(), NULL, 10);
   if (_blockSize == 0)
     Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
   clear();
-  //printf("allocate linear system petsc \n");
-  //_try(PetscOptionsInsertString("-ksp_monitor_true_residual -ksp_rtol 1e-10"));
-  _try(MatCreate(PETSC_COMM_WORLD, &_a)); 
-  _try(MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE));
+  MatCreate(PETSC_COMM_WORLD, &_a); 
+  MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE);
   if (Msg::GetCommSize() > 1) {
-    _try(MatSetType(_a, MATMPIBAIJ));
+    MatSetType(_a, MATMPIBAIJ);
   } else {
-    _try(MatSetType(_a, MATSEQBAIJ));
+    MatSetType(_a, MATSEQBAIJ);
   }
   if (_parameters.count("petscPrefix"))
-    _try(MatAppendOptionsPrefix(_a, _parameters["petscPrefix"].c_str()));
-  _try(MatSetFromOptions(_a));
-  _try(MatGetOwnershipRange(_a, &_localRowStart, &_localRowEnd));
-  _try(MatGetSize(_a, &_globalSize, &_localSize));
+    MatAppendOptionsPrefix(_a, _parameters["petscPrefix"].c_str());
+  MatSetFromOptions(_a);
+  MatGetOwnershipRange(_a, &_localRowStart, &_localRowEnd);
+  MatGetSize(_a, &_globalSize, &_localSize);
   _globalSize /= _blockSize;
   _localSize /= _blockSize;
   _localRowStart /= _blockSize;
   _localRowEnd /= _blockSize;
   // override the default options with the ones from the option
   // database (if any)
-  _try(VecCreate(PETSC_COMM_WORLD, &_x));
-  _try(VecSetSizes(_x, nbRows * _blockSize, PETSC_DETERMINE));
+  VecCreate(PETSC_COMM_WORLD, &_x);
+  VecSetSizes(_x, nbRows * _blockSize, PETSC_DETERMINE);
   // override the default options with the ones from the option
   // database (if any)
   if (_parameters.count("petscPrefix"))
-    _try(VecAppendOptionsPrefix(_x, _parameters["petscPrefix"].c_str()));
-  _try(VecSetFromOptions(_x));
-  _try(VecDuplicate(_x, &_b));
+    VecAppendOptionsPrefix(_x, _parameters["petscPrefix"].c_str());
+  VecSetFromOptions(_x);
+  VecDuplicate(_x, &_b);
   _isAllocated = true;
 }
 
+bool linearSystemPETScBlockDouble::isAllocated() const
+{
+  return _isAllocated; 
+}
+
+void linearSystemPETScBlockDouble::clear()
+{
+  if(_isAllocated){
+    MatDestroy(_a);
+    VecDestroy(_x);
+    VecDestroy(_b);
+  }
+  _isAllocated = false;
+}
+
+int linearSystemPETScBlockDouble::systemSolve()
+{
+  if (!_kspAllocated)
+    _kspCreate();
+  if (_parameters["matrix_reuse"] == "same_sparsity")
+    KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN);
+  else if (_parameters["matrix_reuse"] == "same_matrix")
+    KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER);
+  else
+    KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN);
+  MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
+  MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
+  VecAssemblyBegin(_b);
+  VecAssemblyEnd(_b);
+  KSPSolve(_ksp, _b, _x);
+  return 1;
+}
+
+void linearSystemPETScBlockDouble::preAllocateEntries() 
+{
+  if (_entriesPreAllocated) return;
+  if (!_isAllocated) Msg::Fatal("system must be allocated first");
+  if (_sparsity.getNbRows() == 0) {
+    PetscInt prealloc = 300;
+    PetscTruth set;
+    PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
+    if (_blockSize == 0) {
+      MatSeqAIJSetPreallocation(_a, prealloc, PETSC_NULL);
+    } else {
+      MatSeqBAIJSetPreallocation(_a, _blockSize, 5, PETSC_NULL);
+    }
+  } else {
+    std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize);
+    for (int i = 0; i < _localSize; i++) {
+      int n;
+      const int *r = _sparsity.getRow(i, n);
+      for (int j = 0; j < n; j++) {
+        if (r[j] >= _localRowStart && r[j] < _localRowEnd)
+          nByRowDiag[i] ++;
+        else
+          nByRowOffDiag[i] ++;
+      }
+    }
+    if (_blockSize == 0) {
+      MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0]);
+      MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]);
+    } else {
+      MatSeqBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0]);
+      MatMPIBAIJSetPreallocation(_a, _blockSize, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0]);
+    }
+    _sparsity.clear();
+  }
+  _entriesPreAllocated = true;
+}
+
+void linearSystemPETScBlockDouble::zeroMatrix()
+{
+  if (_isAllocated && _entriesPreAllocated) {
+    MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY);
+    MatZeroEntries(_a);
+  }
+}
+
+void linearSystemPETScBlockDouble::zeroRightHandSide()
+{
+  if (_isAllocated) {
+    VecAssemblyBegin(_b);
+    VecAssemblyEnd(_b);
+    VecZeroEntries(_b);
+  }
+}
+
+linearSystemPETScBlockDouble::linearSystemPETScBlockDouble()
+{
+  _entriesPreAllocated = false;
+  _isAllocated = false;
+  _kspAllocated = false;
+  printf("init\n");
+}
+
+double linearSystemPETScBlockDouble::normInfRightHandSide() const
+{
+  PetscReal nor;
+  VecNorm(_b, NORM_INFINITY, &nor);
+  return nor;
+}
 
 #include "Bindings.h"
 void linearSystemPETScRegisterBindings(binding *b) 
@@ -136,12 +219,12 @@ void linearSystemPETScRegisterBindings(binding *b)
   cm->setDescription ("A new PETSc<PetscScalar> solver");
   cb->setParentClass<linearSystem<PetscScalar> >();
   cm->setArgNames(NULL);
-  cb = b->addClass<linearSystemPETSc<fullMatrix<PetscScalar> > >("linearSystemPETScBlock");
+/*  cb = b->addClass<linearSystemPETScBlockDouble >("linearSystemPETScBlock");
   cb->setDescription("A linear system solver, based on PETSc");
-  cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> > >();
-  cm->setDescription ("A new PETScBlock<PetscScalar> solver");
-  cb->setParentClass<linearSystem<fullMatrix<PetscScalar> > >();
-#endif // FIXME
+  cm = cb->setConstructor<linearSystemPETScBlockDouble >();
+  cm->setDescription ("A new PETScBlockDouble solver");
+  cb->setParentClass<linearSystem<fullMatrix<PetscScalar> > >();*/
+  #endif
 }
 
 #endif // HAVE_PETSC
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index 123ba98f67c517937d3809ddd4396db88da45cb9..3afe836e92a9ada09ed5ba0fa6b87a4f5b009b8c 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -236,7 +236,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
       _try(VecZeroEntries(_b));
     }
   }
-  int systemSolve()
+  virtual int systemSolve()
   {
     if (!_kspAllocated)
       _kspCreate();
@@ -266,10 +266,34 @@ class linearSystemPETSc : public linearSystem<scalar> {
 class binding;
 void linearSystemPETScRegisterBindings(binding *b);
 
+class linearSystemPETScBlockDouble : public linearSystem<fullMatrix<double> > {
+  bool _entriesPreAllocated, _isAllocated, _kspAllocated;
+  sparsityPattern _sparsity;
+  Mat _a;
+  Vec _b, _x;
+  KSP _ksp;
+  int _blockSize, _localRowStart, _localRowEnd, _globalSize, _localSize;
+  public:
+  void _kspCreate();
+  virtual void addToMatrix(int row, int col, const fullMatrix<double> &val);
+  virtual void addToRightHandSide(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;
+  void allocate(int nbRows);
+  bool isAllocated() const;
+  int systemSolve();
+  void preAllocateEntries();
+  void clear();
+  void zeroMatrix();
+  void zeroRightHandSide();
+  double normInfRightHandSide() const;
+  linearSystemPETScBlockDouble();
+};
 
 #else
 
-template <class scalar>
+/*template <class scalar>
 class linearSystemPETSc : public linearSystem<scalar> {
  public :
   linearSystemPETSc()
@@ -288,8 +312,10 @@ class linearSystemPETSc : public linearSystem<scalar> {
   virtual void zeroRightHandSide() {}
   virtual int systemSolve() { return 0; }
   virtual double normInfRightHandSide() const{return 0;}
-};
+};*/
+
 
 #endif
 
+
 #endif