diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 59f611c8a5bcbfd5e2403a4f75d59a6d9a9ca64b..a32b9d2795ecb15ae0d4173518ae1e2adec8a369 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -96,7 +96,7 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*>  & all)
 {
   gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
   gmshAssembler<double> myAssembler(lsys);
-  gmshElasticityTerm El(0,1.0,.333,getTag());     
+  gmshElasticityTerm El(0, 1.0, .333, getTag());
   
   std::vector<MElement*> v, layer;
 
@@ -436,7 +436,7 @@ void localHarmonicMapping(GModel *gm,
   gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
   gmshAssembler<double> myAssembler(lsys);
   gmshFunction f(1.0);
-  gmshLaplaceTerm Laplace (gm,&f,0);     
+  gmshLaplaceTerm Laplace(gm, &f, 0);
   
   myAssembler.fixVertex ( n1 , 0 , 0 , -1.0);
   myAssembler.fixVertex ( n2 , 0 , 0 , -1.0);
@@ -450,7 +450,7 @@ void localHarmonicMapping(GModel *gm,
 
   gmshLinearSystemGmm<double> *lsys1 = new gmshLinearSystemGmm<double>;
   gmshAssembler<double> myAssembler1(lsys1);
-  gmshLaplaceTerm Laplace1 (gm,&f,1);     
+  gmshLaplaceTerm Laplace1(gm, &f, 1);
   
   myAssembler1.fixVertex ( n2 , 0 , 1 , -1.0);
   myAssembler1.fixVertex ( n3 , 0 , 1 , -1.0);
diff --git a/Numeric/gmshElasticity.h b/Numeric/gmshElasticity.h
index a7dc483e049e1e28aee065b523638f9c8ac1b80c..b0089f2f23195283ee067376aabbc975cc930391 100644
--- a/Numeric/gmshElasticity.h
+++ b/Numeric/gmshElasticity.h
@@ -12,7 +12,7 @@
 #include "MElement.h"
 #include "GmshMatrix.h"
 
-class gmshElasticityTerm : public gmshNodalFemTerm {
+class gmshElasticityTerm : public gmshNodalFemTerm<double> {
   double _E,_nu;
   int _iField;
  protected:
@@ -27,7 +27,7 @@ class gmshElasticityTerm : public gmshNodalFemTerm {
   }
  public:
   gmshElasticityTerm(GModel *gm, double E, double nu, int iField = 1) : 
-    gmshNodalFemTerm(gm), _E(E), _nu(nu), _iField(iField){}
+    gmshNodalFemTerm<double>(gm), _E(E), _nu(nu), _iField(iField){}
   void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
 };
 
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index d1ab8d90e148908ed61450ef05006c71aad94aaa..41188e8dc319edf926bdcc867d74f326ceaecbcc 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -13,7 +13,7 @@
 #include "MElement.h"
 #include "GmshMatrix.h"
 
-class gmshLaplaceTerm : public gmshNodalFemTerm {
+class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
  private:
   const gmshFunction *_diffusivity;
   const int _iField ;
@@ -28,7 +28,7 @@ class gmshLaplaceTerm : public gmshNodalFemTerm {
   }
  public:
   gmshLaplaceTerm(GModel *gm, gmshFunction *diffusivity, int iField = 0) : 
-    gmshNodalFemTerm(gm), _diffusivity(diffusivity), _iField(iField){}
+    gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField){}
   virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
 };
 
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index d7b9fc8dbb802401808819101120801eb8ad9892..2285bbd35a852ec930f1f31b37e0d89cb7f11b62 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -14,18 +14,20 @@
 #include "GModel.h"
 #include "MElement.h"
 
+template<class scalar>
 class gmshTermOfFormulation {  
  protected:
   GModel *_gm;
  public:
   gmshTermOfFormulation(GModel *gm) : _gm(gm) {}
   virtual ~gmshTermOfFormulation(){}
-  virtual void addToMatrix(gmshAssembler<double> &) const = 0;
+  virtual void addToMatrix(gmshAssembler<scalar> &) const = 0;
 };
 
 // a nodal finite element term : variables are always defined at nodes
 // of the mesh
-class gmshNodalFemTerm : public gmshTermOfFormulation {
+template<class scalar>
+class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
  protected:
   // return the number of columns of the element matrix
   virtual int sizeOfC(MElement*) const = 0;
@@ -41,43 +43,44 @@ class gmshNodalFemTerm : public gmshTermOfFormulation {
     getLocalDofR(e, iCol, vC, iCompC, iFieldC);
   }
  public:
-  gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation(gm) {}
+  gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation<scalar>(gm) {}
   virtual ~gmshNodalFemTerm (){}
-  virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const = 0;
-  void addToMatrix(gmshAssembler<double> &lsys) const
+  virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const = 0;
+  void addToMatrix(gmshAssembler<scalar> &lsys) const
   {
-    if (_gm->getNumRegions()){
-      for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
+    GModel *m = gmshTermOfFormulation<scalar>::_gm;
+    if (m->getNumRegions()){
+      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
         addToMatrix(lsys, *it);
       }
     }
-    else if(_gm->getNumFaces()){
-      for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
+    else if(m->getNumFaces()){
+      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
         addToMatrix(lsys, *it);
       }
     }  
   }
-  void addToMatrix(gmshAssembler<double> &lsys, GEntity *ge) const
+  void addToMatrix(gmshAssembler<scalar> &lsys, GEntity *ge) const
   {
     for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
       MElement *e = ge->getMeshElement(i);
       addToMatrix(lsys, e);
     }
   }
-  void addToMatrix(gmshAssembler<double> &lsys, MElement *e) const
+  void addToMatrix(gmshAssembler<scalar> &lsys, MElement *e) const
   {
     const int nbR = sizeOfR(e);
     const int nbC = sizeOfC(e);
-    gmshMatrix<double> localMatrix (nbR, nbC);
+    gmshMatrix<scalar> localMatrix(nbR, nbC);
     elementMatrix(e, localMatrix);
     addToMatrix(lsys, localMatrix, e);
   }
-  void addToMatrix(gmshAssembler<double> &lsys, const std::vector<MElement*> &v) const
+  void addToMatrix(gmshAssembler<scalar> &lsys, const std::vector<MElement*> &v) const
   {
     for (unsigned int i = 0; i < v.size(); i++)
       addToMatrix(lsys, v[i]);
   }
-  void addToMatrix(gmshAssembler<double> &lsys, gmshMatrix<double> &localMatrix, 
+  void addToMatrix(gmshAssembler<scalar> &lsys, gmshMatrix<scalar> &localMatrix, 
                    MElement *e) const
   {
     const int nbR = sizeOfR(e);