diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 973133fa9514b3bed64d0757c1c569012f59b8a3..6afed657898f1ae4a607d076a06f3ef7ac442f77 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -16,7 +16,7 @@
 #include "Octree.h"
 #include "SBoundingBox3d.h"
 #include "SPoint3.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
 #include "robustPredicates.h"
 #include "MElementCut.h"
 #include "GEntity.h"
@@ -34,7 +34,7 @@
 #include "discreteFace.h"
 #include "eigenSolver.h"
 
-static void fixEdgeToValue(GEdge *ed, double value, dofManager<double, double> &myAssembler)
+static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
 {
   myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value);
   myAssembler.fixVertex(ed->getEndVertex()->mesh_vertices[0], 0, 1, value);
@@ -869,7 +869,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
 void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
 {
    
-  dofManager<double, double> myAssembler(_lsys);
+  dofManager<double> myAssembler(_lsys);
   simpleFunction<double> ONE(1.0);
 
   if(_type == UNITCIRCLE){
@@ -976,7 +976,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
 void GFaceCompound::parametrize_conformal() const
 {
 
-  dofManager<double, double> myAssembler(_lsys);
+  dofManager<double> myAssembler(_lsys);
 
   std::vector<MVertex*> ordered;
   std::vector<double> coords;  
@@ -1059,7 +1059,7 @@ void GFaceCompound::compute_distance() const
   double L = norm(SVector3(bbox.max(), bbox.min())); 
   double mu = L/28;
   simpleFunction<double> DIFF(mu * mu), MONE(1.0);
-  dofManager<double, double> myAssembler(_lsys);
+  dofManager<double> myAssembler(_lsys);
   distanceTerm distance(model(), 1, &DIFF, &MONE);
 
   std::vector<MVertex*> ordered;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 01b7398403aff03306f94ec4c8159d10c3b689fe..ca3ba8a384cd62f1b9ef13008723e48ae2085183 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -110,7 +110,7 @@ double MElement::rhoShapeMeasure()
 
 void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
 {
-  const functionSpace* fs = getFunctionSpace(o);
+  const polynomialBasis* fs = getFunctionSpace(o);
   if(fs) fs->f(u, v, w, s);
   else Msg::Error("Function space not implemented for this type of element");
 }
@@ -118,7 +118,7 @@ void MElement::getShapeFunctions(double u, double v, double w, double s[], int o
 void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3], 
                                      int o)
 {
-  const functionSpace* fs = getFunctionSpace(o);
+  const polynomialBasis* fs = getFunctionSpace(o);
   if(fs) fs->df(u, v, w, s);
   else Msg::Error("Function space not implemented for this type of element");
 }
diff --git a/Geo/MElement.h b/Geo/MElement.h
index a35759b733606fc3b04d390da54fe0b12179a30b..44f6c2b73c3504d5ec3de4a8a437303045fa4049 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -14,7 +14,7 @@
 #include "MVertex.h"
 #include "MEdge.h"
 #include "MFace.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
 #include "Gauss.h"
 
 class GFace;
@@ -173,7 +173,7 @@ class MElement
   virtual std::string getInfoString();
 
   // get the function space for the element
-  virtual const functionSpace* getFunctionSpace(int order=-1) const { return 0; }
+  virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
   
   // return the interpolating nodal shape functions evaluated at point
   // (u,v,w) in parametric coordinates (if order == -1, use the
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index fa37382c5a5d12df1e5202e392f494840078f6e6..a9f355c15d4665256e53526cbd90b13e2430b16f 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -105,7 +105,7 @@ class MPolyhedron : public MElement {
     return vol;
   }
   virtual int getVolumeSign() { return (getVolume() >= 0) ? 1 : -1; }
-  virtual const functionSpace* getFunctionSpace(int order=-1) const 
+  virtual const polynomialBasis* getFunctionSpace(int order=-1) const 
   {
     return _orig->getFunctionSpace(order);
   }
@@ -207,7 +207,7 @@ class MPolygon : public MElement {
   }
   virtual int getType() const { return TYPE_POLYG; }
   virtual int getTypeForMSH() const { return MSH_POLYG_; }
-  virtual const functionSpace* getFunctionSpace(int order=-1) const 
+  virtual const polynomialBasis* getFunctionSpace(int order=-1) const 
   {
     return _orig->getFunctionSpace(order);
   }
@@ -251,7 +251,7 @@ class MTriangleBorder : public MTriangle {
     if(_domains[1]) return _domains[1]->getParent();
     return NULL;
   }
-  virtual const functionSpace* getFunctionSpace(int order=-1) const 
+  virtual const polynomialBasis* getFunctionSpace(int order=-1) const 
   {
     return getParent()->getFunctionSpace(order);
   }
@@ -287,7 +287,7 @@ class MLineBorder : public MLine {
     if(_domains[1]) return _domains[1]->getParent();
     return NULL;
   }
-  virtual const functionSpace* getFunctionSpace(int order=-1) const 
+  virtual const polynomialBasis* getFunctionSpace(int order=-1) const 
   {
     return getParent()->getFunctionSpace(order);
   }
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index cdd91eb117e667e31bc0270a87901eb1c59271bc..a4f906ffa2d6625a70db2832551844d2ca96053e 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -8,16 +8,16 @@
 #include "Context.h"
 #include "qualityMeasures.h"
 
-const functionSpace* MLine::getFunctionSpace(int o) const
+const polynomialBasis* MLine::getFunctionSpace(int o) const
 {
   int order = (o == -1) ? getPolynomialOrder() : o;
   
   switch (order) {
-  case 1: return &functionSpaces::find(MSH_LIN_2);
-  case 2: return &functionSpaces::find(MSH_LIN_3);
-  case 3: return &functionSpaces::find(MSH_LIN_4);
-  case 4: return &functionSpaces::find(MSH_LIN_5);
-  case 5: return &functionSpaces::find(MSH_LIN_6);
+  case 1: return &polynomialBases::find(MSH_LIN_2);
+  case 2: return &polynomialBases::find(MSH_LIN_3);
+  case 3: return &polynomialBases::find(MSH_LIN_4);
+  case 4: return &polynomialBases::find(MSH_LIN_5);
+  case 5: return &polynomialBases::find(MSH_LIN_6);
   default: Msg::Error("Order %d line function space not implemented", order);
   }
   return 0;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index c78fe34652fa691f3a1ab38a0f7209d8f3d3a66e..74542824aadfa3c4a0c9cf799eeb9e8dfc59fd50 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -63,7 +63,7 @@ class MLine : public MElement {
   {
     MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
   }
-  virtual const functionSpace* getFunctionSpace(int o=-1) const;
+  virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 3daf1324c1f5ad2b7a5b98cf18d0fae633725fed..aec951b0c4d3bd912b6418fbf9256b1ec268eb67 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -8,7 +8,7 @@
 #include "Context.h"
 #include "qualityMeasures.h"
 #define SQU(a)      ((a)*(a))
-const functionSpace* MQuadrangle::getFunctionSpace(int o) const
+const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
 {
   int order = (o == -1) ? getPolynomialOrder() : o;
 
@@ -16,19 +16,19 @@ const functionSpace* MQuadrangle::getFunctionSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-      case 1: return &functionSpaces::find(MSH_QUA_4);
-      case 2: return &functionSpaces::find(MSH_QUA_8);
-      case 3: return &functionSpaces::find(MSH_QUA_12);
-      case 4: return &functionSpaces::find(MSH_QUA_16I);
-      case 5: return &functionSpaces::find(MSH_QUA_20);
+      case 1: return &polynomialBases::find(MSH_QUA_4);
+      case 2: return &polynomialBases::find(MSH_QUA_8);
+      case 3: return &polynomialBases::find(MSH_QUA_12);
+      case 4: return &polynomialBases::find(MSH_QUA_16I);
+      case 5: return &polynomialBases::find(MSH_QUA_20);
     }
   }
   switch (order) {
-    case 1: return &functionSpaces::find(MSH_QUA_4);
-    case 2: return &functionSpaces::find(MSH_QUA_9);
-    case 3: return &functionSpaces::find(MSH_QUA_16);
-    case 4: return &functionSpaces::find(MSH_QUA_25);
-    case 5: return &functionSpaces::find(MSH_QUA_36);
+    case 1: return &polynomialBases::find(MSH_QUA_4);
+    case 2: return &polynomialBases::find(MSH_QUA_9);
+    case 3: return &polynomialBases::find(MSH_QUA_16);
+    case 4: return &polynomialBases::find(MSH_QUA_25);
+    case 5: return &polynomialBases::find(MSH_QUA_36);
     default: Msg::Error("Order %d triangle function space not implemented", order);
   }
   return 0;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index d6457574b0787d6b698561b9d928599ff3970f65..e4e5cbb3734f26c03384b878084219c01a5e3a18 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -115,7 +115,7 @@ class MQuadrangle : public MElement {
   virtual const char *getStringForPOS() const { return "SQ"; }
   virtual const char *getStringForBDF() const { return "CQUAD4"; }
   virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
-  virtual const functionSpace* getFunctionSpace(int o=-1) const;
+  virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual void revert() 
   {
     MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 9875212bffbcac1617a5b0b72f34400286df1a09..e79da9bb508f2df41efaee2289454b87ae6fb96d 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -83,7 +83,7 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3])
   sys3x3(mat, b, uvw, &det);
 }
 
-const functionSpace* MTetrahedron::getFunctionSpace(int o) const
+const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
 {
   int order = (o == -1) ? getPolynomialOrder() : o;
 
@@ -91,21 +91,21 @@ const functionSpace* MTetrahedron::getFunctionSpace(int o) const
   
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 1: return &functionSpaces::find(MSH_TET_4);
-    case 2: return &functionSpaces::find(MSH_TET_10);
-    case 3: return &functionSpaces::find(MSH_TET_20);
-    case 4: return &functionSpaces::find(MSH_TET_34);
-    case 5: return &functionSpaces::find(MSH_TET_52);
+    case 1: return &polynomialBases::find(MSH_TET_4);
+    case 2: return &polynomialBases::find(MSH_TET_10);
+    case 3: return &polynomialBases::find(MSH_TET_20);
+    case 4: return &polynomialBases::find(MSH_TET_34);
+    case 5: return &polynomialBases::find(MSH_TET_52);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
     }
   }
   else { 
     switch (order) {
-    case 1: return &functionSpaces::find(MSH_TET_4);
-    case 2: return &functionSpaces::find(MSH_TET_10);
-    case 3: return &functionSpaces::find(MSH_TET_20);
-    case 4: return &functionSpaces::find(MSH_TET_35);
-    case 5: return &functionSpaces::find(MSH_TET_56);
+    case 1: return &polynomialBases::find(MSH_TET_4);
+    case 2: return &polynomialBases::find(MSH_TET_10);
+    case 3: return &polynomialBases::find(MSH_TET_20);
+    case 4: return &polynomialBases::find(MSH_TET_35);
+    case 5: return &polynomialBases::find(MSH_TET_56);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
     }
   }
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index dec9e6e20a47ae64b3b8118880fc5c8aa2ed4ac3..bc98392d9aa00aefe1dfa6af257514863d338a22 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -130,7 +130,7 @@ class MTetrahedron : public MElement {
   virtual double distoShapeMeasure();
   virtual double etaShapeMeasure();
   void xyz2uvw(double xyz[3], double uvw[3]);
-  virtual const functionSpace* getFunctionSpace(int o=-1) const;
+  virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index cd97bce754fca60bf46622861a0fa3c603d3311a..45b1b501c7e1fbea9de5a532a033a3612c6b1e04 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -41,7 +41,7 @@ double MTriangle::gammaShapeMeasure()
 #endif
 }
 
-const functionSpace* MTriangle::getFunctionSpace(int o) const
+const polynomialBasis* MTriangle::getFunctionSpace(int o) const
 {
   int order = (o == -1) ? getPolynomialOrder() : o;
 
@@ -49,21 +49,21 @@ const functionSpace* MTriangle::getFunctionSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-    case 1: return &functionSpaces::find(MSH_TRI_3);
-    case 2: return &functionSpaces::find(MSH_TRI_6);
-    case 3: return &functionSpaces::find(MSH_TRI_9);
-    case 4: return &functionSpaces::find(MSH_TRI_12);
-    case 5: return &functionSpaces::find(MSH_TRI_15I);
+    case 1: return &polynomialBases::find(MSH_TRI_3);
+    case 2: return &polynomialBases::find(MSH_TRI_6);
+    case 3: return &polynomialBases::find(MSH_TRI_9);
+    case 4: return &polynomialBases::find(MSH_TRI_12);
+    case 5: return &polynomialBases::find(MSH_TRI_15I);
     default: Msg::Error("Order %d triangle function space not implemented", order);
     }
   }
   else { 
     switch (order) {
-    case 1: return &functionSpaces::find(MSH_TRI_3);
-    case 2: return &functionSpaces::find(MSH_TRI_6);
-    case 3: return &functionSpaces::find(MSH_TRI_10);
-    case 4: return &functionSpaces::find(MSH_TRI_15);
-    case 5: return &functionSpaces::find(MSH_TRI_21);
+    case 1: return &polynomialBases::find(MSH_TRI_3);
+    case 2: return &polynomialBases::find(MSH_TRI_6);
+    case 3: return &polynomialBases::find(MSH_TRI_10);
+    case 4: return &polynomialBases::find(MSH_TRI_15);
+    case 5: return &polynomialBases::find(MSH_TRI_21);
     default: Msg::Error("Order %d triangle function space not implemented", order);
     }
   }
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index e59c1b159a91510807db523518289e8e8e5d6bb7..67172ff3581b9e3ca7b213f6d151579eb3b0554f 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -121,7 +121,7 @@ class MTriangle : public MElement {
   {
     MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
   }
-  virtual const functionSpace* getFunctionSpace(int o=-1) const;
+  virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 6c361c397fe026cf67996a717bfd2193813143d6..eecf866a62143a2ae9963675e78f198528dc6155 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -21,7 +21,7 @@
 #include "Numeric.h"
 #include "Context.h"
 #include "fullMatrix.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
 
 #define SQU(a)      ((a)*(a))
 
@@ -521,15 +521,15 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
   
   switch (nPts){
   case 2 :
-    points = functionSpaces::find(MSH_TRI_10).points;
+    points = polynomialBases::find(MSH_TRI_10).points;
     start = 9;
     break;
   case 3 :
-    points = functionSpaces::find(MSH_TRI_15).points;
+    points = polynomialBases::find(MSH_TRI_15).points;
     start = 12;
     break;
   case 4 :
-    points = functionSpaces::find(MSH_TRI_21).points;
+    points = polynomialBases::find(MSH_TRI_21).points;
     start = 15;
     break;
   default :  
@@ -614,11 +614,11 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
 
   switch (nPts){
   case 3 :
-    points = functionSpaces::find(MSH_TET_35).points;
+    points = polynomialBases::find(MSH_TET_35).points;
     start = 34;
     break;
   case 4 :
-    points = functionSpaces::find(MSH_TET_56).points;
+    points = polynomialBases::find(MSH_TET_56).points;
     start = 52;
     break;
   default :  
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index dc56882ddd2a68ab793d8150597e613422ee34ac..e6e5168e5325ad08b6436ab5f67ef4bccf5e996e 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -434,7 +434,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
   lsys->setGmres(1);
   lsys->setPrec(5.e-8);
 #endif
-  dofManager<double,double> myAssembler(lsys);
+  dofManager<double> myAssembler(lsys);
   elasticityTerm El(0, 1.0, .333, getTag());
   
   std::vector<MElement*> layer, v;
@@ -540,7 +540,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*>  & all, GFace *gf)
 
 double highOrderSmoother::smooth_metric_(std::vector<MElement*>  & v, 
                                              GFace *gf, 
-                                             dofManager<double,double> &myAssembler,
+                                             dofManager<double> &myAssembler,
                                              std::set<MVertex*> &verticesToMove,
                                              elasticityTerm &El)
 {
@@ -609,7 +609,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
   lsys->setGmres(1);
   lsys->setPrec(5.e-8);
 #endif
-  dofManager<double,double> myAssembler(lsys);
+  dofManager<double> myAssembler(lsys);
   elasticityTerm El(0, 1.0, .333, getTag());
   
   std::vector<MElement*> layer, v;
diff --git a/Mesh/highOrderSmoother.h b/Mesh/highOrderSmoother.h
index 8b16f060d211544d9d68d4d977a1146424a54b32..3b4379fbc00a6e4678dc1174fb8c6e2c42c3b046 100644
--- a/Mesh/highOrderSmoother.h
+++ b/Mesh/highOrderSmoother.h
@@ -44,7 +44,7 @@ public:
   }  
   void smooth(std::vector<MElement*> & );
   double smooth_metric_(std::vector<MElement*> &, GFace *gf,
-                        dofManager<double,double> &myAssembler,
+                        dofManager<double> &myAssembler,
                         std::set<MVertex*> &verticesToMove,
                         elasticityTerm &El);
   void smooth_metric(std::vector<MElement*> &, GFace *gf );
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index c35eb7284805c1bc7deb48cf53fd676023bdcd32..a4bc0b5b97e4826fd920ede441b377ce7723d548 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -9,7 +9,7 @@
 #include "MTriangle.h"
 #include "MTetrahedron.h"
 #include "Numeric.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
 #include "GmshMessage.h"
 
 double qmTriangle(const BDS_Point *p1, const BDS_Point *p2, const BDS_Point *p3, 
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index bd8066edd235bc240406c0019ccbc42f6ca5a3d4..e2e2451b4daa8854655f3f1acb07b7ac059b9ffc 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -6,7 +6,7 @@
 set(SRC
   Numeric.cpp
     fullMatrix.cpp
-    functionSpace.cpp
+    polynomialBasis.cpp
   GaussQuadratureLin.cpp
     GaussQuadratureTri.cpp
     GaussQuadratureQuad.cpp
diff --git a/Numeric/functionSpace.cpp b/Numeric/polynomialBasis.cpp
similarity index 94%
rename from Numeric/functionSpace.cpp
rename to Numeric/polynomialBasis.cpp
index fa2f604eb4fc22a3285ff04f872e917b49f3b1b4..6a82c3809b1086887eb129742460b5d302665174 100644
--- a/Numeric/functionSpace.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -9,7 +9,7 @@
 
 #include "GmshDefines.h"
 #include "GmshMessage.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
 
 static fullMatrix<double> generate1DMonomials(int order)
 {
@@ -733,7 +733,7 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
   }
 } 
 
-static void generate3dFaceClosure(functionSpace::clCont &closure, int order)
+static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
 {
   for (int iRotate = 0; iRotate < 3; iRotate++){
     for (int iSign = 1; iSign >= -1; iSign -= 2){
@@ -746,7 +746,7 @@ static void generate3dFaceClosure(functionSpace::clCont &closure, int order)
   }
 }
 
-static void generate2dEdgeClosure(functionSpace::clCont &closure, int order)
+static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order)
 {
   closure.clear();
   // first give edge nodes of the three edges in direct order
@@ -764,14 +764,14 @@ static void generate2dEdgeClosure(functionSpace::clCont &closure, int order)
   closure.push_back(c3); closure.push_back(c4); closure.push_back(c5);
 }
 
-std::map<int, functionSpace> functionSpaces::fs;
+std::map<int, polynomialBasis> polynomialBases::fs;
 
-const functionSpace &functionSpaces::find(int tag) 
+const polynomialBasis &polynomialBases::find(int tag) 
 {
-  std::map<int, functionSpace>::const_iterator it = fs.find(tag);
+  std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
   if (it != fs.end()) return it->second;
   
-  functionSpace F;
+  polynomialBasis F;
   
   switch (tag){
   case MSH_PNT:
@@ -801,77 +801,77 @@ const functionSpace &functionSpaces::find(int tag)
   case MSH_TRI_3 :
     F.monomials = generatePascalTriangle(1);
     F.points =    gmshGeneratePointsTriangle(1, false);
-    generate2dEdgeClosure (F.edgeClosure,1);
+    generate2dEdgeClosure(F.edgeClosure, 1);
     break;
   case MSH_TRI_6 :
     F.monomials = generatePascalTriangle(2);
     F.points =    gmshGeneratePointsTriangle(2, false);
-    generate2dEdgeClosure (F.edgeClosure,2);
+    generate2dEdgeClosure(F.edgeClosure, 2);
     break;
   case MSH_TRI_9 :
     F.monomials = generatePascalSerendipityTriangle(3);
     F.points =    gmshGeneratePointsTriangle(3, true);
-    generate2dEdgeClosure (F.edgeClosure,3);
+    generate2dEdgeClosure(F.edgeClosure, 3);
     break;
   case MSH_TRI_10 :
     F.monomials = generatePascalTriangle(3);
     F.points =    gmshGeneratePointsTriangle(3, false);
-    generate2dEdgeClosure (F.edgeClosure,3);
+    generate2dEdgeClosure(F.edgeClosure, 3);
     break;
   case MSH_TRI_12 :
     F.monomials = generatePascalSerendipityTriangle(4);
     F.points =    gmshGeneratePointsTriangle(4, true);
-    generate2dEdgeClosure (F.edgeClosure,4);
+    generate2dEdgeClosure(F.edgeClosure, 4);
     break;
   case MSH_TRI_15 :
     F.monomials = generatePascalTriangle(4);
     F.points =    gmshGeneratePointsTriangle(4, false);
-    generate2dEdgeClosure (F.edgeClosure,4);
+    generate2dEdgeClosure(F.edgeClosure, 4);
     break;
   case MSH_TRI_15I :
     F.monomials = generatePascalSerendipityTriangle(5);
     F.points =    gmshGeneratePointsTriangle(5, true);
-    generate2dEdgeClosure (F.edgeClosure,5);
+    generate2dEdgeClosure(F.edgeClosure, 5);
     break;
   case MSH_TRI_21 :
     F.monomials = generatePascalTriangle(5);
     F.points =    gmshGeneratePointsTriangle(5, false);
-    generate2dEdgeClosure (F.edgeClosure,5);
+    generate2dEdgeClosure(F.edgeClosure, 5);
     break;
   case MSH_TET_4 :
     F.monomials = generatePascalTetrahedron(1);
     F.points =    gmshGeneratePointsTetrahedron(1, false);
-    generate3dFaceClosure (F.faceClosure,1);
+    generate3dFaceClosure(F.faceClosure, 1);
     break;
   case MSH_TET_10 :
     F.monomials = generatePascalTetrahedron(2);
     F.points =    gmshGeneratePointsTetrahedron(2, false);
-    generate3dFaceClosure (F.faceClosure,2);
+    generate3dFaceClosure(F.faceClosure, 2);
     break;
   case MSH_TET_20 :
     F.monomials = generatePascalTetrahedron(3);
     F.points =    gmshGeneratePointsTetrahedron(3, false);
-    generate3dFaceClosure (F.faceClosure,3);
+    generate3dFaceClosure(F.faceClosure, 3);
     break;
   case MSH_TET_35 :
     F.monomials = generatePascalTetrahedron(4);
     F.points =    gmshGeneratePointsTetrahedron(4, false);
-    generate3dFaceClosure (F.faceClosure,4);
+    generate3dFaceClosure(F.faceClosure, 4);
     break;
   case MSH_TET_34 :
     F.monomials = generatePascalSerendipityTetrahedron(4);
     F.points =    gmshGeneratePointsTetrahedron(4, true);
-    generate3dFaceClosure (F.faceClosure,4);
+    generate3dFaceClosure(F.faceClosure, 4);
     break;
   case MSH_TET_52 :
     F.monomials = generatePascalSerendipityTetrahedron(5);
     F.points =    gmshGeneratePointsTetrahedron(5, true);
-    generate3dFaceClosure (F.faceClosure,5);
+    generate3dFaceClosure(F.faceClosure, 5);
     break;
   case MSH_TET_56 :
     F.monomials = generatePascalTetrahedron(5);
     F.points =    gmshGeneratePointsTetrahedron(5, false);
-    generate3dFaceClosure (F.faceClosure,5);
+    generate3dFaceClosure(F.faceClosure, 5);
     break;
   case MSH_QUA_4 :
     F.monomials = generatePascalQuad(1);
@@ -913,7 +913,7 @@ const functionSpace &functionSpaces::find(int tag)
     Msg::Error("Unknown function space %d: reverting to TET_4", tag);
     F.monomials = generatePascalTetrahedron(1);
     F.points =    gmshGeneratePointsTetrahedron(1, false);
-    generate3dFaceClosure (F.faceClosure,1);
+    generate3dFaceClosure(F.faceClosure, 1);
     break;
   }  
   F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
@@ -921,18 +921,18 @@ const functionSpace &functionSpaces::find(int tag)
   return fs[tag];
 }
 
-std::map<std::pair<int, int>, fullMatrix<double> > functionSpaces::injector;
+std::map<std::pair<int, int>, fullMatrix<double> > polynomialBases::injector;
 
-const fullMatrix<double> &functionSpaces::findInjector(int tag1, int tag2)
+const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
 {
   std::pair<int,int> key(tag1,tag2);
   std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it = injector.find(key);
   if (it != injector.end()) return it->second;
 
-  const functionSpace& fs1 = find(tag1);
-  const functionSpace& fs2 = find(tag2);
+  const polynomialBasis& fs1 = find(tag1);
+  const polynomialBasis& fs2 = find(tag2);
 
-  fullMatrix<double> inj(fs1.points.size1(),fs2.points.size1());
+  fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
   
   double sf[256];
   
diff --git a/Numeric/functionSpace.h b/Numeric/polynomialBasis.h
similarity index 83%
rename from Numeric/functionSpace.h
rename to Numeric/polynomialBasis.h
index 6f8ed6be74bb9caaad2ce0c17cca030be3b6074f..01b8f846c300c42300a0eaac3d59d15672e839c0 100644
--- a/Numeric/functionSpace.h
+++ b/Numeric/polynomialBasis.h
@@ -3,19 +3,19 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
-#ifndef _FUNCTION_SPACE_H_
-#define _FUNCTION_SPACE_H_
+#ifndef _POLYNOMIAL_BASIS_H_
+#define _POLYNOMIAL_BASIS_H_
 
 #include <math.h>
 #include <map>
 #include <vector>
 #include "fullMatrix.h"
 
-// presently thos function spaces are only for simplices
-// should be extended to other elements like quads and hexes
-struct functionSpace 
+// presently those function spaces are only for simplices and quads;
+// should be extended to other elements like hexes
+struct polynomialBasis 
 {
-  typedef  std::vector<std::vector<int> > clCont;
+  typedef std::vector<std::vector<int> > clCont;
   clCont faceClosure;
   clCont edgeClosure;
   fullMatrix<double> points;
@@ -23,11 +23,13 @@ struct functionSpace
   fullMatrix<double> coefficients;
   // for a given face/edge, with both a sign and a rotation,
   // give an ordered list of nodes on this face/edge
-  inline const std::vector<int> & getFaceClosure (int iFace, int iSign, int iRot)const{
-    return faceClosure[iFace+4*(iSign==1?0:1)+8*iRot];
+  inline const std::vector<int> &getFaceClosure(int iFace, int iSign, int iRot) const
+  {
+    return faceClosure[iFace + 4 * (iSign == 1 ? 0 : 1) + 8 * iRot];
   }
-  inline const std::vector<int> & getEdgeClosure (int iEdge, int iSign) const{
-    return edgeClosure[iSign == 1 ? iEdge : 3+iEdge];
+  inline const std::vector<int> &getEdgeClosure(int iEdge, int iSign) const
+  {
+    return edgeClosure[iSign == 1 ? iEdge : 3 + iEdge];
   }
   inline void evaluateMonomials(double u, double v, double w, double p[]) const 
   {
@@ -108,13 +110,13 @@ struct functionSpace
   }
 };
 
-class functionSpaces 
+class polynomialBases
 {
  private:
-  static std::map<int, functionSpace> fs;
+  static std::map<int, polynomialBasis> fs;
   static std::map<std::pair<int, int>, fullMatrix<double> > injector;
  public :
-  static const functionSpace &find(int);
+  static const polynomialBasis &find(int);
   static const fullMatrix<double> &findInjector(int, int);
 };
 
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index c82b6eb39635ee862cdbf05ee9f209259c494b54..4f1bf69a30a9356d790d16e588a9610bb9b5f106 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -93,7 +93,7 @@ bool PViewDataGModel::finalize()
 void PViewDataGModel::_addInterpolationMatricesForElement(MElement *e)
 {
   int type = e->getType();
-  const functionSpace *fs = e->getFunctionSpace();
+  const polynomialBasis *fs = e->getFunctionSpace();
   if(fs){
     if(e->getPolynomialOrder() > 1)
       setInterpolationMatrices(type, fs->coefficients, fs->monomials,
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 1a9f11a5f647a65f7390ef4ee05975541df13734..31d3d44b31880c8f34d062076c6d598a25fc31fe 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -6,7 +6,7 @@
 #include "PViewDataList.h"
 #include "GmshMessage.h"
 #include "GmshDefines.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
 #include "Numeric.h"
 #include "SmoothData.h"
 #include "Context.h"
@@ -784,7 +784,7 @@ void PViewDataList::setOrder2(int type)
   case TYPE_PRI: typeMSH = MSH_PRI_18; break;
   case TYPE_PYR: typeMSH = MSH_PYR_14; break;
   }
-  const functionSpace *fs = &functionSpaces::find(typeMSH);
+  const polynomialBasis *fs = &polynomialBases::find(typeMSH);
   if(!fs){
     Msg::Error("Could not find function space for element type %d", typeMSH);
     return;
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index 6a2088756ceee61cf064688b6c745ac7008f6e34..2e3e068f5b2ef777f4fe23069ee6f4302c1da016 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -15,13 +15,13 @@
 #include "fullMatrix.h"
 #include "Numeric.h"
 
-class convexCombinationTerm : public femTerm<double,double> {
+class convexCombinationTerm : public femTerm<double> {
  protected:
   const simpleFunction<double> *_k;
   const int _iField;
  public:
   convexCombinationTerm(GModel *gm, int iField, simpleFunction<double> *k)
-    : femTerm<double,double>(gm), _iField(iField), _k(k) {}
+    : femTerm<double>(gm), _iField(iField), _k(k) {}
   virtual int sizeOfR(SElement *se) const
   {
     return se->getMeshElement()->getNumVertices();
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index 52521de11778c88648c7ece82fd7b84aeb5f0c44..8f9ce458f9fa15933918f96101b5363fd203b2b3 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -14,7 +14,7 @@
 #include "fullMatrix.h"
 #include "Numeric.h"
 
-class crossConfTerm : public femTerm<double, double> {
+class crossConfTerm : public femTerm<double> {
  protected:
   const simpleFunction<double> *_diffusivity;
   const int _iFieldR;
@@ -22,7 +22,7 @@ class crossConfTerm : public femTerm<double, double> {
  public:
   crossConfTerm(GModel *gm, int iFieldR, int iFieldC, 
                 simpleFunction<double> *diffusivity)
-    : femTerm<double, double>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC), 
+    : femTerm<double>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC), 
       _diffusivity(diffusivity) {}
   virtual int sizeOfR(SElement *se) const 
   {
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 30c269ab42fc75f257b5f0ad4a11c6a12c0463ef..de8bfe3ddb0ed0699672a47a8bea0153fa8e0424 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -1,6 +1,6 @@
 #include "dgGroupOfElements.h"
 #include "MElement.h"
-#include "functionSpace.h"
+#include "polynomialBasis.h"
 #include "Numeric.h"
 #include "MTriangle.h"
 #include "MLine.h"
@@ -20,9 +20,9 @@ static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
 }
 
 static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
-      const functionSpace *fsFace,
+      const polynomialBasis *fsFace,
       const fullMatrix<double> &intgFace,
-      const functionSpace *fsElement,
+      const polynomialBasis *fsElement,
       const std::vector <int> *closure) {
   int npts=intgFace.size1();
   fullMatrix<double> *m = new fullMatrix<double> (npts, 4);
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 6277359e3412fb8dc8addcff47660360c252c59a..f00b80a7b4802b3d9c0f365e8d14456b526cf6d2 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -16,7 +16,7 @@
 class MElement;
 class MFace;
 class MEdge;
-class functionSpace;
+class polynomialBasis;
 
 class dgElement {
   MElement *_element;
@@ -42,7 +42,7 @@ class dgGroupOfElements {
   // the ONLY function space that is used to 
   // inerpolated the fields (may be different to the 
   // one of the elements)
-  const functionSpace &_fs;
+  const polynomialBasis &_fs;
   // Ni integration points, matrix of size Ni x 4 (u,v,w,weight)
   fullMatrix<double> *_integration;
   // collocation matrix that maps vertices to integration points.
@@ -108,9 +108,9 @@ class dgGroupOfFaces {
   void addFace(const MFace &topoFace, int iElLeft, int iElRight);
   void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
   void computeFaceNormals();
-  // Two functionSpaces for left and right elements
+  // Two polynomialBases for left and right elements
   // the group has always the same types for left and right
-  const functionSpace *_fsLeft,*_fsRight, *_fsFace;
+  const polynomialBasis *_fsLeft,*_fsRight, *_fsFace;
   // N elements in the group
   std::vector<int>_left, _right;
   std::vector<MElement *>_faces;
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 9ff945e0ff3bcec912af9da80f0d5b56fe6b63e7..fd2d2fccb46ff77f0837dee5183a5bdda69eeca2 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -8,6 +8,7 @@
 
 #include <vector>
 #include <string>
+#include <complex>
 #include <map>
 #include "MVertex.h"
 #include "linearSystem.h"
@@ -43,21 +44,42 @@ class Dof{
   }
 };
 
-template<class dataVec, class dataMat>
+template<class T> struct dofTraits
+{
+  typedef T VecType;
+  typedef T MatType;
+};
+
+template<> struct dofTraits<fullVector<double> >
+{
+  typedef fullVector<double> VecType;
+  typedef fullMatrix<double> MatType;
+};
+
+template<> struct dofTraits<fullVector<std::complex<double> > >
+{
+  typedef fullVector<std::complex<double> > VecType;
+  typedef fullMatrix<std::complex<double> > MatType;
+};
+
+template<class T>
 class DofAffineConstraint{
  public:
-  std::vector<std::pair<Dof, dataMat> > linear;
-  dataVec shift;
+  std::vector<std::pair<Dof, typename dofTraits<T>::MatType> > linear;
+  typename dofTraits<T>::VecType shift;
 };
   
 // data=float, double, complex<double>, smallVec<double>, smallVec<complex<double> >...
-template <class dataVec, class dataMat>
+template <class T>
 class dofManager{
  private:
+  typedef typename dofTraits<T>::VecType dataVec;
+  typedef typename dofTraits<T>::MatType dataMat;
+
   // general affine constraint on sub-blocks, treated by adding
   // equations:
   //   dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
-  std::map<std::pair<Dof, dataMat>, DofAffineConstraint<dataVec, dataMat> > constraints;
+  std::map<std::pair<Dof, dataMat>, DofAffineConstraint<dataVec> > constraints;
     
   // fixations on full blocks, treated by eliminating equations:
   //   DofVec = dataVec    
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index ddfade0e246ddda22fac6f7a45714aadab383ef5..8f771bdbbcfe3b67a44bba01692e44124c7d0083 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -9,7 +9,7 @@
 
 #include <slepceps.h>
 
-eigenSolver::eigenSolver(dofManager<double, double> *manager, std::string A, 
+eigenSolver::eigenSolver(dofManager<double> *manager, std::string A, 
                          std::string B) : _A(0), _B(0)
 {
   if(A.size()){
diff --git a/Solver/eigenSolver.h b/Solver/eigenSolver.h
index 0ca589d14e5c9d8907d96ad908affb0ccb953ea6..ee3e8d9ff8910c0bc133386ea2a79f437cbdd837 100644
--- a/Solver/eigenSolver.h
+++ b/Solver/eigenSolver.h
@@ -24,7 +24,7 @@ class eigenSolver{
   std::vector<std::vector<std::complex<double> > > _eigenVectors;
   void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
  public:
-  eigenSolver(dofManager<double, double> *manager, std::string A, 
+  eigenSolver(dofManager<double> *manager, std::string A, 
               std::string B="");
   void solve(int numEigenValues=0, std::string which="");
   int getNumEigenValues(){ return _eigenValues.size(); }
@@ -38,7 +38,7 @@ class eigenSolver{
  private:
   std::vector<std::complex<double> > _dummy;
  public:
-  eigenSolver(dofManager<double, double> *manager, std::string A, 
+  eigenSolver(dofManager<double> *manager, std::string A, 
               std::string B=""){}
   void solve(int numEigenValues=0, std::string which="")
   {
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 73f99fe40c248f6cfa18c061b202913cbd773b81..4358067e653f4476ebbec577e5b517ff4158be90 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -59,7 +59,7 @@ PView* elasticitySolver::buildDisplacementView  (const std::string &postFileName
 #endif
 
 
-static double vonMises(dofManager<double,double> *a, MElement *e, 
+static double vonMises(dofManager<double> *a, MElement *e, 
                        double u, double v, double w, 
                        double E, double nu, int _tag)
 {
@@ -195,7 +195,7 @@ void elasticitySolver::solve()
   linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
   lsys->setNoisy(2);
 #endif
-  pAssembler = new dofManager<double, double>(lsys);
+  pAssembler = new dofManager<double>(lsys);
 
   std::map<int, std::vector<GEntity*> > groups[4];
   pModel->getPhysicalGroups(groups);
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 10fa60f3429b42069bd697c63afac87ed6e211a8..c0fedc0bb25cdacd6079b3bffcf2764701081061 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -29,7 +29,7 @@ class elasticitySolver{
  protected:
   GModel *pModel;
   int _dim, _tag;
-  dofManager<double, double> *pAssembler;
+  dofManager<double> *pAssembler;
   // young modulus and poisson coefficient per physical
   std::vector<elasticField> elasticFields;
   // imposed nodal forces
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index 1c0292cf884b8a18ec9315963f3d0958b77ae42f..09e32429018ecc1662c562c524b11daac8211b78 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -12,7 +12,7 @@
 #include "SElement.h"
 #include "fullMatrix.h"
 
-class elasticityTerm : public femTerm<double, double> {
+class elasticityTerm : public femTerm<double> {
  protected:
   double _E, _nu;
   int _iFieldR, _iFieldC;
@@ -49,9 +49,9 @@ class elasticityTerm : public femTerm<double, double> {
   }
  public:
  elasticityTerm(GModel *gm, double E, double nu, int fieldr, int fieldc)
-   : femTerm<double, double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldc) {}
+   : femTerm<double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldc) {}
  elasticityTerm(GModel *gm, double E, double nu, int fieldr)
-   : femTerm<double, double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldr) {}
+   : femTerm<double>(gm), _E(E), _nu(nu), _iFieldR(fieldr),_iFieldC(fieldr) {}
   void setVector(const SVector3 &f) { _volumeForce = f; }
   void elementMatrix(SElement *se, fullMatrix<double> &m) const;
   void elementVector(SElement *se, fullVector<double> &m) const;
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 5c35bbe4c00a4ee1397d63d765d7a3fb438b2c74..3536058089a0b5867eb88b4d00a24bbee9c2cd94 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -18,8 +18,11 @@
 
 // a nodal finite element term : variables are always defined at nodes
 // of the mesh
-template<class dataVec, class dataMat>
+template<class T>
 class femTerm {
+ private:
+  typedef typename dofTraits<T>::VecType dataVec;
+  typedef typename dofTraits<T>::MatType dataMat;
  protected:
   GModel *_gm;
  public:
@@ -46,7 +49,7 @@ class femTerm {
 
   // add the contribution from all the elements in the intersection
   // of two element groups L and C
-  void addToMatrix(dofManager<dataVec, dataMat> &dm,
+  void addToMatrix(dofManager<dataVec> &dm,
                    groupOfElements &L,
                    groupOfElements &C) const
   {
@@ -61,7 +64,7 @@ class femTerm {
   }
 
   // add the contribution from a single element to the dof manager
-  void addToMatrix(dofManager<dataVec, dataMat> &dm, SElement *se) const
+  void addToMatrix(dofManager<dataVec> &dm, SElement *se) const
   {
     const int nbR = sizeOfR(se);
     const int nbC = sizeOfC(se);
@@ -69,7 +72,7 @@ class femTerm {
     elementMatrix(se, localMatrix);
     addToMatrix(dm, localMatrix, se);
   }
-  void addToMatrix(dofManager<dataVec, dataMat> &dm,
+  void addToMatrix(dofManager<dataVec> &dm,
                    fullMatrix<dataMat> &localMatrix,
                    SElement *se) const
   {
@@ -116,7 +119,7 @@ class femTerm {
   }
   void dirichletNodalBC(int physical, int dim, int comp, int field,
                         const simpleFunction<dataVec> &e,
-                        dofManager<dataVec, dataMat> &dm)
+                        dofManager<dataVec> &dm)
   {
     std::vector<MVertex *> v;
     GModel *m = _gm;
@@ -126,7 +129,7 @@ class femTerm {
   }
   void neumannNodalBC(int physical, int dim, int comp, int field,
                       const simpleFunction<dataVec> &fct,
-                      dofManager<dataVec, dataMat> &dm)
+                      dofManager<dataVec> &dm)
   {
     std::map<int, std::vector<GEntity*> > groups[4];
     GModel *m = _gm;
@@ -161,7 +164,7 @@ class femTerm {
     }
   }
 
-  void addToRightHandSide(dofManager<dataVec, dataMat> &dm, groupOfElements &C) const
+  void addToRightHandSide(dofManager<dataVec> &dm, groupOfElements &C) const
   {
     groupOfElements::elementContainer::const_iterator it = C.begin();
     for ( ; it != C.end(); ++it){
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 09af7cf50fed1aec619fdb8c2dbae6b24b5af760..9d4401e203a3df2682fa44803abfc7ddc31cf06f 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -17,7 +17,7 @@
 
 // \nabla \cdot k \nabla U - a U 
 template<class scalar>
-class helmholtzTerm : public femTerm<scalar, scalar> {
+class helmholtzTerm : public femTerm<scalar> {
  protected:
   const simpleFunction<scalar> *_k, *_a;
   const int _iFieldR;
@@ -25,7 +25,7 @@ class helmholtzTerm : public femTerm<scalar, scalar> {
  public:
   helmholtzTerm(GModel *gm, int iFieldR, int iFieldC, simpleFunction<scalar> *k,
                 simpleFunction<scalar> *a) 
-    : femTerm<scalar, scalar>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
+    : femTerm<scalar>(gm), _iFieldR(iFieldR), _iFieldC(iFieldC),
       _k(k), _a(a) {}
   // one dof per vertex (nodal fem)
   virtual int sizeOfR(SElement *se) const 
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index e1d4b7b93cc7f106fa35c23d718f36ba6ab5d2b9..0f5c09f5795a7ccd5705658a0c73e11fc950df2e 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -120,7 +120,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
   std::map<MVertex*,SPoint2> solution;
   // PARAMETRIZE ---------------------------------
   for (int step =0 ; step<2 ; step++){
-    dofManager<double, double> myAssembler(_lsys);
+    dofManager<double> myAssembler(_lsys);
     for(std::map<MVertex*,SPoint2>::iterator it = level.coordinates.begin();
 	it != level.coordinates.end(); ++it){
       MVertex *v = it->first;