diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index c204d9f390cf1e03816bdd08a0a0881ff577ddbd..3458b4ce7a39f0d815fb6aabb0ca1c38223fc4fa 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -47,7 +47,7 @@ double MHexahedron::angleShapeMeasure()
      std::vector<MVertex*> vv;
      vv.push_back(getFace(i).getVertex(0));
      vv.push_back(getFace(i).getVertex(1));
-     vv.push_back(getFace(i).getVertex(2)); 
+     vv.push_back(getFace(i).getVertex(2));
      vv.push_back(getFace(i).getVertex(3));
      // MVertex *v0 = new MVertex(0, 0, 0); vv.push_back(v0);
      // MVertex *v1 = new MVertex(1., 0, 0);vv.push_back(v1);
@@ -63,7 +63,7 @@ double MHexahedron::angleShapeMeasure()
      //printf("angle max =%g min =%g \n", angleMax*180/M_PI, angleMin*180/M_PI);
    }
    zeta = 1.-std::max((angleMax-0.5*M_PI)/(0.5*M_PI),(0.5*M_PI-angleMin)/(0.5*M_PI));
-   return zeta; 
+   return zeta;
 #else
    return 1.;
 #endif
@@ -189,30 +189,12 @@ int MHexahedronN::getNumEdgesRep()
 {
   return 12 * CTX::instance()->mesh.numSubEdges;
 }
- 
-const nodalBasis* MHexahedron::getFunctionSpace(int o) const
-{
-  int order = (o == -1) ? getPolynomialOrder() : o;
 
-  int nv = getNumVolumeVertices();
+const nodalBasis* MHexahedron::getFunctionSpace(int order) const
+{
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
 
-  if ((nv == 0) && (o == -1)) {
-    switch (order) {
-    case 0: return BasisFactory::getNodalBasis(MSH_HEX_1);
-    case 1: return BasisFactory::getNodalBasis(MSH_HEX_8);
-    case 2: return BasisFactory::getNodalBasis(MSH_HEX_20);
-    case 3: return BasisFactory::getNodalBasis(MSH_HEX_32);
-    case 4: return BasisFactory::getNodalBasis(MSH_HEX_44);
-    case 5: return BasisFactory::getNodalBasis(MSH_HEX_56);
-    case 6: return BasisFactory::getNodalBasis(MSH_HEX_68);
-    case 7: return BasisFactory::getNodalBasis(MSH_HEX_80);
-    case 8: return BasisFactory::getNodalBasis(MSH_HEX_92);
-    case 9: return BasisFactory::getNodalBasis(MSH_HEX_104);
-    default: Msg::Error("Order %d hex function space not implemented", order); break;
-    }
-  }
-  else {
-    switch (order) {
+  switch (order) {
     case 0: return BasisFactory::getNodalBasis(MSH_HEX_1);
     case 1: return BasisFactory::getNodalBasis(MSH_HEX_8);
     case 2: return BasisFactory::getNodalBasis(MSH_HEX_27);
@@ -224,34 +206,15 @@ const nodalBasis* MHexahedron::getFunctionSpace(int o) const
     case 8: return BasisFactory::getNodalBasis(MSH_HEX_729);
     case 9: return BasisFactory::getNodalBasis(MSH_HEX_1000);
     default: Msg::Error("Order %d hex function space not implemented", order); break;
-    }
   }
-  return 0;
+  return NULL;
 }
 
-const JacobianBasis* MHexahedron::getJacobianFuncSpace(int o) const
+const JacobianBasis* MHexahedron::getJacobianFuncSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
-
-  int nv = getNumVolumeVertices();
+  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
 
-  if ((nv == 0) && (o == -1)) {
-    switch (order) {
-    case 0: return BasisFactory::getJacobianBasis(MSH_HEX_1);
-    case 1: return BasisFactory::getJacobianBasis(MSH_HEX_8);
-    case 2: return BasisFactory::getJacobianBasis(MSH_HEX_20);
-    case 3: return BasisFactory::getJacobianBasis(MSH_HEX_32);
-    case 4: return BasisFactory::getJacobianBasis(MSH_HEX_44);
-    case 5: return BasisFactory::getJacobianBasis(MSH_HEX_56);
-    case 6: return BasisFactory::getJacobianBasis(MSH_HEX_68);
-    case 7: return BasisFactory::getJacobianBasis(MSH_HEX_80);
-    case 8: return BasisFactory::getJacobianBasis(MSH_HEX_92);
-    case 9: return BasisFactory::getJacobianBasis(MSH_HEX_104);
-    default: Msg::Error("Order %d hex incomplete Jacobian function space not implemented", order); break;
-    }
-  }
-  else {
-    switch (order) {
+  switch (order) {
     case 0: return BasisFactory::getJacobianBasis(MSH_HEX_1);
     case 1: return BasisFactory::getJacobianBasis(MSH_HEX_8);
     case 2: return BasisFactory::getJacobianBasis(MSH_HEX_27);
@@ -263,9 +226,8 @@ const JacobianBasis* MHexahedron::getJacobianFuncSpace(int o) const
     case 8: return BasisFactory::getJacobianBasis(MSH_HEX_729);
     case 9: return BasisFactory::getJacobianBasis(MSH_HEX_1000);
     default: Msg::Error("Order %d hex Jacobian function space not implemented", order); break;
-    }
   }
-  return 0;
+  return NULL;
 }
 
 
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 508bf117d7641ec1d98cd441b6b66914ca59045e..49456e2a6a16a0835d5922164f71ead5303880ad 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -9,7 +9,7 @@
 #include "Context.h"
 
 int MPrism::getVolumeSign()
-{ 
+{
   double mat[3][3];
   mat[0][0] = _v[1]->x() - _v[0]->x();
   mat[0][1] = _v[2]->x() - _v[0]->x();
@@ -32,52 +32,44 @@ void MPrism::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *pts = getGQPriPts(pOrder);
 }
 
-const nodalBasis* MPrism::getFunctionSpace(int o) const
+const nodalBasis* MPrism::getFunctionSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
-
-  int nv = getNumVolumeVertices();
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
 
-  if ((nv == 0) && (o == -1)) {
-    switch (order) {
+  switch (order) {
     case 0: return BasisFactory::getNodalBasis(MSH_PRI_1);
     case 1: return BasisFactory::getNodalBasis(MSH_PRI_6);
     case 2: return BasisFactory::getNodalBasis(MSH_PRI_18);
+    case 3: return BasisFactory::getNodalBasis(MSH_PRI_40);
+    case 4: return BasisFactory::getNodalBasis(MSH_PRI_75);
+    case 5: return BasisFactory::getNodalBasis(MSH_PRI_126);
+    case 6: return BasisFactory::getNodalBasis(MSH_PRI_196);
+    case 7: return BasisFactory::getNodalBasis(MSH_PRI_288);
+    case 8: return BasisFactory::getNodalBasis(MSH_PRI_405);
+    case 9: return BasisFactory::getNodalBasis(MSH_PRI_550);
     default: Msg::Error("Order %d prism function space not implemented", order);
-    }
-  }
-  else {
-    switch (order) {
-    case 0: return BasisFactory::getNodalBasis(MSH_PRI_1);
-    case 1: return BasisFactory::getNodalBasis(MSH_PRI_6);
-    case 2: return BasisFactory::getNodalBasis(MSH_PRI_18);
-    default: Msg::Error("Order %d prism function space not implemented", order);
-    }
   }
-  return 0;
+  return NULL;
 }
 
-const JacobianBasis* MPrism::getJacobianFuncSpace(int o) const
+const JacobianBasis* MPrism::getJacobianFuncSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
+  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
 
-  int nv = getNumVolumeVertices();
-  
-  if ((nv == 0) && (o == -1)) {
-    switch (order) {
+  switch (order) {
+    case 0: return BasisFactory::getJacobianBasis(MSH_PRI_1);
     case 1: return BasisFactory::getJacobianBasis(MSH_PRI_6);
     case 2: return BasisFactory::getJacobianBasis(MSH_PRI_18);
+    case 3: return BasisFactory::getJacobianBasis(MSH_PRI_40);
+    case 4: return BasisFactory::getJacobianBasis(MSH_PRI_75);
+    case 5: return BasisFactory::getJacobianBasis(MSH_PRI_126);
+    case 6: return BasisFactory::getJacobianBasis(MSH_PRI_196);
+    case 7: return BasisFactory::getJacobianBasis(MSH_PRI_288);
+    case 8: return BasisFactory::getJacobianBasis(MSH_PRI_405);
+    case 9: return BasisFactory::getJacobianBasis(MSH_PRI_550);
     default: Msg::Error("Order %d prism function space not implemented", order);
-    }
-  }
-  else { 
-    switch (order) {
-    case 1: return BasisFactory::getJacobianBasis(MSH_PRI_6);
-    case 2: return BasisFactory::getJacobianBasis(MSH_PRI_18);
-    default: Msg::Error("Order %d prism function space not implemented", order);
-    }
   }
-  return 0;
+  return NULL;
 }
 
 double MPrism::getInnerRadius()
@@ -123,7 +115,7 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
     }
     else {
       MVertex *v3 = _v[faces_prism(ithFace, 3)];
-      if (v0 == face.getVertex(0) && v1 == face.getVertex(1) && 
+      if (v0 == face.getVertex(0) && v1 == face.getVertex(1) &&
           v2 == face.getVertex(2) && v3 == face.getVertex(3)){
         sign = 1; rot = 0; return;
       }
@@ -135,7 +127,7 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
           v2 == face.getVertex(0) && v3 == face.getVertex(1)){
         sign = 1; rot = 2; return;
       }
-      if (v0 == face.getVertex(3) && v1 == face.getVertex(0) && 
+      if (v0 == face.getVertex(3) && v1 == face.getVertex(0) &&
           v2 == face.getVertex(1) && v3 == face.getVertex(2)){
         sign = 1; rot = 3; return;
       }
@@ -147,7 +139,7 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
           v2 == face.getVertex(3) && v3 == face.getVertex(2)){
         sign = -1; rot = 1; return;
       }
-      if (v0 == face.getVertex(2) && v1 == face.getVertex(1) && 
+      if (v0 == face.getVertex(2) && v1 == face.getVertex(1) &&
           v2 == face.getVertex(0) && v3 == face.getVertex(3)){
         sign = -1; rot = 2; return;
       }
@@ -156,7 +148,7 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
         sign = -1; rot = 3; return;
       }
     }
-      
+
   }
   Msg::Error("Could not get face information for prism %d", getNum());
 }
@@ -451,4 +443,4 @@ int MPrism15::getNumFacesRep() {
 
 int MPrism18::getNumFacesRep() {
   return 4 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
-}
\ No newline at end of file
+}
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index bd553a6e60f9b557a84c41f6cbc5cd1d3be7edc8..fb41b100d7e4eeedbf587f9dfb6a41ffbd8a8b64 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -33,29 +33,11 @@ void MPyramid::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *pts = getGQPyrPts(pOrder);
 }
 
-const nodalBasis* MPyramid::getFunctionSpace(int o) const
+const nodalBasis* MPyramid::getFunctionSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
-
-  int nv = getNumVolumeVertices();
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
 
-  if ((nv == 0) && (o == -1)) {
-    switch (order) {
-    case 1: return BasisFactory::getNodalBasis(MSH_PYR_5);
-    case 2: return BasisFactory::getNodalBasis(MSH_PYR_14);
-    case 0: return BasisFactory::getNodalBasis(MSH_PYR_1);
-    case 3: return BasisFactory::getNodalBasis(MSH_PYR_21);
-    case 4: return BasisFactory::getNodalBasis(MSH_PYR_29);
-    case 5: return BasisFactory::getNodalBasis(MSH_PYR_37);
-    case 6: return BasisFactory::getNodalBasis(MSH_PYR_45);
-    case 7: return BasisFactory::getNodalBasis(MSH_PYR_53);
-    case 8: return BasisFactory::getNodalBasis(MSH_PYR_61);
-    case 9: return BasisFactory::getNodalBasis(MSH_PYR_69);
-    default: Msg::Error("Order %d pyramid function space not implemented", order);
-    }
-  }
-  else {
-    switch (order) {
+  switch (order) {
     case 0: return BasisFactory::getNodalBasis(MSH_PYR_1);
     case 1: return BasisFactory::getNodalBasis(MSH_PYR_5);
     case 2: return BasisFactory::getNodalBasis(MSH_PYR_14);
@@ -67,13 +49,13 @@ const nodalBasis* MPyramid::getFunctionSpace(int o) const
     case 8: return BasisFactory::getNodalBasis(MSH_PYR_285);
     case 9: return BasisFactory::getNodalBasis(MSH_PYR_385);
     default: Msg::Error("Order %d pyramid function space not implemented", order);
-    }
   }
-  return 0;
+  return NULL;
 }
 
 const JacobianBasis* MPyramid::getJacobianFuncSpace(int o) const
 {
+  // FIXME add other order and see MPyramid::getFunctionSpace for 'design'
   int order = (o == -1) ? getPolynomialOrder() : o;
 
   switch (order) {
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 0f6b54877c77af7848f6e4f49cf5e1b9e9f503c5..221208c787d0c67b2b72957dd852715ebd86f852 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -17,27 +17,10 @@
 
 #define SQU(a)      ((a)*(a))
 
-const nodalBasis* MQuadrangle::getFunctionSpace(int o) const
+const nodalBasis* MQuadrangle::getFunctionSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
-
-  int nf = getNumFaceVertices();
-
-  if ((nf == 0) && (o == -1)) {
-    switch (order) {
-      case 0: return BasisFactory::getNodalBasis(MSH_QUA_1);
-      case 1: return BasisFactory::getNodalBasis(MSH_QUA_4);
-      case 2: return BasisFactory::getNodalBasis(MSH_QUA_8);
-      case 3: return BasisFactory::getNodalBasis(MSH_QUA_12);
-      case 4: return BasisFactory::getNodalBasis(MSH_QUA_16I);
-      case 5: return BasisFactory::getNodalBasis(MSH_QUA_20);
-      case 6: return BasisFactory::getNodalBasis(MSH_QUA_24);
-      case 7: return BasisFactory::getNodalBasis(MSH_QUA_28);
-      case 8: return BasisFactory::getNodalBasis(MSH_QUA_32);
-      case 9: return BasisFactory::getNodalBasis(MSH_QUA_36I);
-      case 10: return BasisFactory::getNodalBasis(MSH_QUA_40);
-    }
-  }
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
+
   switch (order) {
     case 0: return BasisFactory::getNodalBasis(MSH_QUA_1);
     case 1: return BasisFactory::getNodalBasis(MSH_QUA_4);
@@ -52,29 +35,12 @@ const nodalBasis* MQuadrangle::getFunctionSpace(int o) const
     case 10: return BasisFactory::getNodalBasis(MSH_QUA_121);
     default: Msg::Error("Order %d quadrangle function space not implemented", order);
   }
-  return 0;
+  return NULL;
 }
 
-const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int o) const
+const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
-
-  int nf = getNumFaceVertices();
-
-  if ((nf == 0) && (o == -1)) {
-    switch (order) {
-      case 1: return BasisFactory::getJacobianBasis(MSH_QUA_4);
-      case 2: return BasisFactory::getJacobianBasis(MSH_QUA_8);
-      case 3: return BasisFactory::getJacobianBasis(MSH_QUA_12);
-      case 4: return BasisFactory::getJacobianBasis(MSH_QUA_16I);
-      case 5: return BasisFactory::getJacobianBasis(MSH_QUA_20);
-      case 6: return BasisFactory::getJacobianBasis(MSH_QUA_24);
-      case 7: return BasisFactory::getJacobianBasis(MSH_QUA_28);
-      case 8: return BasisFactory::getJacobianBasis(MSH_QUA_32);
-      case 9: return BasisFactory::getJacobianBasis(MSH_QUA_36I);
-      case 10: return BasisFactory::getJacobianBasis(MSH_QUA_40);
-    }
-  }
+  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
   switch (order) {
     case 1: return BasisFactory::getJacobianBasis(MSH_QUA_4);
     case 2: return BasisFactory::getJacobianBasis(MSH_QUA_9);
@@ -234,7 +200,7 @@ void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 double  MQuadrangle::etaShapeMeasure()
 {
   double AR = 1;//(minEdge()/maxEdge());
-  
+
   SVector3 v01 (_v[1]->x()-_v[0]->x(),_v[1]->y()-_v[0]->y(),_v[1]->z()-_v[0]->z());
   SVector3 v12 (_v[2]->x()-_v[1]->x(),_v[2]->y()-_v[1]->y(),_v[2]->z()-_v[1]->z());
   SVector3 v23 (_v[3]->x()-_v[2]->x(),_v[3]->y()-_v[2]->y(),_v[3]->z()-_v[2]->z());
@@ -268,12 +234,12 @@ double  MQuadrangle::etaShapeMeasure()
 }
 
 /// a shape measure for quadrangles
-/// assume (for now) 2D elements -- 
+/// assume (for now) 2D elements --
 ///  sf = (1 \pm xi) (1 \pm eta) / 4
 ///  dsf_xi  =  \pm (1 \pm  eta) / 4
-///             1 + eta , -(1+eta) , -(1-eta), 1-eta  
+///             1 + eta , -(1+eta) , -(1-eta), 1-eta
 ///  dsf_eta =  \pm (1 \pm  xi)  / 4
-///             1 + xi , 1 - xi ,  -(1-xi), -(1+xi)    
+///             1 + xi , 1 - xi ,  -(1-xi), -(1+xi)
 double MQuadrangle::gammaShapeMeasure(){
   return etaShapeMeasure();
 }
@@ -366,10 +332,10 @@ double MQuadrangle::getInnerRadius()
   }
   return R;
 #else // HAVE_LAPACK
-  // Default implementation. Not sure that the following give exactly 
+  // Default implementation. Not sure that the following give exactly
   // the same value as the HAVE_LAPACK case !
   // but same value for a square
-  
+
   // Mid-point of each edge of the quadrangle
   SPoint3 A(_v[0]->x()+_v[1]->x(),_v[0]->y()+_v[1]->y(),_v[0]->z()+_v[1]->z());
   SPoint3 B(_v[1]->x()+_v[2]->x(),_v[1]->y()+_v[2]->y(),_v[1]->z()+_v[2]->z());
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 2fd31d1f447b1ccd500ff782316494301ec1e9e7..6b4e3c7b16b012a637198abf88de6f9b4d2bf418 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -100,30 +100,11 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
   sys3x3(mat, b, uvw, &det);
 }
 
-const nodalBasis* MTetrahedron::getFunctionSpace(int o) const
+const nodalBasis* MTetrahedron::getFunctionSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
 
-  int nv = getNumVolumeVertices();
-
-  if ((nv == 0) && (o == -1)) {
-    switch (order) {
-    case 0: return BasisFactory::getNodalBasis(MSH_TET_1);
-    case 1: return BasisFactory::getNodalBasis(MSH_TET_4);
-    case 2: return BasisFactory::getNodalBasis(MSH_TET_10);
-    case 3: return BasisFactory::getNodalBasis(MSH_TET_16); // not just nv==0
-    case 4: return BasisFactory::getNodalBasis(MSH_TET_22);
-    case 5: return BasisFactory::getNodalBasis(MSH_TET_28);
-    case 6: return BasisFactory::getNodalBasis(MSH_TET_34);
-    case 7: return BasisFactory::getNodalBasis(MSH_TET_40);
-    case 8: return BasisFactory::getNodalBasis(MSH_TET_46);
-    case 9: return BasisFactory::getNodalBasis(MSH_TET_52);
-    case 10: return BasisFactory::getNodalBasis(MSH_TET_58);
-    default: Msg::Error("Order %d tetrahedron function space not implemented", order);
-    }
-  }
-  else {
-    switch (order) {
+  switch (order) {
     case 0: return BasisFactory::getNodalBasis(MSH_TET_1);
     case 1: return BasisFactory::getNodalBasis(MSH_TET_4);
     case 2: return BasisFactory::getNodalBasis(MSH_TET_10);
@@ -136,34 +117,15 @@ const nodalBasis* MTetrahedron::getFunctionSpace(int o) const
     case 9: return BasisFactory::getNodalBasis(MSH_TET_220);
     case 10: return BasisFactory::getNodalBasis(MSH_TET_286);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
-    }
   }
-  return 0;
+  return NULL;
 }
 
-const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int o) const
+const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
+  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
 
-  int nv = getNumVolumeVertices();
-
-  if ((nv == 0) && (o == -1)) {
-    switch (order) {
-    case 1: return BasisFactory::getJacobianBasis(MSH_TET_4);
-    case 2: return BasisFactory::getJacobianBasis(MSH_TET_10);
-    case 3: return BasisFactory::getJacobianBasis(MSH_TET_16); // not just nv==0
-    case 4: return BasisFactory::getJacobianBasis(MSH_TET_22);
-    case 5: return BasisFactory::getJacobianBasis(MSH_TET_28);
-    case 6: return BasisFactory::getJacobianBasis(MSH_TET_34);
-    case 7: return BasisFactory::getJacobianBasis(MSH_TET_40);
-    case 8: return BasisFactory::getJacobianBasis(MSH_TET_46);
-    case 9: return BasisFactory::getJacobianBasis(MSH_TET_52);
-    case 10: return BasisFactory::getJacobianBasis(MSH_TET_58);
-    default: Msg::Error("Order %d tetrahedron function space not implemented", order);
-    }
-  }
-  else {
-    switch (order) {
+  switch (order) {
     case 1: return BasisFactory::getJacobianBasis(MSH_TET_4);
     case 2: return BasisFactory::getJacobianBasis(MSH_TET_10);
     case 3: return BasisFactory::getJacobianBasis(MSH_TET_20);
@@ -175,9 +137,8 @@ const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int o) const
     case 9: return BasisFactory::getJacobianBasis(MSH_TET_220);
     case 10: return BasisFactory::getJacobianBasis(MSH_TET_286);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
-    }
   }
-  return 0;
+  return NULL;
 }
 
 int MTetrahedron10::getNumEdgesRep(){ return 6 * CTX::instance()->mesh.numSubEdges; }
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 55ac74da6749c1de9c3cb717231fb77a7c0c8ea0..121307e6aca2c049f4df1a78b6171cabed32cbfc 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -119,30 +119,11 @@ double MTriangle::gammaShapeMeasure()
    uvw[2] = 0.;
 }
 
-const nodalBasis* MTriangle::getFunctionSpace(int o) const
+const nodalBasis* MTriangle::getFunctionSpace(int order) const
 {
-  int order = (o == -1) ? getPolynomialOrder() : o;
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
 
-  int nf = getNumFaceVertices();
-
-  if ((nf == 0) && (o == -1)) {
-    switch (order) {
-    case 0: return BasisFactory::getNodalBasis(MSH_TRI_1);
-    case 1: return BasisFactory::getNodalBasis(MSH_TRI_3);
-    case 2: return BasisFactory::getNodalBasis(MSH_TRI_6);
-    case 3: return BasisFactory::getNodalBasis(MSH_TRI_9);
-    case 4: return BasisFactory::getNodalBasis(MSH_TRI_12);
-    case 5: return BasisFactory::getNodalBasis(MSH_TRI_15I);
-    case 6: return BasisFactory::getNodalBasis(MSH_TRI_18);
-    case 7: return BasisFactory::getNodalBasis(MSH_TRI_21I);
-    case 8: return BasisFactory::getNodalBasis(MSH_TRI_24);
-    case 9: return BasisFactory::getNodalBasis(MSH_TRI_27);
-    case 10: return BasisFactory::getNodalBasis(MSH_TRI_30);
-    default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
-    }
-  }
-  else {
-    switch (order) {
+  switch (order) {
     case 0: return BasisFactory::getNodalBasis(MSH_TRI_1);
     case 1: return BasisFactory::getNodalBasis(MSH_TRI_3);
     case 2: return BasisFactory::getNodalBasis(MSH_TRI_6);
@@ -155,36 +136,15 @@ const nodalBasis* MTriangle::getFunctionSpace(int o) const
     case 9: return BasisFactory::getNodalBasis(MSH_TRI_55);
     case 10: return BasisFactory::getNodalBasis(MSH_TRI_66);
     default: Msg::Error("Order %d triangle function space not implemented", order);
-    }
   }
-  return 0;
+  return NULL;
 }
 
-const JacobianBasis* MTriangle::getJacobianFuncSpace(int o) const
+const JacobianBasis* MTriangle::getJacobianFuncSpace(int order) const
 {
+  if (order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
 
-  int order = (o == -1) ? getPolynomialOrder() : o;
-
-  int nf = getNumFaceVertices();
-
-  if ((nf == 0) && (o == -1)) {
-    switch (order) {
-    case 0: return BasisFactory::getJacobianBasis(MSH_TRI_1);
-    case 1: return BasisFactory::getJacobianBasis(MSH_TRI_3);
-    case 2: return BasisFactory::getJacobianBasis(MSH_TRI_6);
-    case 3: return BasisFactory::getJacobianBasis(MSH_TRI_9);
-    case 4: return BasisFactory::getJacobianBasis(MSH_TRI_12);
-    case 5: return BasisFactory::getJacobianBasis(MSH_TRI_15I);
-    case 6: return BasisFactory::getJacobianBasis(MSH_TRI_18);
-    case 7: return BasisFactory::getJacobianBasis(MSH_TRI_21I);
-    case 8: return BasisFactory::getJacobianBasis(MSH_TRI_24);
-    case 9: return BasisFactory::getJacobianBasis(MSH_TRI_27);
-    case 10: return BasisFactory::getJacobianBasis(MSH_TRI_30);
-    default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
-    }
-  }
-  else {
-    switch (order) {
+  switch (order) {
     case 1: return BasisFactory::getJacobianBasis(MSH_TRI_3);
     case 2: return BasisFactory::getJacobianBasis(MSH_TRI_6);
     case 3: return BasisFactory::getJacobianBasis(MSH_TRI_10);
@@ -196,9 +156,8 @@ const JacobianBasis* MTriangle::getJacobianFuncSpace(int o) const
     case 9: return BasisFactory::getJacobianBasis(MSH_TRI_55);
     case 10: return BasisFactory::getJacobianBasis(MSH_TRI_66);
     default: Msg::Error("Order %d triangle function space not implemented", order);
-    }
   }
-  return 0;
+  return NULL;
 }
 
 int MTriangleN::getNumEdgesRep(){ return 3 * CTX::instance()->mesh.numSubEdges; }