diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 52d0a067f5c9d9a0a415d05ab1f84a1a35b56624..29cfc61496c87f9cf79c21d5ac148d0cfebe1da6 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -198,31 +198,31 @@ const nodalBasis* MHexahedron::getFunctionSpace(int o) const
 
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_HEX_1);
-    case 1: return BasisFactory::create(MSH_HEX_8);
-    case 2: return BasisFactory::create(MSH_HEX_20);
-    case 3: return BasisFactory::create(MSH_HEX_56);
-    case 4: return BasisFactory::create(MSH_HEX_98);
-    case 5: return BasisFactory::create(MSH_HEX_152);
-    case 6: return BasisFactory::create(MSH_HEX_218);
-    case 7: return BasisFactory::create(MSH_HEX_296);
-    case 8: return BasisFactory::create(MSH_HEX_386);
-    case 9: return BasisFactory::create(MSH_HEX_488);
+    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_56);
+    case 4: return BasisFactory::getNodalBasis(MSH_HEX_98);
+    case 5: return BasisFactory::getNodalBasis(MSH_HEX_152);
+    case 6: return BasisFactory::getNodalBasis(MSH_HEX_218);
+    case 7: return BasisFactory::getNodalBasis(MSH_HEX_296);
+    case 8: return BasisFactory::getNodalBasis(MSH_HEX_386);
+    case 9: return BasisFactory::getNodalBasis(MSH_HEX_488);
     default: Msg::Error("Order %d hex function space not implemented", order); break;
     }
   }
   else {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_HEX_1);
-    case 1: return BasisFactory::create(MSH_HEX_8);
-    case 2: return BasisFactory::create(MSH_HEX_27);
-    case 3: return BasisFactory::create(MSH_HEX_64);
-    case 4: return BasisFactory::create(MSH_HEX_125);
-    case 5: return BasisFactory::create(MSH_HEX_216);
-    case 6: return BasisFactory::create(MSH_HEX_343);
-    case 7: return BasisFactory::create(MSH_HEX_512);
-    case 8: return BasisFactory::create(MSH_HEX_729);
-    case 9: return BasisFactory::create(MSH_HEX_1000);
+    case 0: return BasisFactory::getNodalBasis(MSH_HEX_1);
+    case 1: return BasisFactory::getNodalBasis(MSH_HEX_8);
+    case 2: return BasisFactory::getNodalBasis(MSH_HEX_27);
+    case 3: return BasisFactory::getNodalBasis(MSH_HEX_64);
+    case 4: return BasisFactory::getNodalBasis(MSH_HEX_125);
+    case 5: return BasisFactory::getNodalBasis(MSH_HEX_216);
+    case 6: return BasisFactory::getNodalBasis(MSH_HEX_343);
+    case 7: return BasisFactory::getNodalBasis(MSH_HEX_512);
+    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;
     }
   }
@@ -237,31 +237,31 @@ const JacobianBasis* MHexahedron::getJacobianFuncSpace(int o) const
 
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 0: return JacobianBasis::find(MSH_HEX_1);
-    case 1: return JacobianBasis::find(MSH_HEX_8);
-    case 2: return JacobianBasis::find(MSH_HEX_20);
-    case 3: return JacobianBasis::find(MSH_HEX_56);
-    case 4: return JacobianBasis::find(MSH_HEX_98);
-    case 5: return JacobianBasis::find(MSH_HEX_152);
-    case 6: return JacobianBasis::find(MSH_HEX_218);
-    case 7: return JacobianBasis::find(MSH_HEX_296);
-    case 8: return JacobianBasis::find(MSH_HEX_386);
-    case 9: return JacobianBasis::find(MSH_HEX_488);
+    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_56);
+    case 4: return BasisFactory::getJacobianBasis(MSH_HEX_98);
+    case 5: return BasisFactory::getJacobianBasis(MSH_HEX_152);
+    case 6: return BasisFactory::getJacobianBasis(MSH_HEX_218);
+    case 7: return BasisFactory::getJacobianBasis(MSH_HEX_296);
+    case 8: return BasisFactory::getJacobianBasis(MSH_HEX_386);
+    case 9: return BasisFactory::getJacobianBasis(MSH_HEX_488);
     default: Msg::Error("Order %d hex incomplete Jacobian function space not implemented", order); break;
     }
   }
   else {
     switch (order) {
-    case 0: return JacobianBasis::find(MSH_HEX_1);
-    case 1: return JacobianBasis::find(MSH_HEX_8);
-    case 2: return JacobianBasis::find(MSH_HEX_27);
-    case 3: return JacobianBasis::find(MSH_HEX_64);
-    case 4: return JacobianBasis::find(MSH_HEX_125);
-    case 5: return JacobianBasis::find(MSH_HEX_216);
-    case 6: return JacobianBasis::find(MSH_HEX_343);
-    case 7: return JacobianBasis::find(MSH_HEX_512);
-    case 8: return JacobianBasis::find(MSH_HEX_729);
-    case 9: return JacobianBasis::find(MSH_HEX_1000);
+    case 0: return BasisFactory::getJacobianBasis(MSH_HEX_1);
+    case 1: return BasisFactory::getJacobianBasis(MSH_HEX_8);
+    case 2: return BasisFactory::getJacobianBasis(MSH_HEX_27);
+    case 3: return BasisFactory::getJacobianBasis(MSH_HEX_64);
+    case 4: return BasisFactory::getJacobianBasis(MSH_HEX_125);
+    case 5: return BasisFactory::getJacobianBasis(MSH_HEX_216);
+    case 6: return BasisFactory::getJacobianBasis(MSH_HEX_343);
+    case 7: return BasisFactory::getJacobianBasis(MSH_HEX_512);
+    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;
     }
   }
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 91d940459530d5e6a6f88d1a82a95bf3d28601d0..7b33b69a68010be237e2cfd8e11513d6a1a8e711 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -16,17 +16,17 @@ const nodalBasis* MLine::getFunctionSpace(int o) const
   int order = (o == -1) ? getPolynomialOrder() : o;
   
   switch (order) {
-  case 0: return BasisFactory::create(MSH_LIN_1);
-  case 1: return BasisFactory::create(MSH_LIN_2);
-  case 2: return BasisFactory::create(MSH_LIN_3);
-  case 3: return BasisFactory::create(MSH_LIN_4);
-  case 4: return BasisFactory::create(MSH_LIN_5);
-  case 5: return BasisFactory::create(MSH_LIN_6);
-  case 6: return BasisFactory::create(MSH_LIN_7);
-  case 7: return BasisFactory::create(MSH_LIN_8);
-  case 8: return BasisFactory::create(MSH_LIN_9);
-  case 9: return BasisFactory::create(MSH_LIN_10);
-  case 10: return BasisFactory::create(MSH_LIN_11);
+  case 0: return BasisFactory::getNodalBasis(MSH_LIN_1);
+  case 1: return BasisFactory::getNodalBasis(MSH_LIN_2);
+  case 2: return BasisFactory::getNodalBasis(MSH_LIN_3);
+  case 3: return BasisFactory::getNodalBasis(MSH_LIN_4);
+  case 4: return BasisFactory::getNodalBasis(MSH_LIN_5);
+  case 5: return BasisFactory::getNodalBasis(MSH_LIN_6);
+  case 6: return BasisFactory::getNodalBasis(MSH_LIN_7);
+  case 7: return BasisFactory::getNodalBasis(MSH_LIN_8);
+  case 8: return BasisFactory::getNodalBasis(MSH_LIN_9);
+  case 9: return BasisFactory::getNodalBasis(MSH_LIN_10);
+  case 10: return BasisFactory::getNodalBasis(MSH_LIN_11);
   default: Msg::Error("Order %d line function space not implemented", order);
   }
   return 0;
@@ -37,16 +37,16 @@ const JacobianBasis* MLine::getJacobianFuncSpace(int o) const
   int order = (o == -1) ? getPolynomialOrder() : o;
   
   switch (order) {
-  case 1: return JacobianBasis::find(MSH_LIN_2);
-  case 2: return JacobianBasis::find(MSH_LIN_3);
-  case 3: return JacobianBasis::find(MSH_LIN_4);
-  case 4: return JacobianBasis::find(MSH_LIN_5);
-  case 5: return JacobianBasis::find(MSH_LIN_6);
-  case 6: return JacobianBasis::find(MSH_LIN_7);
-  case 7: return JacobianBasis::find(MSH_LIN_8);
-  case 8: return JacobianBasis::find(MSH_LIN_9);
-  case 9: return JacobianBasis::find(MSH_LIN_10);
-  case 10: return JacobianBasis::find(MSH_LIN_11);
+  case 1: return BasisFactory::getJacobianBasis(MSH_LIN_2);
+  case 2: return BasisFactory::getJacobianBasis(MSH_LIN_3);
+  case 3: return BasisFactory::getJacobianBasis(MSH_LIN_4);
+  case 4: return BasisFactory::getJacobianBasis(MSH_LIN_5);
+  case 5: return BasisFactory::getJacobianBasis(MSH_LIN_6);
+  case 6: return BasisFactory::getJacobianBasis(MSH_LIN_7);
+  case 7: return BasisFactory::getJacobianBasis(MSH_LIN_8);
+  case 8: return BasisFactory::getJacobianBasis(MSH_LIN_9);
+  case 9: return BasisFactory::getJacobianBasis(MSH_LIN_10);
+  case 10: return BasisFactory::getJacobianBasis(MSH_LIN_11);
   default: Msg::Error("Order %d line function space not implemented", order);
   }
   return 0;
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 66ae8582c12e137e95b91cb6c11905dc9f277d8e..08ccc7fbc11f2a4f2616dc505821e29cf0f93270 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -63,11 +63,11 @@ class MPoint : public MElement {
   }
   virtual const nodalBasis* getFunctionSpace(int o) const
   {
-    return BasisFactory::create(MSH_PNT);
+    return BasisFactory::getNodalBasis(MSH_PNT);
   }
   virtual const JacobianBasis* getJacobianFuncSpace(int o) const
   {
-    return JacobianBasis::find(MSH_PNT);
+    return BasisFactory::getJacobianBasis(MSH_PNT);
   }
   virtual bool isInside(double u, double v, double w) const
   {
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index f364c55674a05e1c43b27b5a237485a80b72c53c..508bf117d7641ec1d98cd441b6b66914ca59045e 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -40,17 +40,17 @@ const nodalBasis* MPrism::getFunctionSpace(int o) const
 
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_PRI_1);
-    case 1: return BasisFactory::create(MSH_PRI_6);
-    case 2: return BasisFactory::create(MSH_PRI_18);
+    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);
     }
   }
   else {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_PRI_1);
-    case 1: return BasisFactory::create(MSH_PRI_6);
-    case 2: return BasisFactory::create(MSH_PRI_18);
+    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);
     }
   }
@@ -65,15 +65,15 @@ const JacobianBasis* MPrism::getJacobianFuncSpace(int o) const
   
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 1: return JacobianBasis::find(MSH_PRI_6);
-    case 2: return JacobianBasis::find(MSH_PRI_18);
+    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);
     }
   }
   else { 
     switch (order) {
-    case 1: return JacobianBasis::find(MSH_PRI_6);
-    case 2: return JacobianBasis::find(MSH_PRI_18);
+    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);
     }
   }
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index f6e3c41307385d07a0dcc5f37443598f632e70cb..1a44f7c4b67845a3e20d77ed584827b0a31a9d10 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -35,31 +35,31 @@ const nodalBasis* MPyramid::getFunctionSpace(int o) const
 
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 1: return BasisFactory::create(MSH_PYR_5);
-    case 2: return BasisFactory::create(MSH_PYR_14);
-    case 0: return BasisFactory::create(MSH_PYR_1);
-    case 3: return BasisFactory::create(MSH_PYR_29);
-    case 4: return BasisFactory::create(MSH_PYR_50);
-    case 5: return BasisFactory::create(MSH_PYR_77);
-    case 6: return BasisFactory::create(MSH_PYR_110);
-    case 7: return BasisFactory::create(MSH_PYR_149);
-    case 8: return BasisFactory::create(MSH_PYR_194);
-    case 9: return BasisFactory::create(MSH_PYR_245);
+    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_29);
+    case 4: return BasisFactory::getNodalBasis(MSH_PYR_50);
+    case 5: return BasisFactory::getNodalBasis(MSH_PYR_77);
+    case 6: return BasisFactory::getNodalBasis(MSH_PYR_110);
+    case 7: return BasisFactory::getNodalBasis(MSH_PYR_149);
+    case 8: return BasisFactory::getNodalBasis(MSH_PYR_194);
+    case 9: return BasisFactory::getNodalBasis(MSH_PYR_245);
     default: Msg::Error("Order %d pyramid function space not implemented", order);
     }
   }
   else {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_PYR_1);
-    case 1: return BasisFactory::create(MSH_PYR_5);
-    case 2: return BasisFactory::create(MSH_PYR_14);
-    case 3: return BasisFactory::create(MSH_PYR_30);
-    case 4: return BasisFactory::create(MSH_PYR_55);
-    case 5: return BasisFactory::create(MSH_PYR_91);
-    case 6: return BasisFactory::create(MSH_PYR_140);
-    case 7: return BasisFactory::create(MSH_PYR_204);
-    case 8: return BasisFactory::create(MSH_PYR_285);
-    case 9: return BasisFactory::create(MSH_PYR_385);
+    case 0: return BasisFactory::getNodalBasis(MSH_PYR_1);
+    case 1: return BasisFactory::getNodalBasis(MSH_PYR_5);
+    case 2: return BasisFactory::getNodalBasis(MSH_PYR_14);
+    case 3: return BasisFactory::getNodalBasis(MSH_PYR_30);
+    case 4: return BasisFactory::getNodalBasis(MSH_PYR_55);
+    case 5: return BasisFactory::getNodalBasis(MSH_PYR_91);
+    case 6: return BasisFactory::getNodalBasis(MSH_PYR_140);
+    case 7: return BasisFactory::getNodalBasis(MSH_PYR_204);
+    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);
     }
   }
@@ -71,9 +71,9 @@ const JacobianBasis* MPyramid::getJacobianFuncSpace(int o) const
   int order = (o == -1) ? getPolynomialOrder() : o;
 
   switch (order) {
-    case 1: return JacobianBasis::find(MSH_PYR_5);
-    case 2: return JacobianBasis::find(MSH_PYR_14);
-    case 3: return JacobianBasis::find(MSH_PYR_30);
+    case 1: return BasisFactory::getJacobianBasis(MSH_PYR_5);
+    case 2: return BasisFactory::getJacobianBasis(MSH_PYR_14);
+    case 3: return BasisFactory::getJacobianBasis(MSH_PYR_30);
     default: Msg::Error("Order %d pyramid function space not implemented", order); break;
   }
 
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 6808e69e23b9fb7b35df7cad93deb902c9ae4e7d..0f6b54877c77af7848f6e4f49cf5e1b9e9f503c5 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -25,31 +25,31 @@ const nodalBasis* MQuadrangle::getFunctionSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-      case 0: return BasisFactory::create(MSH_QUA_1);
-      case 1: return BasisFactory::create(MSH_QUA_4);
-      case 2: return BasisFactory::create(MSH_QUA_8);
-      case 3: return BasisFactory::create(MSH_QUA_12);
-      case 4: return BasisFactory::create(MSH_QUA_16I);
-      case 5: return BasisFactory::create(MSH_QUA_20);
-      case 6: return BasisFactory::create(MSH_QUA_24);
-      case 7: return BasisFactory::create(MSH_QUA_28);
-      case 8: return BasisFactory::create(MSH_QUA_32);
-      case 9: return BasisFactory::create(MSH_QUA_36I);
-      case 10: return BasisFactory::create(MSH_QUA_40);
+      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);
     }
   }
   switch (order) {
-    case 0: return BasisFactory::create(MSH_QUA_1);
-    case 1: return BasisFactory::create(MSH_QUA_4);
-    case 2: return BasisFactory::create(MSH_QUA_9);
-    case 3: return BasisFactory::create(MSH_QUA_16);
-    case 4: return BasisFactory::create(MSH_QUA_25);
-    case 5: return BasisFactory::create(MSH_QUA_36);
-    case 6: return BasisFactory::create(MSH_QUA_49);
-    case 7: return BasisFactory::create(MSH_QUA_64);
-    case 8: return BasisFactory::create(MSH_QUA_81);
-    case 9: return BasisFactory::create(MSH_QUA_100);
-    case 10: return BasisFactory::create(MSH_QUA_121);
+    case 0: return BasisFactory::getNodalBasis(MSH_QUA_1);
+    case 1: return BasisFactory::getNodalBasis(MSH_QUA_4);
+    case 2: return BasisFactory::getNodalBasis(MSH_QUA_9);
+    case 3: return BasisFactory::getNodalBasis(MSH_QUA_16);
+    case 4: return BasisFactory::getNodalBasis(MSH_QUA_25);
+    case 5: return BasisFactory::getNodalBasis(MSH_QUA_36);
+    case 6: return BasisFactory::getNodalBasis(MSH_QUA_49);
+    case 7: return BasisFactory::getNodalBasis(MSH_QUA_64);
+    case 8: return BasisFactory::getNodalBasis(MSH_QUA_81);
+    case 9: return BasisFactory::getNodalBasis(MSH_QUA_100);
+    case 10: return BasisFactory::getNodalBasis(MSH_QUA_121);
     default: Msg::Error("Order %d quadrangle function space not implemented", order);
   }
   return 0;
@@ -63,29 +63,29 @@ const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-      case 1: return JacobianBasis::find(MSH_QUA_4);
-      case 2: return JacobianBasis::find(MSH_QUA_8);
-      case 3: return JacobianBasis::find(MSH_QUA_12);
-      case 4: return JacobianBasis::find(MSH_QUA_16I);
-      case 5: return JacobianBasis::find(MSH_QUA_20);
-      case 6: return JacobianBasis::find(MSH_QUA_24);
-      case 7: return JacobianBasis::find(MSH_QUA_28);
-      case 8: return JacobianBasis::find(MSH_QUA_32);
-      case 9: return JacobianBasis::find(MSH_QUA_36I);
-      case 10: return JacobianBasis::find(MSH_QUA_40);
+      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);
     }
   }
   switch (order) {
-    case 1: return JacobianBasis::find(MSH_QUA_4);
-    case 2: return JacobianBasis::find(MSH_QUA_9);
-    case 3: return JacobianBasis::find(MSH_QUA_16);
-    case 4: return JacobianBasis::find(MSH_QUA_25);
-    case 5: return JacobianBasis::find(MSH_QUA_36);
-    case 6: return JacobianBasis::find(MSH_QUA_49);
-    case 7: return JacobianBasis::find(MSH_QUA_64);
-    case 8: return JacobianBasis::find(MSH_QUA_81);
-    case 9: return JacobianBasis::find(MSH_QUA_100);
-    case 10: return JacobianBasis::find(MSH_QUA_121);
+    case 1: return BasisFactory::getJacobianBasis(MSH_QUA_4);
+    case 2: return BasisFactory::getJacobianBasis(MSH_QUA_9);
+    case 3: return BasisFactory::getJacobianBasis(MSH_QUA_16);
+    case 4: return BasisFactory::getJacobianBasis(MSH_QUA_25);
+    case 5: return BasisFactory::getJacobianBasis(MSH_QUA_36);
+    case 6: return BasisFactory::getJacobianBasis(MSH_QUA_49);
+    case 7: return BasisFactory::getJacobianBasis(MSH_QUA_64);
+    case 8: return BasisFactory::getJacobianBasis(MSH_QUA_81);
+    case 9: return BasisFactory::getJacobianBasis(MSH_QUA_100);
+    case 10: return BasisFactory::getJacobianBasis(MSH_QUA_121);
     default: Msg::Error("Order %d quadrangle function space not implemented", order);
   }
   return 0;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 5d1d29a0b1f168d18d9ac9161da197ce502d15f3..af351888636a1c2b6c943304b476c8e647994dc5 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -8,7 +8,7 @@
 #include "Numeric.h"
 #include "Context.h"
 #include "BasisFactory.h"
-#include "JacobianBasis.h"
+
 
 #if defined(HAVE_MESH)
 #include "qualityMeasures.h"
@@ -108,33 +108,33 @@ const nodalBasis* MTetrahedron::getFunctionSpace(int o) const
 
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_TET_1);
-    case 1: return BasisFactory::create(MSH_TET_4);
-    case 2: return BasisFactory::create(MSH_TET_10);
-    case 3: return BasisFactory::create(MSH_TET_20);
-    case 4: return BasisFactory::create(MSH_TET_34);
-    case 5: return BasisFactory::create(MSH_TET_52);
-    case 6: return BasisFactory::create(MSH_TET_74);
-    case 7: return BasisFactory::create(MSH_TET_100);
-    case 8: return BasisFactory::create(MSH_TET_130);
-    case 9: return BasisFactory::create(MSH_TET_164);
-    case 10: return BasisFactory::create(MSH_TET_202);
+    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_20);
+    case 4: return BasisFactory::getNodalBasis(MSH_TET_34);
+    case 5: return BasisFactory::getNodalBasis(MSH_TET_52);
+    case 6: return BasisFactory::getNodalBasis(MSH_TET_74);
+    case 7: return BasisFactory::getNodalBasis(MSH_TET_100);
+    case 8: return BasisFactory::getNodalBasis(MSH_TET_130);
+    case 9: return BasisFactory::getNodalBasis(MSH_TET_164);
+    case 10: return BasisFactory::getNodalBasis(MSH_TET_202);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
     }
   }
   else {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_TET_1);
-    case 1: return BasisFactory::create(MSH_TET_4);
-    case 2: return BasisFactory::create(MSH_TET_10);
-    case 3: return BasisFactory::create(MSH_TET_20);
-    case 4: return BasisFactory::create(MSH_TET_35);
-    case 5: return BasisFactory::create(MSH_TET_56);
-    case 6: return BasisFactory::create(MSH_TET_84);
-    case 7: return BasisFactory::create(MSH_TET_120);
-    case 8: return BasisFactory::create(MSH_TET_165);
-    case 9: return BasisFactory::create(MSH_TET_220);
-    case 10: return BasisFactory::create(MSH_TET_286);
+    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_20);
+    case 4: return BasisFactory::getNodalBasis(MSH_TET_35);
+    case 5: return BasisFactory::getNodalBasis(MSH_TET_56);
+    case 6: return BasisFactory::getNodalBasis(MSH_TET_84);
+    case 7: return BasisFactory::getNodalBasis(MSH_TET_120);
+    case 8: return BasisFactory::getNodalBasis(MSH_TET_165);
+    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);
     }
   }
@@ -149,26 +149,26 @@ const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int o) const
 
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 1: return JacobianBasis::find(MSH_TET_4);
-    case 2: return JacobianBasis::find(MSH_TET_10);
-    case 3: return JacobianBasis::find(MSH_TET_20);
-    case 4: return JacobianBasis::find(MSH_TET_34);
-    case 5: return JacobianBasis::find(MSH_TET_52);
+    case 1: return BasisFactory::getJacobianBasis(MSH_TET_4);
+    case 2: return BasisFactory::getJacobianBasis(MSH_TET_10);
+    case 3: return BasisFactory::getJacobianBasis(MSH_TET_20);
+    case 4: return BasisFactory::getJacobianBasis(MSH_TET_34);
+    case 5: return BasisFactory::getJacobianBasis(MSH_TET_52);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
     }
   }
   else {
     switch (order) {
-    case 1: return JacobianBasis::find(MSH_TET_4);
-    case 2: return JacobianBasis::find(MSH_TET_10);
-    case 3: return JacobianBasis::find(MSH_TET_20);
-    case 4: return JacobianBasis::find(MSH_TET_35);
-    case 5: return JacobianBasis::find(MSH_TET_56);
-    case 6: return JacobianBasis::find(MSH_TET_84);
-    case 7: return JacobianBasis::find(MSH_TET_120);
-    case 8: return JacobianBasis::find(MSH_TET_165);
-    case 9: return JacobianBasis::find(MSH_TET_220);
-    case 10: return JacobianBasis::find(MSH_TET_286);
+    case 1: return BasisFactory::getJacobianBasis(MSH_TET_4);
+    case 2: return BasisFactory::getJacobianBasis(MSH_TET_10);
+    case 3: return BasisFactory::getJacobianBasis(MSH_TET_20);
+    case 4: return BasisFactory::getJacobianBasis(MSH_TET_35);
+    case 5: return BasisFactory::getJacobianBasis(MSH_TET_56);
+    case 6: return BasisFactory::getJacobianBasis(MSH_TET_84);
+    case 7: return BasisFactory::getJacobianBasis(MSH_TET_120);
+    case 8: return BasisFactory::getJacobianBasis(MSH_TET_165);
+    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);
     }
   }
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 49d68dee7b8735233c4a220eb8a2fc81c92765a7..3c73eb8de45409d0096bfc6f37af017bca679e52 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -127,33 +127,33 @@ const nodalBasis* MTriangle::getFunctionSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-    case 0: return BasisFactory::create(MSH_TRI_1);
-    case 1: return BasisFactory::create(MSH_TRI_3);
-    case 2: return BasisFactory::create(MSH_TRI_6);
-    case 3: return BasisFactory::create(MSH_TRI_9);
-    case 4: return BasisFactory::create(MSH_TRI_12);
-    case 5: return BasisFactory::create(MSH_TRI_15I);
-    case 6: return BasisFactory::create(MSH_TRI_18);
-    case 7: return BasisFactory::create(MSH_TRI_21I);
-    case 8: return BasisFactory::create(MSH_TRI_24);
-    case 9: return BasisFactory::create(MSH_TRI_27);
-    case 10: return BasisFactory::create(MSH_TRI_30);
+    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) {
-    case 0: return BasisFactory::create(MSH_TRI_1);
-    case 1: return BasisFactory::create(MSH_TRI_3);
-    case 2: return BasisFactory::create(MSH_TRI_6);
-    case 3: return BasisFactory::create(MSH_TRI_10);
-    case 4: return BasisFactory::create(MSH_TRI_15);
-    case 5: return BasisFactory::create(MSH_TRI_21);
-    case 6: return BasisFactory::create(MSH_TRI_28);
-    case 7: return BasisFactory::create(MSH_TRI_36);
-    case 8: return BasisFactory::create(MSH_TRI_45);
-    case 9: return BasisFactory::create(MSH_TRI_55);
-    case 10: return BasisFactory::create(MSH_TRI_66);
+    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_10);
+    case 4: return BasisFactory::getNodalBasis(MSH_TRI_15);
+    case 5: return BasisFactory::getNodalBasis(MSH_TRI_21);
+    case 6: return BasisFactory::getNodalBasis(MSH_TRI_28);
+    case 7: return BasisFactory::getNodalBasis(MSH_TRI_36);
+    case 8: return BasisFactory::getNodalBasis(MSH_TRI_45);
+    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);
     }
   }
@@ -169,32 +169,32 @@ const JacobianBasis* MTriangle::getJacobianFuncSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-    case 0: return JacobianBasis::find(MSH_TRI_1);
-    case 1: return JacobianBasis::find(MSH_TRI_3);
-    case 2: return JacobianBasis::find(MSH_TRI_6);
-    case 3: return JacobianBasis::find(MSH_TRI_9);
-    case 4: return JacobianBasis::find(MSH_TRI_12);
-    case 5: return JacobianBasis::find(MSH_TRI_15I);
-    case 6: return JacobianBasis::find(MSH_TRI_18);
-    case 7: return JacobianBasis::find(MSH_TRI_21I);
-    case 8: return JacobianBasis::find(MSH_TRI_24);
-    case 9: return JacobianBasis::find(MSH_TRI_27);
-    case 10: return JacobianBasis::find(MSH_TRI_30);
+    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) {
-    case 1: return JacobianBasis::find(MSH_TRI_3);
-    case 2: return JacobianBasis::find(MSH_TRI_6);
-    case 3: return JacobianBasis::find(MSH_TRI_10);
-    case 4: return JacobianBasis::find(MSH_TRI_15);
-    case 5: return JacobianBasis::find(MSH_TRI_21);
-    case 6: return JacobianBasis::find(MSH_TRI_28);
-    case 7: return JacobianBasis::find(MSH_TRI_36);
-    case 8: return JacobianBasis::find(MSH_TRI_45);
-    case 9: return JacobianBasis::find(MSH_TRI_55);
-    case 10: return JacobianBasis::find(MSH_TRI_66);
+    case 1: return BasisFactory::getJacobianBasis(MSH_TRI_3);
+    case 2: return BasisFactory::getJacobianBasis(MSH_TRI_6);
+    case 3: return BasisFactory::getJacobianBasis(MSH_TRI_10);
+    case 4: return BasisFactory::getJacobianBasis(MSH_TRI_15);
+    case 5: return BasisFactory::getJacobianBasis(MSH_TRI_21);
+    case 6: return BasisFactory::getJacobianBasis(MSH_TRI_28);
+    case 7: return BasisFactory::getJacobianBasis(MSH_TRI_36);
+    case 8: return BasisFactory::getJacobianBasis(MSH_TRI_45);
+    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);
     }
   }
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index c79bd36304258270cd7edd6daeaa81e5d1991d45..7ce3f119ff83d066e448635b42b60aec112894f5 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -501,7 +501,7 @@ double discreteEdge::curvature(double par) const
   double c0, c1;
   Curvature& curvature  = Curvature::getInstance();
   if( !Curvature::valueAlreadyComputed() ) {
-    std::cout << "Need to compute discrete curvature (in discreteEdge)" << std::endl;
+    Msg::Warning("Need to compute discrete curvature (in discreteEdge)");
     Curvature::typeOfCurvature type = Curvature::RUSIN; //RUSIN; //RBF
     curvature.computeCurvature(model(), type); 
   }
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index dbc93f24e1fd1069001adec09438a5ac8ddd1f62..59f54a2bb052d762b817af7be708a26cbae48e35 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -472,14 +472,14 @@ static int getNewFacePoints(int numPrimaryFacePoints, int nPts, fullMatrix<doubl
     switch (nPts){
     case 0 : break;
     case 1 : break;
-    case 2 : points = BasisFactory::create(MSH_TRI_10)->points; start = 9; break;
-    case 3 : points = BasisFactory::create(MSH_TRI_15)->points; start = 12; break;
-    case 4 : points = BasisFactory::create(MSH_TRI_21)->points; start = 15; break;
-    case 5 : points = BasisFactory::create(MSH_TRI_28)->points; start = 18; break;
-    case 6 : points = BasisFactory::create(MSH_TRI_36)->points; start = 21; break;
-    case 7 : points = BasisFactory::create(MSH_TRI_45)->points; start = 24; break;
-    case 8 : points = BasisFactory::create(MSH_TRI_55)->points; start = 27; break;
-    case 9 : points = BasisFactory::create(MSH_TRI_66)->points; start = 30; break;
+    case 2 : points = BasisFactory::getNodalBasis(MSH_TRI_10)->points; start = 9; break;
+    case 3 : points = BasisFactory::getNodalBasis(MSH_TRI_15)->points; start = 12; break;
+    case 4 : points = BasisFactory::getNodalBasis(MSH_TRI_21)->points; start = 15; break;
+    case 5 : points = BasisFactory::getNodalBasis(MSH_TRI_28)->points; start = 18; break;
+    case 6 : points = BasisFactory::getNodalBasis(MSH_TRI_36)->points; start = 21; break;
+    case 7 : points = BasisFactory::getNodalBasis(MSH_TRI_45)->points; start = 24; break;
+    case 8 : points = BasisFactory::getNodalBasis(MSH_TRI_55)->points; start = 27; break;
+    case 9 : points = BasisFactory::getNodalBasis(MSH_TRI_66)->points; start = 30; break;
     default :
       Msg::Error("getFaceVertices not implemented for order %i", nPts +1);
       break;
@@ -488,14 +488,14 @@ static int getNewFacePoints(int numPrimaryFacePoints, int nPts, fullMatrix<doubl
   case 4:
     switch (nPts){
     case 0 : break;
-    case 1 : points = BasisFactory::create(MSH_QUA_9)->points; break;
-    case 2 : points = BasisFactory::create(MSH_QUA_16)->points; break;
-    case 3 : points = BasisFactory::create(MSH_QUA_25)->points; break;
-    case 4 : points = BasisFactory::create(MSH_QUA_36)->points; break;
-    case 5 : points = BasisFactory::create(MSH_QUA_49)->points; break;
-    case 6 : points = BasisFactory::create(MSH_QUA_64)->points; break;
-    case 7 : points = BasisFactory::create(MSH_QUA_81)->points; break;
-    case 8 : points = BasisFactory::create(MSH_QUA_100)->points; break;
+    case 1 : points = BasisFactory::getNodalBasis(MSH_QUA_9)->points; break;
+    case 2 : points = BasisFactory::getNodalBasis(MSH_QUA_16)->points; break;
+    case 3 : points = BasisFactory::getNodalBasis(MSH_QUA_25)->points; break;
+    case 4 : points = BasisFactory::getNodalBasis(MSH_QUA_36)->points; break;
+    case 5 : points = BasisFactory::getNodalBasis(MSH_QUA_49)->points; break;
+    case 6 : points = BasisFactory::getNodalBasis(MSH_QUA_64)->points; break;
+    case 7 : points = BasisFactory::getNodalBasis(MSH_QUA_81)->points; break;
+    case 8 : points = BasisFactory::getNodalBasis(MSH_QUA_100)->points; break;
     default :
       Msg::Error("getFaceVertices not implemented for order %i", nPts +1);
       break;
@@ -599,15 +599,15 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
     case 0: return;
     case 1: return;
     case 2:
-      BasisFactory::create(MSH_TET_20)->points.print();
-      points = BasisFactory::create(MSH_TET_20)->points; break;
-    case 3: points = BasisFactory::create(MSH_TET_35)->points; break;
-    case 4: points = BasisFactory::create(MSH_TET_56)->points; break;
-    case 5: points = BasisFactory::create(MSH_TET_84)->points; break;
-    case 6: points = BasisFactory::create(MSH_TET_120)->points; break;
-    case 7: points = BasisFactory::create(MSH_TET_165)->points; break;
-    case 8: points = BasisFactory::create(MSH_TET_220)->points; break;
-    case 9: points = BasisFactory::create(MSH_TET_286)->points; break;
+      BasisFactory::getNodalBasis(MSH_TET_20)->points.print();
+      points = BasisFactory::getNodalBasis(MSH_TET_20)->points; break;
+    case 3: points = BasisFactory::getNodalBasis(MSH_TET_35)->points; break;
+    case 4: points = BasisFactory::getNodalBasis(MSH_TET_56)->points; break;
+    case 5: points = BasisFactory::getNodalBasis(MSH_TET_84)->points; break;
+    case 6: points = BasisFactory::getNodalBasis(MSH_TET_120)->points; break;
+    case 7: points = BasisFactory::getNodalBasis(MSH_TET_165)->points; break;
+    case 8: points = BasisFactory::getNodalBasis(MSH_TET_220)->points; break;
+    case 9: points = BasisFactory::getNodalBasis(MSH_TET_286)->points; break;
     default:
       Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
       break;
@@ -617,14 +617,14 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
   case TYPE_HEX :
     switch (nPts){
     case 0: return;
-    case 1: points = BasisFactory::create(MSH_HEX_27)->points; break;
-    case 2: points = BasisFactory::create(MSH_HEX_64)->points; break;
-    case 3: points = BasisFactory::create(MSH_HEX_125)->points; break;
-    case 4: points = BasisFactory::create(MSH_HEX_216)->points; break;
-    case 5: points = BasisFactory::create(MSH_HEX_343)->points; break;
-    case 6: points = BasisFactory::create(MSH_HEX_512)->points; break;
-    case 7: points = BasisFactory::create(MSH_HEX_729)->points; break;
-    case 8: points = BasisFactory::create(MSH_HEX_1000)->points; break;
+    case 1: points = BasisFactory::getNodalBasis(MSH_HEX_27)->points; break;
+    case 2: points = BasisFactory::getNodalBasis(MSH_HEX_64)->points; break;
+    case 3: points = BasisFactory::getNodalBasis(MSH_HEX_125)->points; break;
+    case 4: points = BasisFactory::getNodalBasis(MSH_HEX_216)->points; break;
+    case 5: points = BasisFactory::getNodalBasis(MSH_HEX_343)->points; break;
+    case 6: points = BasisFactory::getNodalBasis(MSH_HEX_512)->points; break;
+    case 7: points = BasisFactory::getNodalBasis(MSH_HEX_729)->points; break;
+    case 8: points = BasisFactory::getNodalBasis(MSH_HEX_1000)->points; break;
     default :
       Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
       break;
@@ -635,13 +635,13 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
     switch (nPts){
     case 0:
     case 1: return;
-    case 2: points = BasisFactory::create(MSH_PYR_30)->points; break;
-    case 3: points = BasisFactory::create(MSH_PYR_55)->points; break;
-    case 4: points = BasisFactory::create(MSH_PYR_91)->points; break;
-    case 5: points = BasisFactory::create(MSH_PYR_140)->points; break;
-    case 6: points = BasisFactory::create(MSH_PYR_204)->points; break;
-    case 7: points = BasisFactory::create(MSH_PYR_285)->points; break;
-    case 8: points = BasisFactory::create(MSH_PYR_385)->points; break;
+    case 2: points = BasisFactory::getNodalBasis(MSH_PYR_30)->points; break;
+    case 3: points = BasisFactory::getNodalBasis(MSH_PYR_55)->points; break;
+    case 4: points = BasisFactory::getNodalBasis(MSH_PYR_91)->points; break;
+    case 5: points = BasisFactory::getNodalBasis(MSH_PYR_140)->points; break;
+    case 6: points = BasisFactory::getNodalBasis(MSH_PYR_204)->points; break;
+    case 7: points = BasisFactory::getNodalBasis(MSH_PYR_285)->points; break;
+    case 8: points = BasisFactory::getNodalBasis(MSH_PYR_385)->points; break;
     default :
       Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
       break;
diff --git a/Mesh/meshGFaceRecombine.h b/Mesh/meshGFaceRecombine.h
index 486f70f3661b5251cd4527820910ddef73748d1f..94c7cc8f61df076f2f3eeec27c8dec6d483d88b6 100644
--- a/Mesh/meshGFaceRecombine.h
+++ b/Mesh/meshGFaceRecombine.h
@@ -416,7 +416,7 @@ class Rec2DChange {
     void *_info;
   
   public :
-    Rec2DChange() {Msg::Error("[Rec2DChange] I should not be created in this manner");}
+    Rec2DChange() {Msg::Fatal("[Rec2DChange] I should not be created in this manner");}
     Rec2DChange(int);
     Rec2DChange(Rec2DEdge*, bool toHide = false);
     Rec2DChange(Rec2DVertex*, bool toHide = false);
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index 46adb2e19f96d142d09490838f34a672a827956b..36730f595046f308653cfaa78cc2de83e15f9a4c 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -11,9 +11,11 @@
 #include "BasisFactory.h"
 
 std::map<int, nodalBasis*> BasisFactory::fs;
+std::map<int, JacobianBasis*> BasisFactory::js;
+std::map<int, bezierBasis*> BasisFactory::bs;
 
-const nodalBasis* BasisFactory::create(int elementType) {
-
+const nodalBasis* BasisFactory::getNodalBasis(int elementType)
+{
   // If the Basis has already been built, return it.
   std::map<int, nodalBasis*>::const_iterator it = fs.find(elementType);
   if (it != fs.end()) {
@@ -22,7 +24,7 @@ const nodalBasis* BasisFactory::create(int elementType) {
   // Get the parent type to see which kind of basis
   // we want to create
   int parentType = MElement::ParentTypeFromTag(elementType);
-  nodalBasis* B = 0;
+  nodalBasis* F = NULL;
 
   switch(parentType) {
     case(TYPE_PNT):
@@ -32,19 +34,40 @@ const nodalBasis* BasisFactory::create(int elementType) {
     case(TYPE_PRI):
     case(TYPE_TET):
     case(TYPE_HEX):
-      B = new polynomialBasis(elementType);
+      F = new polynomialBasis(elementType);
       break;
     case(TYPE_PYR):
-      B = new pyramidalBasis(elementType);
+      F = new pyramidalBasis(elementType);
       break;
     default:
       Msg::Error("Unknown type of element.");
-      return 0;
+      return NULL;
   }
 
   // FIXME: check if already exists to deallocate if necessary
-  fs.insert(std::make_pair(elementType, B));
+  fs.insert(std::make_pair(elementType, F));
 
   return fs[elementType];
+}
+
+const bezierBasis* BasisFactory::getBezierBasis(int elementType)
+{
+  std::map<int, bezierBasis*>::const_iterator it = bs.find(elementType);
+  if (it != bs.end())
+    return it->second;
+
+  bezierBasis* B = new bezierBasis(elementType);
+  if (B) bs.insert(std::make_pair(elementType, B));
+  return B;
+}
+
+const JacobianBasis* BasisFactory::getJacobianBasis(int elementType)
+{
+  std::map<int, JacobianBasis*>::const_iterator it = js.find(elementType);
+  if (it != js.end())
+    return it->second;
 
+  JacobianBasis* J = new JacobianBasis(elementType);
+  if (J) js.insert(std::make_pair(elementType, J));
+  return J;
 }
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 7932c3eae36b1fbfdc0017ae8b0d3b1ce441808f..af4f6d7e980a72033f912bb05d10f5d0741bbe5b 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -9,16 +9,21 @@
 #include "MElement.h"
 #include "MPyramid.h"
 #include "nodalBasis.h"
+#include "JacobianBasis.h"
+
 
 class BasisFactory
 {
   private:
     static std::map<int, nodalBasis*> fs;
+    static std::map<int, bezierBasis*> bs;
+    static std::map<int, JacobianBasis*> js;
 
   public:
-    // Caution: the returned pointer should be
-    // checked against 0.
-    static const nodalBasis* create(int elementType);
+    // Caution: the returned pointer can be NULL
+    static const nodalBasis* getNodalBasis(int elementType);
+    static const bezierBasis* getBezierBasis(int elementType);
+    static const JacobianBasis* getJacobianBasis(int elementType);
 };
 
 #endif
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 0931e656452fa2b07ab2e42cf05a6cedc55ab4b4..231bfe602cdffc0a6525889903140f41d5362a9e 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -8,12 +8,13 @@ set(SRC
     fullMatrix.cpp
   BasisFactory.cpp
     nodalBasis.cpp
-    polynomialBasis.cpp
+	polynomialBasis.cpp
+	pyramidalBasis.cpp
+	BergotBasis.cpp
+	jacobiPolynomials.cpp
+	legendrePolynomials.cpp
+    bezierBasis.cpp
     JacobianBasis.cpp
-    BergotBasis.cpp
-    pyramidalBasis.cpp
-    jacobiPolynomials.cpp
-    legendrePolynomials.cpp
   pointsGenerators.cpp
   GaussIntegration.cpp
     GaussQuadratureLin.cpp
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 8ac87aa574172efe02dab879fe0095988e14f68e..4b4e6610327f72b9d8cfe484bbe4cec8b8f12625 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -5,7 +5,7 @@
 
 #include "GmshDefines.h"
 #include "GmshMessage.h"
-#include "JacobianBasis.h"
+
 #include <vector>
 #include "polynomialBasis.h"
 #include "pyramidalBasis.h"
@@ -13,1068 +13,6 @@
 #include "BasisFactory.h"
 #include "Numeric.h"
 
-// Bezier Exponents
-static fullMatrix<double> generate1DExponents(int order)
-{
-  fullMatrix<double> exponents(order + 1, 1);
-  exponents(0, 0) = 0;
-  if (order > 0) {
-    exponents(1, 0) = order;
-    for (int i = 2; i < order + 1; i++) exponents(i, 0) = i - 1;
-  }
-
-  return exponents;
-}
-
-static fullMatrix<double> generateExponentsTriangle(int order)
-{
-  int nbPoints = (order + 1) * (order + 2) / 2;
-  fullMatrix<double> exponents(nbPoints, 2);
-
-  exponents(0, 0) = 0;
-  exponents(0, 1) = 0;
-
-  if (order > 0) {
-    exponents(1, 0) = order;
-    exponents(1, 1) = 0;
-    exponents(2, 0) = 0;
-    exponents(2, 1) = order;
-
-    if (order > 1) {
-      int index = 3, ksi = 0, eta = 0;
-
-      for (int i = 0; i < order - 1; i++, index++) {
-        ksi++;
-        exponents(index, 0) = ksi;
-        exponents(index, 1) = eta;
-      }
-
-      ksi = order;
-
-      for (int i = 0; i < order - 1; i++, index++) {
-        ksi--;
-        eta++;
-        exponents(index, 0) = ksi;
-        exponents(index, 1) = eta;
-      }
-
-      eta = order;
-      ksi = 0;
-
-      for (int i = 0; i < order - 1; i++, index++) {
-        eta--;
-        exponents(index, 0) = ksi;
-        exponents(index, 1) = eta;
-      }
-
-      if (order > 2) {
-        fullMatrix<double> inner = generateExponentsTriangle(order - 3);
-        inner.add(1);
-        exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
-      }
-    }
-  }
-
-  return exponents;
-}
-static fullMatrix<double> generateExponentsQuad(int order)
-{
-  int nbPoints = (order+1)*(order+1);
-  fullMatrix<double> exponents(nbPoints, 2);
-
-  exponents(0, 0) = 0;
-  exponents(0, 1) = 0;
-
-  if (order > 0) {
-    exponents(1, 0) = order;
-    exponents(1, 1) = 0;
-    exponents(2, 0) = order;
-    exponents(2, 1) = order;
-    exponents(3, 0) = 0;
-    exponents(3, 1) = order;
-
-    if (order > 1) {
-      int index = 4;
-      const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
-      for (int iedge=0; iedge<4; iedge++) {
-        int p0 = edges[iedge][0];
-        int p1 = edges[iedge][1];
-        for (int i = 1; i < order; i++, index++) {
-          exponents(index, 0) = exponents(p0, 0) + i*(exponents(p1,0)-exponents(p0,0))/order;
-          exponents(index, 1) = exponents(p0, 1) + i*(exponents(p1,1)-exponents(p0,1))/order;
-        }
-      }
-      if (order > 2) {
-        fullMatrix<double> inner = generateExponentsQuad(order - 2);
-        inner.add(1);
-        exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
-      }
-    }
-  }
-
-  //  exponents.print("expq");
-
-  return exponents;
-}
-static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-
-static void nodepositionface0(int order, double *u, double *v, double *w)
-{ // uv surface - orientation v0-v2-v1
-  int ndofT = nbdoftriangle(order);
-  if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
-  u[0]= 0.;    v[0]= 0.;    w[0] = 0.;
-  u[1]= 0.;    v[1]= order; w[1] = 0.;
-  u[2]= order; v[2]= 0.;    w[2] = 0.;
-
-  // edges
-  for (int k = 0; k < (order - 1); k++){
-    u[3 + k] = 0.;
-    v[3 + k] = k + 1;
-    w[3 + k] = 0.;
-
-    u[3 + order - 1 + k] = k + 1;
-    v[3 + order - 1 + k] = order - 1 - k ;
-    w[3 + order - 1 + k] = 0.;
-
-    u[3 + 2 * (order - 1) + k] = order - 1 - k;
-    v[3 + 2 * (order - 1) + k] = 0.;
-    w[3 + 2 * (order - 1) + k] = 0.;
-  }
-
-  if (order > 2){
-    int nbdoftemp = nbdoftriangle(order - 3);
-    nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
-                      &w[3 + 3* (order - 1)]);
-    for (int k = 0; k < nbdoftemp; k++){
-      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
-    }
-  }
-  for (int k = 0; k < ndofT; k++){
-    u[k] = u[k] / order;
-    v[k] = v[k] / order;
-    w[k] = w[k] / order;
-  }
-}
-
-static void nodepositionface1(int order, double *u, double *v, double *w)
-{ // uw surface - orientation v0-v1-v3
-  int ndofT = nbdoftriangle(order);
-  if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
-  u[0] = 0.;    v[0]= 0.;  w[0] = 0.;
-  u[1] = order; v[1]= 0.;  w[1] = 0.;
-  u[2] = 0.;    v[2]= 0.;  w[2] = order;
-  // edges
-  for (int k = 0; k < (order - 1); k++){
-    u[3 + k] = k + 1;
-    v[3 + k] = 0.;
-    w[3 + k] = 0.;
-
-    u[3 + order - 1 + k] = order - 1 - k;
-    v[3 + order - 1 + k] = 0.;
-    w[3 + order - 1+ k ] = k + 1;
-
-    u[3 + 2 * (order - 1) + k] = 0. ;
-    v[3 + 2 * (order - 1) + k] = 0.;
-    w[3 + 2 * (order - 1) + k] = order - 1 - k;
-  }
-  if (order > 2){
-    int nbdoftemp = nbdoftriangle(order - 3);
-    nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
-      &w[3 + 3 * (order - 1)]);
-    for (int k = 0; k < nbdoftemp; k++){
-      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
-      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-    }
-  }
-  for (int k = 0; k < ndofT; k++){
-    u[k] = u[k] / order;
-    v[k] = v[k] / order;
-    w[k] = w[k] / order;
-  }
-}
-
-static void nodepositionface2(int order, double *u, double *v, double *w)
-{ // vw surface - orientation v0-v3-v2
-   int ndofT = nbdoftriangle(order);
-   if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
-
-   u[0]= 0.; v[0]= 0.;    w[0] = 0.;
-   u[1]= 0.; v[1]= 0.;    w[1] = order;
-   u[2]= 0.; v[2]= order; w[2] = 0.;
-   // edges
-   for (int k = 0; k < (order - 1); k++){
-
-     u[3 + k] = 0.;
-     v[3 + k] = 0.;
-     w[3 + k] = k + 1;
-
-     u[3 + order - 1 + k] = 0.;
-     v[3 + order - 1 + k] = k + 1;
-     w[3 + order - 1 + k] = order - 1 - k;
-
-     u[3 + 2 * (order - 1) + k] = 0.;
-     v[3 + 2 * (order - 1) + k] = order - 1 - k;
-     w[3 + 2 * (order - 1) + k] = 0.;
-   }
-   if (order > 2){
-     int nbdoftemp = nbdoftriangle(order - 3);
-     nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
-                       &w[3 + 3 * (order - 1)]);
-     for (int k = 0; k < nbdoftemp; k++){
-       u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
-       v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-       w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-     }
-   }
-   for (int k = 0; k < ndofT; k++){
-     u[k] = u[k] / order;
-     v[k] = v[k] / order;
-     w[k] = w[k] / order;
-   }
-}
-
-static void nodepositionface3(int order,  double *u,  double *v,  double *w)
-{ // uvw surface  - orientation v3-v1-v2
-   int ndofT = nbdoftriangle(order);
-   if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
-   u[0]= 0.;    v[0]= 0.;    w[0] = order;
-   u[1]= order; v[1]= 0.;    w[1] = 0.;
-   u[2]= 0.;    v[2]= order; w[2] = 0.;
-   // edges
-   for (int k = 0; k < (order - 1); k++){
-
-     u[3 + k] = k + 1;
-     v[3 + k] = 0.;
-     w[3 + k] = order - 1 - k;
-
-     u[3 + order - 1 + k] = order - 1 - k;
-     v[3 + order - 1 + k] = k + 1;
-     w[3 + order - 1 + k] = 0.;
-
-     u[3 + 2 * (order - 1) + k] = 0.;
-     v[3 + 2 * (order - 1) + k] = order - 1 - k;
-     w[3 + 2 * (order - 1) + k] = k + 1;
-   }
-   if (order > 2){
-     int nbdoftemp = nbdoftriangle(order - 3);
-     nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
-                       &w[3 + 3 * (order - 1)]);
-     for (int k = 0; k < nbdoftemp; k++){
-       u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-       v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-       w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
-     }
-   }
-   for (int k = 0; k < ndofT; k++){
-     u[k] = u[k] / order;
-     v[k] = v[k] / order;
-     w[k] = w[k] / order;
-   }
-}
-
-static fullMatrix<double> generateExponentsTetrahedron(int order)
-{
-  int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
-
-  fullMatrix<double> exponents(nbPoints, 3);
-
-  exponents(0, 0) = 0;
-  exponents(0, 1) = 0;
-  exponents(0, 2) = 0;
-
-  if (order > 0) {
-    exponents(1, 0) = order;
-    exponents(1, 1) = 0;
-    exponents(1, 2) = 0;
-
-    exponents(2, 0) = 0;
-    exponents(2, 1) = order;
-    exponents(2, 2) = 0;
-
-    exponents(3, 0) = 0;
-    exponents(3, 1) = 0;
-    exponents(3, 2) = order;
-
-    // edges e5 and e6 switched in original version, opposite direction
-    // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
-
-    if (order > 1) {
-      for (int k = 0; k < (order - 1); k++) {
-        exponents(4 + k, 0) = k + 1;
-        exponents(4 +      order - 1  + k, 0) = order - 1 - k;
-        exponents(4 + 2 * (order - 1) + k, 0) = 0;
-        exponents(4 + 3 * (order - 1) + k, 0) = 0;
-        // exponents(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
-        // exponents(4 + 5 * (order - 1) + k, 0) = 0.;
-        exponents(4 + 4 * (order - 1) + k, 0) = 0;
-        exponents(4 + 5 * (order - 1) + k, 0) = k+1;
-
-        exponents(4 + k, 1) = 0;
-        exponents(4 +      order - 1  + k, 1) = k + 1;
-        exponents(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
-        exponents(4 + 3 * (order - 1) + k, 1) = 0;
-        //         exponents(4 + 4 * (order - 1) + k, 1) = 0.;
-        //         exponents(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
-        exponents(4 + 4 * (order - 1) + k, 1) = k+1;
-        exponents(4 + 5 * (order - 1) + k, 1) = 0;
-
-        exponents(4 + k, 2) = 0;
-        exponents(4 +      order - 1  + k, 2) = 0;
-        exponents(4 + 2 * (order - 1) + k, 2) = 0;
-        exponents(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
-        exponents(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
-        exponents(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
-      }
-
-      if (order > 2) {
-        int ns = 4 + 6 * (order - 1);
-        int nbdofface = nbdoftriangle(order - 3);
-
-        double *u = new double[nbdofface];
-        double *v = new double[nbdofface];
-        double *w = new double[nbdofface];
-
-        nodepositionface0(order - 3, u, v, w);
-
-        // u-v plane
-
-        for (int i = 0; i < nbdofface; i++){
-          exponents(ns + i, 0) = u[i] * (order - 3) + 1;
-          exponents(ns + i, 1) = v[i] * (order - 3) + 1;
-          exponents(ns + i, 2) = w[i] * (order - 3);
-        }
-
-        ns = ns + nbdofface;
-
-        // u-w plane
-
-        nodepositionface1(order - 3, u, v, w);
-
-        for (int i=0; i < nbdofface; i++){
-          exponents(ns + i, 0) = u[i] * (order - 3) + 1;
-          exponents(ns + i, 1) = v[i] * (order - 3) ;
-          exponents(ns + i, 2) = w[i] * (order - 3) + 1;
-        }
-
-        // v-w plane
-
-        ns = ns + nbdofface;
-
-        nodepositionface2(order - 3, u, v, w);
-
-        for (int i = 0; i < nbdofface; i++){
-          exponents(ns + i, 0) = u[i] * (order - 3);
-          exponents(ns + i, 1) = v[i] * (order - 3) + 1;
-          exponents(ns + i, 2) = w[i] * (order - 3) + 1;
-        }
-
-        // u-v-w plane
-
-        ns = ns + nbdofface;
-
-        nodepositionface3(order - 3, u, v, w);
-
-        for (int i = 0; i < nbdofface; i++){
-          exponents(ns + i, 0) = u[i] * (order - 3) + 1;
-          exponents(ns + i, 1) = v[i] * (order - 3) + 1;
-          exponents(ns + i, 2) = w[i] * (order - 3) + 1;
-        }
-
-        ns = ns + nbdofface;
-
-        delete [] u;
-        delete [] v;
-        delete [] w;
-
-        if (order > 3) {
-
-          fullMatrix<double> interior = generateExponentsTetrahedron(order - 4);
-          for (int k = 0; k < interior.size1(); k++) {
-            exponents(ns + k, 0) = 1 + interior(k, 0);
-            exponents(ns + k, 1) = 1 + interior(k, 1);
-            exponents(ns + k, 2) = 1 + interior(k, 2);
-          }
-        }
-      }
-    }
-  }
-
-  return exponents;
-}
-
-static fullMatrix<double> generateExponentsPrism(int order)
-{
-  int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
-  fullMatrix<double> exponents(nbPoints, 3);
-
-  int index = 0;
-  fullMatrix<double> triExp = generateExponentsTriangle(order);
-  fullMatrix<double> lineExp = generate1DExponents(order);
-  // First exponents (i < 3) must relate to corner
-  for (int i = 0; i < triExp.size1(); i++) {
-    exponents(index,0) = triExp(i,0);
-    exponents(index,1) = triExp(i,1);
-    exponents(index,2) = lineExp(0,0);
-    index ++;
-    exponents(index,0) = triExp(i,0);
-    exponents(index,1) = triExp(i,1);
-    exponents(index,2) = lineExp(1,0);
-    index ++;
-  }
-  for (int j = 2; j <lineExp.size1() ; j++) {
-    for (int i = 0; i < triExp.size1(); i++) {
-      exponents(index,0) = triExp(i,0);
-      exponents(index,1) = triExp(i,1);
-      exponents(index,2) = lineExp(j,0);
-      index ++;
-    }
-  }
-
-  return exponents;
-}
-
-static fullMatrix<double> generateExponentsHex(int order)
-{
-  int nbPoints = (order+1) * (order+1) * (order+1);
-  fullMatrix<double> exponents(nbPoints, 3);
-
-  if (order == 0) {
-    exponents(0, 0) = .0;
-    return exponents;
-  }
-
-  int index = 0;
-  fullMatrix<double> lineExp = generate1DExponents(order);
-
-  for (int i = 0; i < 2; ++i) {
-    for (int j = 0; j < 2; ++j) {
-      for (int k = 0; k < 2; ++k) {
-        exponents(index, 0) = lineExp(i, 0);
-        exponents(index, 1) = lineExp(j, 0);
-        exponents(index, 2) = lineExp(k, 0);
-        ++index;
-      }
-    }
-  }
-
-  for (int i = 2; i < lineExp.size1(); ++i) {
-    for (int j = 0; j < 2; ++j) {
-      for (int k = 0; k < 2; ++k) {
-        exponents(index, 0) = lineExp(i, 0);
-        exponents(index, 1) = lineExp(j, 0);
-        exponents(index, 2) = lineExp(k, 0);
-        ++index;
-      }
-    }
-  }
-  for (int i = 0; i < lineExp.size1(); ++i) {
-    for (int j = 2; j < lineExp.size1(); ++j) {
-      for (int k = 0; k < 2; ++k) {
-        exponents(index, 0) = lineExp(i, 0);
-        exponents(index, 1) = lineExp(j, 0);
-        exponents(index, 2) = lineExp(k, 0);
-        ++index;
-      }
-    }
-  }
-  for (int i = 0; i < lineExp.size1(); ++i) {
-    for (int j = 0; j < lineExp.size1(); ++j) {
-      for (int k = 2; k < lineExp.size1(); ++k) {
-        exponents(index, 0) = lineExp(i, 0);
-        exponents(index, 1) = lineExp(j, 0);
-        exponents(index, 2) = lineExp(k, 0);
-        ++index;
-      }
-    }
-  }
-
-  return exponents;
-}
-
-// Sampling Points
-static fullMatrix<double> generate1DPoints(int order)
-{
-  fullMatrix<double> line(order + 1, 1);
-  line(0,0) = .0;
-  if (order > 0) {
-    line(1, 0) = 1.;
-    double dd = 1. / order;
-    for (int i = 2; i < order + 1; i++) line(i, 0) = dd * (i - 1);
-  }
-
-  return line;
-}
-
-static fullMatrix<double> generatePointsQuadRecur(int order, bool serendip)
-{
-  int nbPoints = serendip ? order*4 : (order+1)*(order+1);
-  fullMatrix<double> point(nbPoints, 2);
-
-  if (order > 0) {
-    point(0, 0) = -1;
-    point(0, 1) = -1;
-    point(1, 0) = 1;
-    point(1, 1) = -1;
-    point(2, 0) = 1;
-    point(2, 1) = 1;
-    point(3, 0) = -1;
-    point(3, 1) = 1;
-
-    if (order > 1) {
-      int index = 4;
-      const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
-      for (int iedge=0; iedge<4; iedge++) {
-        int p0 = edges[iedge][0];
-        int p1 = edges[iedge][1];
-        for (int i = 1; i < order; i++, index++) {
-          point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
-          point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
-        }
-      }
-      if (order >= 2 && !serendip) {
-        fullMatrix<double> inner = generatePointsQuadRecur(order - 2, false);
-        inner.scale(1. - 2./order);
-        point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
-      }
-    }
-  }
-  else {
-    point(0, 0) = 0;
-    point(0, 1) = 0;
-  }
-
-
-  return point;
-}
-
-static fullMatrix<double> generatePointsQuad(int order, bool serendip)
-{
-  fullMatrix<double>  point = generatePointsQuadRecur(order, serendip);
-  // rescale to [0,1] x [0,1]
-  for (int i=0;i<point.size1();i++){
-    point(i,0) = (1.+point(i,0))*.5;
-    point(i,1) = (1.+point(i,1))*.5;
-  }
-  return point;
-}
-
-static fullMatrix<double> generatePointsPrism(int order, bool serendip)
-{
-  const double prism18Pts[18][3] = {
-    {0, 0, 0}, // 0
-    {1, 0, 0}, // 1
-    {0, 1, 0}, // 2
-    {0, 0, 1},  // 3
-    {1, 0, 1},  // 4
-    {0, 1, 1},  // 5
-    {.5, 0, 0},  // 6
-    {0, .5, 0},  // 7
-    {0, 0, .5},  // 8
-    {.5, .5, 0},  // 9
-    {1, 0, .5},  // 10
-    {0, 1, .5},  // 11
-    {.5, 0, 1},  // 12
-    {0, .5, 1},  // 13
-    {.5, .5, 1},  // 14
-    {.5, 0, .5},  // 15
-    {0, .5, .5},  // 16
-    {.5, .5, .5},  // 17
-  };
-
-  int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
-  fullMatrix<double> point(nbPoints, 3);
-
-  int index = 0;
-
-  if (order == 2)
-    for (int i =0; i<18; i++)
-      for (int j=0; j<3;j++)
-        point(i,j) = prism18Pts[i][j];
-  else {
-    fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
-    fullMatrix<double> linePoints = generate1DPoints(order);
-    for (int j = 0; j <linePoints.size1() ; j++) {
-      for (int i = 0; i < triPoints.size1(); i++) {
-        point(index,0) = triPoints(i,0);
-        point(index,1) = triPoints(i,1);
-        point(index,2) = linePoints(j,0);
-        index ++;
-      }
-    }
-  }
-  return point;
-}
-
-static fullMatrix<double> generatePointsHex(int order, bool serendip)
-{
-  int nbPoints = (order+1) * (order+1) * (order+1);
-  fullMatrix<double> point(nbPoints, 3);
-
-  fullMatrix<double> linePoints = generate1DPoints(order);
-  int index = 0;
-
-  for (int i = 0; i < linePoints.size1(); ++i) {
-    for (int j = 0; j < linePoints.size1(); ++j) {
-      for (int k = 0; k < linePoints.size1(); ++k) {
-        point(index, 0) = linePoints(i, 0);
-        point(index, 1) = linePoints(j, 0);
-        point(index, 2) = linePoints(k, 0);
-        ++index;
-      }
-    }
-  }
-
-  return point;
-}
-
-// Sub Control Points
-static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
-{
-  std::vector< fullMatrix<double> > subPoints(2);
-
-  subPoints[0] = generate1DPoints(order);
-  subPoints[0].scale(.5);
-
-  subPoints[1].copy(subPoints[0]);
-  subPoints[1].add(.5);
-
-  return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
-{
-  std::vector< fullMatrix<double> > subPoints(4);
-  fullMatrix<double> prox;
-
-  subPoints[0] = gmshGeneratePointsTriangle(order, false);
-  subPoints[0].scale(.5);
-
-  subPoints[1].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[1], 0, 1);
-  prox.add(.5);
-
-  subPoints[2].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[2], 1, 1);
-  prox.add(.5);
-
-  subPoints[3].copy(subPoints[0]);
-  subPoints[3].scale(-1.);
-  subPoints[3].add(.5);
-
-  return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
-{
-  std::vector< fullMatrix<double> > subPoints(4);
-  fullMatrix<double> prox;
-
-  subPoints[0] = generatePointsQuad(order, false);
-  subPoints[0].scale(.5);
-
-  subPoints[1].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[1], 0, 1);
-  prox.add(.5);
-
-  subPoints[2].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[2], 1, 1);
-  prox.add(.5);
-
-  subPoints[3].copy(subPoints[1]);
-  prox.setAsProxy(subPoints[3], 1, 1);
-  prox.add(.5);
-
-  return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
-{
-  std::vector< fullMatrix<double> > subPoints(8);
-  fullMatrix<double> prox1;
-  fullMatrix<double> prox2;
-
-  subPoints[0] = gmshGeneratePointsTetrahedron(order, false);
-  subPoints[0].scale(.5);
-
-  subPoints[1].copy(subPoints[0]);
-  prox1.setAsProxy(subPoints[1], 0, 1);
-  prox1.add(.5);
-
-  subPoints[2].copy(subPoints[0]);
-  prox1.setAsProxy(subPoints[2], 1, 1);
-  prox1.add(.5);
-
-  subPoints[3].copy(subPoints[0]);
-  prox1.setAsProxy(subPoints[3], 2, 1);
-  prox1.add(.5);
-
-  // u := .5-u-w
-  // v := .5-v-w
-  // w := w
-  subPoints[4].copy(subPoints[0]);
-  prox1.setAsProxy(subPoints[4], 0, 2);
-  prox1.scale(-1.);
-  prox1.add(.5);
-  prox1.setAsProxy(subPoints[4], 0, 1);
-  prox2.setAsProxy(subPoints[4], 2, 1);
-  prox1.add(prox2, -1.);
-  prox1.setAsProxy(subPoints[4], 1, 1);
-  prox1.add(prox2, -1.);
-
-  // u := u
-  // v := .5-v
-  // w := w+v
-  subPoints[5].copy(subPoints[0]);
-  prox1.setAsProxy(subPoints[5], 2, 1);
-  prox2.setAsProxy(subPoints[5], 1, 1);
-  prox1.add(prox2);
-  prox2.scale(-1.);
-  prox2.add(.5);
-
-  // u := .5-u
-  // v := v
-  // w := w+u
-  subPoints[6].copy(subPoints[0]);
-  prox1.setAsProxy(subPoints[6], 2, 1);
-  prox2.setAsProxy(subPoints[6], 0, 1);
-  prox1.add(prox2);
-  prox2.scale(-1.);
-  prox2.add(.5);
-
-  // u := u+w
-  // v := v+w
-  // w := .5-w
-  subPoints[7].copy(subPoints[0]);
-  prox1.setAsProxy(subPoints[7], 0, 1);
-  prox2.setAsProxy(subPoints[7], 2, 1);
-  prox1.add(prox2);
-  prox1.setAsProxy(subPoints[7], 1, 1);
-  prox1.add(prox2);
-  prox2.scale(-1.);
-  prox2.add(.5);
-
-
-  return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
-{
-  std::vector< fullMatrix<double> > subPoints(8);
-  fullMatrix<double> prox;
-
-  subPoints[0] = generatePointsPrism(order, false);
-  subPoints[0].scale(.5);
-
-  subPoints[1].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[1], 0, 1);
-  prox.add(.5);
-
-  subPoints[2].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[2], 1, 1);
-  prox.add(.5);
-
-  subPoints[3].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[3], 0, 2);
-  prox.scale(-1.);
-  prox.add(.5);
-
-  subPoints[4].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[4], 2, 1);
-  prox.add(.5);
-
-  subPoints[5].copy(subPoints[1]);
-  prox.setAsProxy(subPoints[5], 2, 1);
-  prox.add(.5);
-
-  subPoints[6].copy(subPoints[2]);
-  prox.setAsProxy(subPoints[6], 2, 1);
-  prox.add(.5);
-
-  subPoints[7].copy(subPoints[3]);
-  prox.setAsProxy(subPoints[7], 2, 1);
-  prox.add(.5);
-
-  return subPoints;
-}
-
-static std::vector< fullMatrix<double> > generateSubPointsHex(int order)
-{
-  std::vector< fullMatrix<double> > subPoints(8);
-  fullMatrix<double> prox;
-
-  subPoints[0] = generatePointsHex(order, false);
-  subPoints[0].scale(.5);
-
-  subPoints[1].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[1], 0, 1);
-  prox.add(.5);
-
-  subPoints[2].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[2], 1, 1);
-  prox.add(.5);
-
-  subPoints[3].copy(subPoints[1]);
-  prox.setAsProxy(subPoints[3], 1, 1);
-  prox.add(.5);
-
-  subPoints[4].copy(subPoints[0]);
-  prox.setAsProxy(subPoints[4], 2, 1);
-  prox.add(.5);
-
-  subPoints[5].copy(subPoints[1]);
-  prox.setAsProxy(subPoints[5], 2, 1);
-  prox.add(.5);
-
-  subPoints[6].copy(subPoints[2]);
-  prox.setAsProxy(subPoints[6], 2, 1);
-  prox.add(.5);
-
-  subPoints[7].copy(subPoints[3]);
-  prox.setAsProxy(subPoints[7], 2, 1);
-  prox.add(.5);
-
-  return subPoints;
-}
-
-// Matrices generation
-static int nChoosek(int n, int k)
-{
-  if (n < k || k < 0) {
-    Msg::Error("Wrong argument for combination.");
-    return 1;
-  }
-
-  if (k > n/2) k = n-k;
-  if (k == 1)
-    return n;
-  if (k == 0)
-    return 1;
-
-  int c = 1;
-  for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
-  return c;
-}
-
-
-static fullMatrix<double> generateBez2LagMatrix
-  (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
-   int order, int dimSimplex)
-{
-
-  if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
-    Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
-      exponent.size1(),point.size1(),
-      exponent.size2(),point.size2());
-    return fullMatrix<double>(1, 1);
-  }
-
-  int ndofs = exponent.size1();
-  int dim = exponent.size2();
-
-  fullMatrix<double> bez2Lag(ndofs, ndofs);
-  for (int i = 0; i < ndofs; i++) {
-    for (int j = 0; j < ndofs; j++) {
-      double dd = 1.;
-
-      {
-        double pointCompl = 1.;
-        int exponentCompl = order;
-        for (int k = 0; k < dimSimplex; k++) {
-          dd *= nChoosek(exponentCompl, (int) exponent(i, k))
-            * pow(point(j, k), exponent(i, k));
-          pointCompl -= point(j, k);
-          exponentCompl -= (int) exponent(i, k);
-        }
-        dd *= pow(pointCompl, exponentCompl);
-      }
-
-      for (int k = dimSimplex; k < dim; k++)
-        dd *= nChoosek(order, (int) exponent(i, k))
-            * pow(point(j, k), exponent(i, k))
-            * pow(1. - point(j, k), order - exponent(i, k));
-
-      bez2Lag(j, i) = dd;
-    }
-  }
-  return bez2Lag;
-}
-
-
-
-static fullMatrix<double> generateSubDivisor
-  (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
-   const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
-{
-  if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
-    Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
-      exponents.size1(), lag2Bez.size1(),
-      exponents.size1(), lag2Bez.size2());
-    return fullMatrix<double>(1, 1);
-  }
-
-  int nbPts = lag2Bez.size1();
-  int nbSubPts = nbPts * subPoints.size();
-
-  fullMatrix<double> intermediate2(nbPts, nbPts);
-  fullMatrix<double> subDivisor(nbSubPts, nbPts);
-
-  for (unsigned int i = 0; i < subPoints.size(); i++) {
-    fullMatrix<double> intermediate1 =
-      generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
-    lag2Bez.mult(intermediate1, intermediate2);
-    subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
-  }
-  return subDivisor;
-}
-
-
-
-// Convert Bezier points to Lagrange points
-static fullMatrix<double> bez2LagPoints(bool dim1, bool dim2, bool dim3, const fullMatrix<double> &bezierPoints)
-{
-
-  // FIXME BUG for 2nd order quads we seem to try to access bezierPoints(i, 2);
-  // but bezierPoints only has 2 columns
-
-  const int numPt = bezierPoints.size1();
-  fullMatrix<double> lagPoints(numPt,3);
-
-  for (int i=0; i<numPt; i++) {
-    lagPoints(i,0) = dim1 ? -1. + 2*bezierPoints(i,0) : bezierPoints(i,0);
-    lagPoints(i,1) = dim2 ? -1. + 2*bezierPoints(i,1) : bezierPoints(i,1);
-    lagPoints(i,2) = dim3 ? -1. + 2*bezierPoints(i,2) : bezierPoints(i,2);
-  }
-
-  return lagPoints;
-
-}
-
-
-
-std::map<int, bezierBasis> bezierBasis::_bbs;
-const bezierBasis *bezierBasis::find(int tag)
-{
-  std::map<int, bezierBasis>::const_iterator it = _bbs.find(tag);
-  if (it != _bbs.end())
-    return &it->second;
-
-  bezierBasis &B = _bbs[tag];
-  B.order = MElement::OrderFromTag(tag);
-
-  if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
-    B.dim = 3;
-    B.numLagPts = 5;
-    B.bezierPoints = gmshGeneratePointsPyramid(B.order,false);
-    B.lagPoints = B.bezierPoints;
-    B.matrixLag2Bez.resize(B.bezierPoints.size1(),B.bezierPoints.size1(),0.);
-    B.matrixBez2Lag.resize(B.bezierPoints.size1(),B.bezierPoints.size1(),0.);
-    for (int i=0; i<B.matrixBez2Lag.size1(); ++i) {
-      B.matrixLag2Bez(i,i) = 1.;
-      B.matrixBez2Lag(i,i) = 1.;
-    }
-    // TODO: Implement subdidivsor
-  }
-  else {
-    int dimSimplex;
-    std::vector< fullMatrix<double> > subPoints;
-    switch (MElement::ParentTypeFromTag(tag)) {
-      case TYPE_PNT :
-        B.dim = 0;
-        B.numLagPts = 1;
-        B.exponents = generate1DExponents(0);
-        B.bezierPoints = generate1DPoints(0);
-        B.lagPoints = B.bezierPoints;
-        dimSimplex = 0;
-        break;
-      case TYPE_LIN : {
-        B.dim = 1;
-        B.numLagPts = 2;
-        B.exponents = generate1DExponents(B.order);
-        B.bezierPoints = generate1DPoints(B.order);
-        B.lagPoints = bez2LagPoints(true,false,false,B.bezierPoints);
-        dimSimplex = 0;
-        subPoints = generateSubPointsLine(0);
-        break;
-      }
-      case TYPE_TRI : {
-        B.dim = 2;
-        B.numLagPts = 3;
-        B.exponents = generateExponentsTriangle(B.order);
-        B.bezierPoints = gmshGeneratePointsTriangle(B.order,false);
-        B.lagPoints = B.bezierPoints;
-        dimSimplex = 2;
-        subPoints = generateSubPointsTriangle(B.order);
-        break;
-      }
-      case TYPE_QUA : {
-        B.dim = 2;
-        B.numLagPts = 4;
-        B.exponents = generateExponentsQuad(B.order);
-        B.bezierPoints = generatePointsQuad(B.order,false);
-        B.lagPoints = bez2LagPoints(true,true,false,B.bezierPoints);
-        dimSimplex = 0;
-        subPoints = generateSubPointsQuad(B.order);
-        //      B.points.print("points");
-        break;
-      }
-      case TYPE_TET : {
-        B.dim = 3;
-        B.numLagPts = 4;
-        B.exponents = generateExponentsTetrahedron(B.order);
-        B.bezierPoints = gmshGeneratePointsTetrahedron(B.order,false);
-        B.lagPoints = B.bezierPoints;
-        dimSimplex = 3;
-        subPoints = generateSubPointsTetrahedron(B.order);
-        break;
-      }
-      case TYPE_PRI : {
-        B.dim = 3;
-        B.numLagPts = 6;
-        B.exponents = generateExponentsPrism(B.order);
-        B.bezierPoints = generatePointsPrism(B.order, false);
-        B.lagPoints = bez2LagPoints(false,false,true,B.bezierPoints);
-        dimSimplex = 2;
-        subPoints = generateSubPointsPrism(B.order);
-        break;
-      }
-      case TYPE_HEX : {
-        B.dim = 3;
-        B.numLagPts = 8;
-        B.exponents = generateExponentsHex(B.order);
-        B.bezierPoints = generatePointsHex(B.order, false);
-        B.lagPoints = bez2LagPoints(true,true,true,B.bezierPoints);
-        dimSimplex = 0;
-        subPoints = generateSubPointsHex(B.order);
-        break;
-      }
-      default : {
-        Msg::Error("Unknown function space %d: reverting to TET_1", tag);
-        B.dim = 3;
-        B.numLagPts = 4;
-        B.exponents = generateExponentsTetrahedron(0);
-        B.bezierPoints = gmshGeneratePointsTetrahedron(0, false);
-        B.lagPoints = B.bezierPoints;
-        dimSimplex = 3;
-        subPoints = generateSubPointsTetrahedron(0);
-        break;
-      }
-    }
-    B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.bezierPoints, B.order, dimSimplex);
-    B.matrixBez2Lag.invert(B.matrixLag2Bez);
-    B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, B.order, dimSimplex);
-    B.numDivisions = (int) pow(2., (int) B.bezierPoints.size2());
-  }
-
-  return &B;
-}
-
-
 
 JacobianBasis::JacobianBasis(int tag)
 {
@@ -1113,10 +51,10 @@ JacobianBasis::JacobianBasis(int tag)
   }
 
   // Store Bezier basis
-  bezier = bezierBasis::find(jacType);
+  bezier = BasisFactory::getBezierBasis(jacType);
 
   // Store shape function gradients of mapping at Jacobian nodes
-  const nodalBasis *mapBasis = BasisFactory::create(tag);
+  const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
   fullMatrix<double> allDPsi;
   mapBasis->df(getPoints(), allDPsi);
   numJacNodes = getPoints().size1();
@@ -1132,7 +70,7 @@ JacobianBasis::JacobianBasis(int tag)
     }
 
   // Compute matrix for lifting from primary Jacobian basis to Jacobian basis
-  const nodalBasis *primJacBasis = BasisFactory::create(primJacType);
+  const nodalBasis *primJacBasis = BasisFactory::getNodalBasis(primJacType);
   numPrimJacNodes = primJacBasis->getNumShapeFunctions();
   matrixPrimJac2Jac.resize(numJacNodes,numPrimJacNodes);
   primJacBasis->f(getPoints(),matrixPrimJac2Jac);
@@ -1140,7 +78,7 @@ JacobianBasis::JacobianBasis(int tag)
   // Compute shape function gradients of primary mapping at barycenter,
   // in order to compute normal to straight element
   const int primMapType = polynomialBasis::getTag(parentType, 1, false);
-  const nodalBasis *primMapBasis = BasisFactory::create(primMapType);
+  const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(primMapType);
   numPrimMapNodes = primMapBasis->getNumShapeFunctions();
   double xBar = 0., yBar = 0., zBar = 0.;
   for (int i=0; i<numPrimMapNodes; i++) {
@@ -1168,21 +106,6 @@ JacobianBasis::JacobianBasis(int tag)
 
 
 
-std::map<int, JacobianBasis*> JacobianBasis::_fs;
-const JacobianBasis *JacobianBasis::find(int tag)
-{
-
-  std::map<int, JacobianBasis*>::const_iterator it = _fs.find(tag);
-  if (it != _fs.end()) return it->second;
-
-  JacobianBasis *B = new JacobianBasis(tag);
-  _fs.insert(std::make_pair(tag, B));
-  return B;
-
-}
-
-
-
 // Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value)
 double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
 {
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 8901c1d9a5a32f265916ec94abb4af1a70d68e43..2e8f52b23c84600c3576e7f5bdfe34e2495bde03 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -8,28 +8,12 @@
 
 #include <map>
 #include <vector>
+#include "bezierBasis.h"
 #include "fullMatrix.h"
 
-class MElement;
-
-class bezierBasis {
- private :
-  static std::map<int, bezierBasis> _bbs;
- public :
-  int dim, order;
-  int numLagPts;
-  int numDivisions;
-  // the 'numLagPts' first exponents are related to 'Lagrangian' values
-  fullMatrix<double> exponents; // exponents of Bezier FF
-  fullMatrix<double> bezierPoints, lagPoints; // sampling point
-  fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
-  fullMatrix<double> subDivisor;
-  static const bezierBasis *find(int);
-};
 
 class JacobianBasis {
  private:
-  static std::map<int, JacobianBasis*> _fs;
   const bezierBasis *bezier;
   fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
   fullVector<double> primGradShapeBarX, primGradShapeBarY, primGradShapeBarZ;
@@ -37,7 +21,6 @@ class JacobianBasis {
   int numJacNodes, numMapNodes, numPrimJacNodes, numPrimMapNodes;
   inline const fullMatrix<double> &getPoints() const { return bezier->lagPoints; }
  public :
-  static const JacobianBasis *find(int);
   JacobianBasis(int tag);
   inline int getNumJacNodes() const { return numJacNodes; }
   inline int getNumMapNodes() const { return numMapNodes; }
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f0592dc93802fe471eb8190332e45f609e7f7bc
--- /dev/null
+++ b/Numeric/bezierBasis.cpp
@@ -0,0 +1,1069 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+#include <vector>
+#include "polynomialBasis.h"
+#include "pyramidalBasis.h"
+#include "pointsGenerators.h"
+#include "BasisFactory.h"
+#include "Numeric.h"
+
+// Bezier Exponents
+static fullMatrix<double> generate1DExponents(int order)
+{
+  fullMatrix<double> exponents(order + 1, 1);
+  exponents(0, 0) = 0;
+  if (order > 0) {
+    exponents(1, 0) = order;
+    for (int i = 2; i < order + 1; i++) exponents(i, 0) = i - 1;
+  }
+
+  return exponents;
+}
+
+static fullMatrix<double> generateExponentsTriangle(int order)
+{
+  int nbPoints = (order + 1) * (order + 2) / 2;
+  fullMatrix<double> exponents(nbPoints, 2);
+
+  exponents(0, 0) = 0;
+  exponents(0, 1) = 0;
+
+  if (order > 0) {
+    exponents(1, 0) = order;
+    exponents(1, 1) = 0;
+    exponents(2, 0) = 0;
+    exponents(2, 1) = order;
+
+    if (order > 1) {
+      int index = 3, ksi = 0, eta = 0;
+
+      for (int i = 0; i < order - 1; i++, index++) {
+        ksi++;
+        exponents(index, 0) = ksi;
+        exponents(index, 1) = eta;
+      }
+
+      ksi = order;
+
+      for (int i = 0; i < order - 1; i++, index++) {
+        ksi--;
+        eta++;
+        exponents(index, 0) = ksi;
+        exponents(index, 1) = eta;
+      }
+
+      eta = order;
+      ksi = 0;
+
+      for (int i = 0; i < order - 1; i++, index++) {
+        eta--;
+        exponents(index, 0) = ksi;
+        exponents(index, 1) = eta;
+      }
+
+      if (order > 2) {
+        fullMatrix<double> inner = generateExponentsTriangle(order - 3);
+        inner.add(1);
+        exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+      }
+    }
+  }
+
+  return exponents;
+}
+static fullMatrix<double> generateExponentsQuad(int order)
+{
+  int nbPoints = (order+1)*(order+1);
+  fullMatrix<double> exponents(nbPoints, 2);
+
+  exponents(0, 0) = 0;
+  exponents(0, 1) = 0;
+
+  if (order > 0) {
+    exponents(1, 0) = order;
+    exponents(1, 1) = 0;
+    exponents(2, 0) = order;
+    exponents(2, 1) = order;
+    exponents(3, 0) = 0;
+    exponents(3, 1) = order;
+
+    if (order > 1) {
+      int index = 4;
+      const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
+      for (int iedge=0; iedge<4; iedge++) {
+        int p0 = edges[iedge][0];
+        int p1 = edges[iedge][1];
+        for (int i = 1; i < order; i++, index++) {
+          exponents(index, 0) = exponents(p0, 0) + i*(exponents(p1,0)-exponents(p0,0))/order;
+          exponents(index, 1) = exponents(p0, 1) + i*(exponents(p1,1)-exponents(p0,1))/order;
+        }
+      }
+      if (order > 2) {
+        fullMatrix<double> inner = generateExponentsQuad(order - 2);
+        inner.add(1);
+        exponents.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+      }
+    }
+  }
+
+  //  exponents.print("expq");
+
+  return exponents;
+}
+static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
+
+static void nodepositionface0(int order, double *u, double *v, double *w)
+{ // uv surface - orientation v0-v2-v1
+  int ndofT = nbdoftriangle(order);
+  if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+  u[0]= 0.;    v[0]= 0.;    w[0] = 0.;
+  u[1]= 0.;    v[1]= order; w[1] = 0.;
+  u[2]= order; v[2]= 0.;    w[2] = 0.;
+
+  // edges
+  for (int k = 0; k < (order - 1); k++){
+    u[3 + k] = 0.;
+    v[3 + k] = k + 1;
+    w[3 + k] = 0.;
+
+    u[3 + order - 1 + k] = k + 1;
+    v[3 + order - 1 + k] = order - 1 - k ;
+    w[3 + order - 1 + k] = 0.;
+
+    u[3 + 2 * (order - 1) + k] = order - 1 - k;
+    v[3 + 2 * (order - 1) + k] = 0.;
+    w[3 + 2 * (order - 1) + k] = 0.;
+  }
+
+  if (order > 2){
+    int nbdoftemp = nbdoftriangle(order - 3);
+    nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+                      &w[3 + 3* (order - 1)]);
+    for (int k = 0; k < nbdoftemp; k++){
+      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3);
+    }
+  }
+  for (int k = 0; k < ndofT; k++){
+    u[k] = u[k] / order;
+    v[k] = v[k] / order;
+    w[k] = w[k] / order;
+  }
+}
+
+static void nodepositionface1(int order, double *u, double *v, double *w)
+{ // uw surface - orientation v0-v1-v3
+  int ndofT = nbdoftriangle(order);
+  if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+  u[0] = 0.;    v[0]= 0.;  w[0] = 0.;
+  u[1] = order; v[1]= 0.;  w[1] = 0.;
+  u[2] = 0.;    v[2]= 0.;  w[2] = order;
+  // edges
+  for (int k = 0; k < (order - 1); k++){
+    u[3 + k] = k + 1;
+    v[3 + k] = 0.;
+    w[3 + k] = 0.;
+
+    u[3 + order - 1 + k] = order - 1 - k;
+    v[3 + order - 1 + k] = 0.;
+    w[3 + order - 1+ k ] = k + 1;
+
+    u[3 + 2 * (order - 1) + k] = 0. ;
+    v[3 + 2 * (order - 1) + k] = 0.;
+    w[3 + 2 * (order - 1) + k] = order - 1 - k;
+  }
+  if (order > 2){
+    int nbdoftemp = nbdoftriangle(order - 3);
+    nodepositionface1(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order -1 )],
+      &w[3 + 3 * (order - 1)]);
+    for (int k = 0; k < nbdoftemp; k++){
+      u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+      v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
+      w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+    }
+  }
+  for (int k = 0; k < ndofT; k++){
+    u[k] = u[k] / order;
+    v[k] = v[k] / order;
+    w[k] = w[k] / order;
+  }
+}
+
+static void nodepositionface2(int order, double *u, double *v, double *w)
+{ // vw surface - orientation v0-v3-v2
+   int ndofT = nbdoftriangle(order);
+   if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
+
+   u[0]= 0.; v[0]= 0.;    w[0] = 0.;
+   u[1]= 0.; v[1]= 0.;    w[1] = order;
+   u[2]= 0.; v[2]= order; w[2] = 0.;
+   // edges
+   for (int k = 0; k < (order - 1); k++){
+
+     u[3 + k] = 0.;
+     v[3 + k] = 0.;
+     w[3 + k] = k + 1;
+
+     u[3 + order - 1 + k] = 0.;
+     v[3 + order - 1 + k] = k + 1;
+     w[3 + order - 1 + k] = order - 1 - k;
+
+     u[3 + 2 * (order - 1) + k] = 0.;
+     v[3 + 2 * (order - 1) + k] = order - 1 - k;
+     w[3 + 2 * (order - 1) + k] = 0.;
+   }
+   if (order > 2){
+     int nbdoftemp = nbdoftriangle(order - 3);
+     nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+                       &w[3 + 3 * (order - 1)]);
+     for (int k = 0; k < nbdoftemp; k++){
+       u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
+       v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+       w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+     }
+   }
+   for (int k = 0; k < ndofT; k++){
+     u[k] = u[k] / order;
+     v[k] = v[k] / order;
+     w[k] = w[k] / order;
+   }
+}
+
+static void nodepositionface3(int order,  double *u,  double *v,  double *w)
+{ // uvw surface  - orientation v3-v1-v2
+   int ndofT = nbdoftriangle(order);
+   if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
+
+   u[0]= 0.;    v[0]= 0.;    w[0] = order;
+   u[1]= order; v[1]= 0.;    w[1] = 0.;
+   u[2]= 0.;    v[2]= order; w[2] = 0.;
+   // edges
+   for (int k = 0; k < (order - 1); k++){
+
+     u[3 + k] = k + 1;
+     v[3 + k] = 0.;
+     w[3 + k] = order - 1 - k;
+
+     u[3 + order - 1 + k] = order - 1 - k;
+     v[3 + order - 1 + k] = k + 1;
+     w[3 + order - 1 + k] = 0.;
+
+     u[3 + 2 * (order - 1) + k] = 0.;
+     v[3 + 2 * (order - 1) + k] = order - 1 - k;
+     w[3 + 2 * (order - 1) + k] = k + 1;
+   }
+   if (order > 2){
+     int nbdoftemp = nbdoftriangle(order - 3);
+     nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+                       &w[3 + 3 * (order - 1)]);
+     for (int k = 0; k < nbdoftemp; k++){
+       u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+       v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+       w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+     }
+   }
+   for (int k = 0; k < ndofT; k++){
+     u[k] = u[k] / order;
+     v[k] = v[k] / order;
+     w[k] = w[k] / order;
+   }
+}
+
+static fullMatrix<double> generateExponentsTetrahedron(int order)
+{
+  int nbPoints = (order + 1) * (order + 2) * (order + 3) / 6;
+
+  fullMatrix<double> exponents(nbPoints, 3);
+
+  exponents(0, 0) = 0;
+  exponents(0, 1) = 0;
+  exponents(0, 2) = 0;
+
+  if (order > 0) {
+    exponents(1, 0) = order;
+    exponents(1, 1) = 0;
+    exponents(1, 2) = 0;
+
+    exponents(2, 0) = 0;
+    exponents(2, 1) = order;
+    exponents(2, 2) = 0;
+
+    exponents(3, 0) = 0;
+    exponents(3, 1) = 0;
+    exponents(3, 2) = order;
+
+    // edges e5 and e6 switched in original version, opposite direction
+    // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
+
+    if (order > 1) {
+      for (int k = 0; k < (order - 1); k++) {
+        exponents(4 + k, 0) = k + 1;
+        exponents(4 +      order - 1  + k, 0) = order - 1 - k;
+        exponents(4 + 2 * (order - 1) + k, 0) = 0;
+        exponents(4 + 3 * (order - 1) + k, 0) = 0;
+        // exponents(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
+        // exponents(4 + 5 * (order - 1) + k, 0) = 0.;
+        exponents(4 + 4 * (order - 1) + k, 0) = 0;
+        exponents(4 + 5 * (order - 1) + k, 0) = k+1;
+
+        exponents(4 + k, 1) = 0;
+        exponents(4 +      order - 1  + k, 1) = k + 1;
+        exponents(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+        exponents(4 + 3 * (order - 1) + k, 1) = 0;
+        //         exponents(4 + 4 * (order - 1) + k, 1) = 0.;
+        //         exponents(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
+        exponents(4 + 4 * (order - 1) + k, 1) = k+1;
+        exponents(4 + 5 * (order - 1) + k, 1) = 0;
+
+        exponents(4 + k, 2) = 0;
+        exponents(4 +      order - 1  + k, 2) = 0;
+        exponents(4 + 2 * (order - 1) + k, 2) = 0;
+        exponents(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
+        exponents(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+        exponents(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
+      }
+
+      if (order > 2) {
+        int ns = 4 + 6 * (order - 1);
+        int nbdofface = nbdoftriangle(order - 3);
+
+        double *u = new double[nbdofface];
+        double *v = new double[nbdofface];
+        double *w = new double[nbdofface];
+
+        nodepositionface0(order - 3, u, v, w);
+
+        // u-v plane
+
+        for (int i = 0; i < nbdofface; i++){
+          exponents(ns + i, 0) = u[i] * (order - 3) + 1;
+          exponents(ns + i, 1) = v[i] * (order - 3) + 1;
+          exponents(ns + i, 2) = w[i] * (order - 3);
+        }
+
+        ns = ns + nbdofface;
+
+        // u-w plane
+
+        nodepositionface1(order - 3, u, v, w);
+
+        for (int i=0; i < nbdofface; i++){
+          exponents(ns + i, 0) = u[i] * (order - 3) + 1;
+          exponents(ns + i, 1) = v[i] * (order - 3) ;
+          exponents(ns + i, 2) = w[i] * (order - 3) + 1;
+        }
+
+        // v-w plane
+
+        ns = ns + nbdofface;
+
+        nodepositionface2(order - 3, u, v, w);
+
+        for (int i = 0; i < nbdofface; i++){
+          exponents(ns + i, 0) = u[i] * (order - 3);
+          exponents(ns + i, 1) = v[i] * (order - 3) + 1;
+          exponents(ns + i, 2) = w[i] * (order - 3) + 1;
+        }
+
+        // u-v-w plane
+
+        ns = ns + nbdofface;
+
+        nodepositionface3(order - 3, u, v, w);
+
+        for (int i = 0; i < nbdofface; i++){
+          exponents(ns + i, 0) = u[i] * (order - 3) + 1;
+          exponents(ns + i, 1) = v[i] * (order - 3) + 1;
+          exponents(ns + i, 2) = w[i] * (order - 3) + 1;
+        }
+
+        ns = ns + nbdofface;
+
+        delete [] u;
+        delete [] v;
+        delete [] w;
+
+        if (order > 3) {
+
+          fullMatrix<double> interior = generateExponentsTetrahedron(order - 4);
+          for (int k = 0; k < interior.size1(); k++) {
+            exponents(ns + k, 0) = 1 + interior(k, 0);
+            exponents(ns + k, 1) = 1 + interior(k, 1);
+            exponents(ns + k, 2) = 1 + interior(k, 2);
+          }
+        }
+      }
+    }
+  }
+
+  return exponents;
+}
+
+static fullMatrix<double> generateExponentsPrism(int order)
+{
+  int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+  fullMatrix<double> exponents(nbPoints, 3);
+
+  int index = 0;
+  fullMatrix<double> triExp = generateExponentsTriangle(order);
+  fullMatrix<double> lineExp = generate1DExponents(order);
+  // First exponents (i < 3) must relate to corner
+  for (int i = 0; i < triExp.size1(); i++) {
+    exponents(index,0) = triExp(i,0);
+    exponents(index,1) = triExp(i,1);
+    exponents(index,2) = lineExp(0,0);
+    index ++;
+    exponents(index,0) = triExp(i,0);
+    exponents(index,1) = triExp(i,1);
+    exponents(index,2) = lineExp(1,0);
+    index ++;
+  }
+  for (int j = 2; j <lineExp.size1() ; j++) {
+    for (int i = 0; i < triExp.size1(); i++) {
+      exponents(index,0) = triExp(i,0);
+      exponents(index,1) = triExp(i,1);
+      exponents(index,2) = lineExp(j,0);
+      index ++;
+    }
+  }
+
+  return exponents;
+}
+
+static fullMatrix<double> generateExponentsHex(int order)
+{
+  int nbPoints = (order+1) * (order+1) * (order+1);
+  fullMatrix<double> exponents(nbPoints, 3);
+
+  if (order == 0) {
+    exponents(0, 0) = .0;
+    return exponents;
+  }
+
+  int index = 0;
+  fullMatrix<double> lineExp = generate1DExponents(order);
+
+  for (int i = 0; i < 2; ++i) {
+    for (int j = 0; j < 2; ++j) {
+      for (int k = 0; k < 2; ++k) {
+        exponents(index, 0) = lineExp(i, 0);
+        exponents(index, 1) = lineExp(j, 0);
+        exponents(index, 2) = lineExp(k, 0);
+        ++index;
+      }
+    }
+  }
+
+  for (int i = 2; i < lineExp.size1(); ++i) {
+    for (int j = 0; j < 2; ++j) {
+      for (int k = 0; k < 2; ++k) {
+        exponents(index, 0) = lineExp(i, 0);
+        exponents(index, 1) = lineExp(j, 0);
+        exponents(index, 2) = lineExp(k, 0);
+        ++index;
+      }
+    }
+  }
+  for (int i = 0; i < lineExp.size1(); ++i) {
+    for (int j = 2; j < lineExp.size1(); ++j) {
+      for (int k = 0; k < 2; ++k) {
+        exponents(index, 0) = lineExp(i, 0);
+        exponents(index, 1) = lineExp(j, 0);
+        exponents(index, 2) = lineExp(k, 0);
+        ++index;
+      }
+    }
+  }
+  for (int i = 0; i < lineExp.size1(); ++i) {
+    for (int j = 0; j < lineExp.size1(); ++j) {
+      for (int k = 2; k < lineExp.size1(); ++k) {
+        exponents(index, 0) = lineExp(i, 0);
+        exponents(index, 1) = lineExp(j, 0);
+        exponents(index, 2) = lineExp(k, 0);
+        ++index;
+      }
+    }
+  }
+
+  return exponents;
+}
+
+// Sampling Points
+static fullMatrix<double> generate1DPoints(int order)
+{
+  fullMatrix<double> line(order + 1, 1);
+  line(0,0) = .0;
+  if (order > 0) {
+    line(1, 0) = 1.;
+    double dd = 1. / order;
+    for (int i = 2; i < order + 1; i++) line(i, 0) = dd * (i - 1);
+  }
+
+  return line;
+}
+
+static fullMatrix<double> generatePointsQuadRecur(int order, bool serendip)
+{
+  int nbPoints = serendip ? order*4 : (order+1)*(order+1);
+  fullMatrix<double> point(nbPoints, 2);
+
+  if (order > 0) {
+    point(0, 0) = -1;
+    point(0, 1) = -1;
+    point(1, 0) = 1;
+    point(1, 1) = -1;
+    point(2, 0) = 1;
+    point(2, 1) = 1;
+    point(3, 0) = -1;
+    point(3, 1) = 1;
+
+    if (order > 1) {
+      int index = 4;
+      const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
+      for (int iedge=0; iedge<4; iedge++) {
+        int p0 = edges[iedge][0];
+        int p1 = edges[iedge][1];
+        for (int i = 1; i < order; i++, index++) {
+          point(index, 0) = point(p0, 0) + i*(point(p1,0)-point(p0,0))/order;
+          point(index, 1) = point(p0, 1) + i*(point(p1,1)-point(p0,1))/order;
+        }
+      }
+      if (order >= 2 && !serendip) {
+        fullMatrix<double> inner = generatePointsQuadRecur(order - 2, false);
+        inner.scale(1. - 2./order);
+        point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
+      }
+    }
+  }
+  else {
+    point(0, 0) = 0;
+    point(0, 1) = 0;
+  }
+
+
+  return point;
+}
+
+static fullMatrix<double> generatePointsQuad(int order, bool serendip)
+{
+  fullMatrix<double>  point = generatePointsQuadRecur(order, serendip);
+  // rescale to [0,1] x [0,1]
+  for (int i=0;i<point.size1();i++){
+    point(i,0) = (1.+point(i,0))*.5;
+    point(i,1) = (1.+point(i,1))*.5;
+  }
+  return point;
+}
+
+static fullMatrix<double> generatePointsPrism(int order, bool serendip)
+{
+  const double prism18Pts[18][3] = {
+    {0, 0, 0}, // 0
+    {1, 0, 0}, // 1
+    {0, 1, 0}, // 2
+    {0, 0, 1},  // 3
+    {1, 0, 1},  // 4
+    {0, 1, 1},  // 5
+    {.5, 0, 0},  // 6
+    {0, .5, 0},  // 7
+    {0, 0, .5},  // 8
+    {.5, .5, 0},  // 9
+    {1, 0, .5},  // 10
+    {0, 1, .5},  // 11
+    {.5, 0, 1},  // 12
+    {0, .5, 1},  // 13
+    {.5, .5, 1},  // 14
+    {.5, 0, .5},  // 15
+    {0, .5, .5},  // 16
+    {.5, .5, .5},  // 17
+  };
+
+  int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+  fullMatrix<double> point(nbPoints, 3);
+
+  int index = 0;
+
+  if (order == 2)
+    for (int i =0; i<18; i++)
+      for (int j=0; j<3;j++)
+        point(i,j) = prism18Pts[i][j];
+  else {
+    fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
+    fullMatrix<double> linePoints = generate1DPoints(order);
+    for (int j = 0; j <linePoints.size1() ; j++) {
+      for (int i = 0; i < triPoints.size1(); i++) {
+        point(index,0) = triPoints(i,0);
+        point(index,1) = triPoints(i,1);
+        point(index,2) = linePoints(j,0);
+        index ++;
+      }
+    }
+  }
+  return point;
+}
+
+static fullMatrix<double> generatePointsHex(int order, bool serendip)
+{
+  int nbPoints = (order+1) * (order+1) * (order+1);
+  fullMatrix<double> point(nbPoints, 3);
+
+  fullMatrix<double> linePoints = generate1DPoints(order);
+  int index = 0;
+
+  for (int i = 0; i < linePoints.size1(); ++i) {
+    for (int j = 0; j < linePoints.size1(); ++j) {
+      for (int k = 0; k < linePoints.size1(); ++k) {
+        point(index, 0) = linePoints(i, 0);
+        point(index, 1) = linePoints(j, 0);
+        point(index, 2) = linePoints(k, 0);
+        ++index;
+      }
+    }
+  }
+
+  return point;
+}
+
+// Sub Control Points
+static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
+{
+  std::vector< fullMatrix<double> > subPoints(2);
+
+  subPoints[0] = generate1DPoints(order);
+  subPoints[0].scale(.5);
+
+  subPoints[1].copy(subPoints[0]);
+  subPoints[1].add(.5);
+
+  return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
+{
+  std::vector< fullMatrix<double> > subPoints(4);
+  fullMatrix<double> prox;
+
+  subPoints[0] = gmshGeneratePointsTriangle(order, false);
+  subPoints[0].scale(.5);
+
+  subPoints[1].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[1], 0, 1);
+  prox.add(.5);
+
+  subPoints[2].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[2], 1, 1);
+  prox.add(.5);
+
+  subPoints[3].copy(subPoints[0]);
+  subPoints[3].scale(-1.);
+  subPoints[3].add(.5);
+
+  return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
+{
+  std::vector< fullMatrix<double> > subPoints(4);
+  fullMatrix<double> prox;
+
+  subPoints[0] = generatePointsQuad(order, false);
+  subPoints[0].scale(.5);
+
+  subPoints[1].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[1], 0, 1);
+  prox.add(.5);
+
+  subPoints[2].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[2], 1, 1);
+  prox.add(.5);
+
+  subPoints[3].copy(subPoints[1]);
+  prox.setAsProxy(subPoints[3], 1, 1);
+  prox.add(.5);
+
+  return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
+{
+  std::vector< fullMatrix<double> > subPoints(8);
+  fullMatrix<double> prox1;
+  fullMatrix<double> prox2;
+
+  subPoints[0] = gmshGeneratePointsTetrahedron(order, false);
+  subPoints[0].scale(.5);
+
+  subPoints[1].copy(subPoints[0]);
+  prox1.setAsProxy(subPoints[1], 0, 1);
+  prox1.add(.5);
+
+  subPoints[2].copy(subPoints[0]);
+  prox1.setAsProxy(subPoints[2], 1, 1);
+  prox1.add(.5);
+
+  subPoints[3].copy(subPoints[0]);
+  prox1.setAsProxy(subPoints[3], 2, 1);
+  prox1.add(.5);
+
+  // u := .5-u-w
+  // v := .5-v-w
+  // w := w
+  subPoints[4].copy(subPoints[0]);
+  prox1.setAsProxy(subPoints[4], 0, 2);
+  prox1.scale(-1.);
+  prox1.add(.5);
+  prox1.setAsProxy(subPoints[4], 0, 1);
+  prox2.setAsProxy(subPoints[4], 2, 1);
+  prox1.add(prox2, -1.);
+  prox1.setAsProxy(subPoints[4], 1, 1);
+  prox1.add(prox2, -1.);
+
+  // u := u
+  // v := .5-v
+  // w := w+v
+  subPoints[5].copy(subPoints[0]);
+  prox1.setAsProxy(subPoints[5], 2, 1);
+  prox2.setAsProxy(subPoints[5], 1, 1);
+  prox1.add(prox2);
+  prox2.scale(-1.);
+  prox2.add(.5);
+
+  // u := .5-u
+  // v := v
+  // w := w+u
+  subPoints[6].copy(subPoints[0]);
+  prox1.setAsProxy(subPoints[6], 2, 1);
+  prox2.setAsProxy(subPoints[6], 0, 1);
+  prox1.add(prox2);
+  prox2.scale(-1.);
+  prox2.add(.5);
+
+  // u := u+w
+  // v := v+w
+  // w := .5-w
+  subPoints[7].copy(subPoints[0]);
+  prox1.setAsProxy(subPoints[7], 0, 1);
+  prox2.setAsProxy(subPoints[7], 2, 1);
+  prox1.add(prox2);
+  prox1.setAsProxy(subPoints[7], 1, 1);
+  prox1.add(prox2);
+  prox2.scale(-1.);
+  prox2.add(.5);
+
+
+  return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
+{
+  std::vector< fullMatrix<double> > subPoints(8);
+  fullMatrix<double> prox;
+
+  subPoints[0] = generatePointsPrism(order, false);
+  subPoints[0].scale(.5);
+
+  subPoints[1].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[1], 0, 1);
+  prox.add(.5);
+
+  subPoints[2].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[2], 1, 1);
+  prox.add(.5);
+
+  subPoints[3].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[3], 0, 2);
+  prox.scale(-1.);
+  prox.add(.5);
+
+  subPoints[4].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[4], 2, 1);
+  prox.add(.5);
+
+  subPoints[5].copy(subPoints[1]);
+  prox.setAsProxy(subPoints[5], 2, 1);
+  prox.add(.5);
+
+  subPoints[6].copy(subPoints[2]);
+  prox.setAsProxy(subPoints[6], 2, 1);
+  prox.add(.5);
+
+  subPoints[7].copy(subPoints[3]);
+  prox.setAsProxy(subPoints[7], 2, 1);
+  prox.add(.5);
+
+  return subPoints;
+}
+
+static std::vector< fullMatrix<double> > generateSubPointsHex(int order)
+{
+  std::vector< fullMatrix<double> > subPoints(8);
+  fullMatrix<double> prox;
+
+  subPoints[0] = generatePointsHex(order, false);
+  subPoints[0].scale(.5);
+
+  subPoints[1].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[1], 0, 1);
+  prox.add(.5);
+
+  subPoints[2].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[2], 1, 1);
+  prox.add(.5);
+
+  subPoints[3].copy(subPoints[1]);
+  prox.setAsProxy(subPoints[3], 1, 1);
+  prox.add(.5);
+
+  subPoints[4].copy(subPoints[0]);
+  prox.setAsProxy(subPoints[4], 2, 1);
+  prox.add(.5);
+
+  subPoints[5].copy(subPoints[1]);
+  prox.setAsProxy(subPoints[5], 2, 1);
+  prox.add(.5);
+
+  subPoints[6].copy(subPoints[2]);
+  prox.setAsProxy(subPoints[6], 2, 1);
+  prox.add(.5);
+
+  subPoints[7].copy(subPoints[3]);
+  prox.setAsProxy(subPoints[7], 2, 1);
+  prox.add(.5);
+
+  return subPoints;
+}
+
+// Matrices generation
+static int nChoosek(int n, int k)
+{
+  if (n < k || k < 0) {
+    Msg::Error("Wrong argument for combination.");
+    return 1;
+  }
+
+  if (k > n/2) k = n-k;
+  if (k == 1)
+    return n;
+  if (k == 0)
+    return 1;
+
+  int c = 1;
+  for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
+  return c;
+}
+
+
+static fullMatrix<double> generateBez2LagMatrix
+  (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
+   int order, int dimSimplex)
+{
+
+  if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
+    Msg::Fatal("Wrong sizes for Bezier coefficients generation %d %d -- %d %d",
+      exponent.size1(),point.size1(),
+      exponent.size2(),point.size2());
+    return fullMatrix<double>(1, 1);
+  }
+
+  int ndofs = exponent.size1();
+  int dim = exponent.size2();
+
+  fullMatrix<double> bez2Lag(ndofs, ndofs);
+  for (int i = 0; i < ndofs; i++) {
+    for (int j = 0; j < ndofs; j++) {
+      double dd = 1.;
+
+      {
+        double pointCompl = 1.;
+        int exponentCompl = order;
+        for (int k = 0; k < dimSimplex; k++) {
+          dd *= nChoosek(exponentCompl, (int) exponent(i, k))
+            * pow(point(j, k), exponent(i, k));
+          pointCompl -= point(j, k);
+          exponentCompl -= (int) exponent(i, k);
+        }
+        dd *= pow(pointCompl, exponentCompl);
+      }
+
+      for (int k = dimSimplex; k < dim; k++)
+        dd *= nChoosek(order, (int) exponent(i, k))
+            * pow(point(j, k), exponent(i, k))
+            * pow(1. - point(j, k), order - exponent(i, k));
+
+      bez2Lag(j, i) = dd;
+    }
+  }
+  return bez2Lag;
+}
+
+
+
+static fullMatrix<double> generateSubDivisor
+  (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
+   const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
+{
+  if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
+    Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+      exponents.size1(), lag2Bez.size1(),
+      exponents.size1(), lag2Bez.size2());
+    return fullMatrix<double>(1, 1);
+  }
+
+  int nbPts = lag2Bez.size1();
+  int nbSubPts = nbPts * subPoints.size();
+
+  fullMatrix<double> intermediate2(nbPts, nbPts);
+  fullMatrix<double> subDivisor(nbSubPts, nbPts);
+
+  for (unsigned int i = 0; i < subPoints.size(); i++) {
+    fullMatrix<double> intermediate1 =
+      generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
+    lag2Bez.mult(intermediate1, intermediate2);
+    subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
+  }
+  return subDivisor;
+}
+
+
+
+// Convert Bezier points to Lagrange points
+static fullMatrix<double> bez2LagPoints(bool dim1, bool dim2, bool dim3, const fullMatrix<double> &bezierPoints)
+{
+
+  // FIXME BUG for 2nd order quads we seem to try to access bezierPoints(i, 2);
+  // but bezierPoints only has 2 columns
+
+  const int numPt = bezierPoints.size1();
+  fullMatrix<double> lagPoints(numPt,3);
+
+  for (int i=0; i<numPt; i++) {
+    lagPoints(i,0) = dim1 ? -1. + 2*bezierPoints(i,0) : bezierPoints(i,0);
+    lagPoints(i,1) = dim2 ? -1. + 2*bezierPoints(i,1) : bezierPoints(i,1);
+    lagPoints(i,2) = dim3 ? -1. + 2*bezierPoints(i,2) : bezierPoints(i,2);
+  }
+
+  return lagPoints;
+
+}
+
+
+bezierBasis::bezierBasis(int tag)
+{
+  order = MElement::OrderFromTag(tag);
+
+  if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
+    dim = 3;
+    numLagPts = 5;
+    bezierPoints = gmshGeneratePointsPyramid(order,false);
+    lagPoints = bezierPoints;
+    matrixLag2Bez.resize(bezierPoints.size1(), bezierPoints.size1(), false);
+    matrixBez2Lag.resize(bezierPoints.size1(), bezierPoints.size1(), false);
+    for (int i=0; i<matrixBez2Lag.size1(); ++i) {
+      matrixLag2Bez(i,i) = 1.;
+      matrixBez2Lag(i,i) = 1.;
+    }
+    // TODO: Implement subdidivsor
+  }
+  else {
+    int dimSimplex;
+    std::vector< fullMatrix<double> > subPoints;
+    switch (MElement::ParentTypeFromTag(tag)) {
+      case TYPE_PNT :
+        dim = 0;
+        numLagPts = 1;
+        exponents = generate1DExponents(0);
+        bezierPoints = generate1DPoints(0);
+        lagPoints = bezierPoints;
+        dimSimplex = 0;
+        break;
+      case TYPE_LIN : {
+        dim = 1;
+        numLagPts = 2;
+        exponents = generate1DExponents(order);
+        bezierPoints = generate1DPoints(order);
+        lagPoints = bez2LagPoints(true,false,false,bezierPoints);
+        dimSimplex = 0;
+        subPoints = generateSubPointsLine(0);
+        break;
+      }
+      case TYPE_TRI : {
+        dim = 2;
+        numLagPts = 3;
+        exponents = generateExponentsTriangle(order);
+        bezierPoints = gmshGeneratePointsTriangle(order,false);
+        lagPoints = bezierPoints;
+        dimSimplex = 2;
+        subPoints = generateSubPointsTriangle(order);
+        break;
+      }
+      case TYPE_QUA : {
+        dim = 2;
+        numLagPts = 4;
+        exponents = generateExponentsQuad(order);
+        bezierPoints = generatePointsQuad(order,false);
+        lagPoints = bez2LagPoints(true,true,false,bezierPoints);
+        dimSimplex = 0;
+        subPoints = generateSubPointsQuad(order);
+        //      points.print("points");
+        break;
+      }
+      case TYPE_TET : {
+        dim = 3;
+        numLagPts = 4;
+        exponents = generateExponentsTetrahedron(order);
+        bezierPoints = gmshGeneratePointsTetrahedron(order,false);
+        lagPoints = bezierPoints;
+        dimSimplex = 3;
+        subPoints = generateSubPointsTetrahedron(order);
+        break;
+      }
+      case TYPE_PRI : {
+        dim = 3;
+        numLagPts = 6;
+        exponents = generateExponentsPrism(order);
+        bezierPoints = generatePointsPrism(order, false);
+        lagPoints = bez2LagPoints(false,false,true,bezierPoints);
+        dimSimplex = 2;
+        subPoints = generateSubPointsPrism(order);
+        break;
+      }
+      case TYPE_HEX : {
+        dim = 3;
+        numLagPts = 8;
+        exponents = generateExponentsHex(order);
+        bezierPoints = generatePointsHex(order, false);
+        lagPoints = bez2LagPoints(true,true,true,bezierPoints);
+        dimSimplex = 0;
+        subPoints = generateSubPointsHex(order);
+        break;
+      }
+      default : {
+        Msg::Error("Unknown function space %d: reverting to TET_1", tag);
+        dim = 3;
+        numLagPts = 4;
+        exponents = generateExponentsTetrahedron(0);
+        bezierPoints = gmshGeneratePointsTetrahedron(0, false);
+        lagPoints = bezierPoints;
+        dimSimplex = 3;
+        subPoints = generateSubPointsTetrahedron(0);
+        break;
+      }
+    }
+    Msg::Info("%d %d %d %d %d %d", exponents.size1(), exponents.size2(), bezierPoints.size1(), bezierPoints.size2(), order, dimSimplex);
+    matrixBez2Lag = generateBez2LagMatrix(exponents, bezierPoints, order, dimSimplex);
+    Msg::Info("%d %d %d %d", matrixBez2Lag.size1(), matrixBez2Lag.size2(), matrixLag2Bez.size1(), matrixLag2Bez.size2());
+    Msg::Info("%f %f %f %f", bezierPoints(0, 0), bezierPoints(0, 1), exponents(0, 0), exponents(0, 1));
+    matrixBez2Lag.invert(matrixLag2Bez);
+    subDivisor = generateSubDivisor(exponents, subPoints, matrixLag2Bez, order, dimSimplex);
+    numDivisions = (int) pow(2., (int) bezierPoints.size2());
+  }
+}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..e8241ff6b7881c360c2d4f6a6bb519b731873876
--- /dev/null
+++ b/Numeric/bezierBasis.h
@@ -0,0 +1,30 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef _BEZIER_BASIS_H_
+#define _BEZIER_BASIS_H_
+
+#include <map>
+#include <vector>
+#include "fullMatrix.h"
+
+class bezierBasis {
+ private :
+ public :
+  int dim, order;
+  int numLagPts;
+  int numDivisions;
+  // the 'numLagPts' first exponents are related to 'Lagrangian' values
+  bezierBasis() {Msg::Fatal("Should not be created this way");};
+  bezierBasis(int tag);
+  fullMatrix<double> exponents; // exponents of Bezier FF
+  fullMatrix<double> bezierPoints, lagPoints; // sampling point
+  fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
+  fullMatrix<double> subDivisor;
+};
+
+#endif
+
+
diff --git a/Numeric/jacobiPolynomials.h b/Numeric/jacobiPolynomials.h
index 9d83f5077f16886f3f7d8c9412cdbffd5d0daa91..b3d510c3412a85fe05617b90eddc2e4860eb5539 100644
--- a/Numeric/jacobiPolynomials.h
+++ b/Numeric/jacobiPolynomials.h
@@ -11,7 +11,7 @@
 class JacobiPolynomials {
 
  public:
-  JacobiPolynomials() {};
+  JacobiPolynomials() : n(0) {};
   JacobiPolynomials(double a, double b, int o);
   ~JacobiPolynomials();
 
diff --git a/Numeric/legendrePolynomials.h b/Numeric/legendrePolynomials.h
index f7fd118ddf0ad5fd037ba139adc182f751fdf49d..381e62e24964d9d3173276da0311964452979e25 100644
--- a/Numeric/legendrePolynomials.h
+++ b/Numeric/legendrePolynomials.h
@@ -11,7 +11,7 @@
 class LegendrePolynomials {
 
  public:
-  LegendrePolynomials() {};
+  LegendrePolynomials() : n(0) {};
   LegendrePolynomials(int o);
   ~LegendrePolynomials();
 
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 81675224b72a588548824ce3d46c8cb596916904..38126f89c91b09b77e0e73ac68f6c7a2cbcdb7ea 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -1141,7 +1141,7 @@ static void generateFaceClosureHex(nodalBasis::clCont &closure, int order,
                                    bool serendip, const fullMatrix<double> &points)
 {
   closure.clear();
-  const nodalBasis &fsFace = *BasisFactory::create(nodalBasis::getTag(TYPE_QUA, order, serendip));
+  const nodalBasis &fsFace = *BasisFactory::getNodalBasis(nodalBasis::getTag(TYPE_QUA, order, serendip));
   for (int iRotate = 0; iRotate < 4; iRotate++){
     for (int iSign = 1; iSign >= -1; iSign -= 2){
       for (int iFace = 0; iFace < 6; iFace++) {
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 4516507d35b6ad11c9f1c3ca450f63ef55bbf9d5..0c09ef633c896f699efc23abe42a2f98e5fc49e6 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -30,8 +30,7 @@ inline double pow_int(const double &a, const int &n)
   case 5 :
     {
       const double a2 = a*a;
-      const double a3 = a*a*a;
-      return a2*a3;
+      return a2*a2*a;
     }
   case 6 :
     {
@@ -51,9 +50,8 @@ inline double pow_int(const double &a, const int &n)
     }
   case 9 :
     {
-      const double a2 = a*a;
-      const double a4 = a2*a2;
-      return a4*a4*a;
+      const double a3 = a*a*a;
+      return a3*a3*a3;
     }
   case 10 :
     {
@@ -93,7 +91,7 @@ class polynomialBasis : public nodalBasis
   virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
   virtual void df(double u, double v, double w, double grads[][3]) const;
   virtual void ddf(double u, double v, double w, double hess[][3][3]) const;
-  virtual  void dddf(double u, double v, double w, double third[][3][3][3]) const;
+  virtual void dddf(double u, double v, double w, double third[][3][3][3]) const;
 
   virtual int getNumShapeFunctions() const;
 
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 05e5653bde94f347bd9309de974aef7b76530fe5..1e4f5438ce8d92dc3a847a0419e1060c50fc2704 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -6,7 +6,7 @@
 #include "Gmsh.h"
 #include "GModel.h"
 #include "GmshMessage.h"
-#include "JacobianBasis.h"
+
 #include "polynomialBasis.h"
 #include "AnalyseCurvedMesh.h"
 #include "Context.h"
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 01feaae9796ca853c04ce2528aea3c6af0aea81a..211bb224677dda8b731d617eef662d62fc3af012 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -903,7 +903,7 @@ void PViewDataList::setOrder2(int type)
   case TYPE_PRI: typeMSH = MSH_PRI_18; break;
   // case TYPE_PYR: typeMSH = MSH_PYR_14; break;
   }
-  const polynomialBasis *fs = (polynomialBasis*)BasisFactory::create(typeMSH);
+  const polynomialBasis *fs = (polynomialBasis*)BasisFactory::getNodalBasis(typeMSH);
   if(!fs){
     Msg::Error("Could not find polynomial function space for element type %d",
                typeMSH);
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 478fd0cace1de5f701bee190b07ffcaa56c202a0..3d2dbaba46031d2323053d8acb03b53f92c68602 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -1045,7 +1045,7 @@ bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2)
 const nodalBasis* DI_Line::getFunctionSpace(int o) const{
   int order = (o == -1) ? getPolynomialOrder() : o;
   int tag = polynomialBasis::getTag(TYPE_LIN, order);
-  return BasisFactory::create(tag);
+  return BasisFactory::getNodalBasis(tag);
 }
 
 void DI_Line::computeIntegral() {
@@ -1066,7 +1066,7 @@ const nodalBasis* DI_Triangle::getFunctionSpace(int o) const
 {
   int order = (o == -1) ? getPolynomialOrder() : o;
   int tag = polynomialBasis::getTag(TYPE_TRI, order);
-  return BasisFactory::create(tag);
+  return BasisFactory::getNodalBasis(tag);
 }
 void DI_Triangle::computeIntegral() {
   integral_ = TriSurf(pt(0), pt(1), pt(2));
@@ -1090,7 +1090,7 @@ double DI_Triangle::quality() const {
 const nodalBasis* DI_Quad::getFunctionSpace(int o) const{
  int order = (o == -1) ? getPolynomialOrder() : o;
   int tag = polynomialBasis::getTag(TYPE_QUA, order);
-  return BasisFactory::create(tag);
+  return BasisFactory::getNodalBasis(tag);
 }
 
 void DI_Quad::computeIntegral() {
@@ -1114,7 +1114,7 @@ void DI_Quad::computeIntegral() {
 const nodalBasis* DI_Tetra::getFunctionSpace(int o) const{
  int order = (o == -1) ? getPolynomialOrder() : o;
   int tag = polynomialBasis::getTag(TYPE_TET, order);
-  return BasisFactory::create(tag);
+  return BasisFactory::getNodalBasis(tag);
 }
 
 void DI_Tetra::computeIntegral() {
@@ -1129,7 +1129,7 @@ double DI_Tetra::quality() const {
 const nodalBasis* DI_Hexa::getFunctionSpace(int o) const{
   int order = (o == -1) ? getPolynomialOrder() : o;
   int tag = polynomialBasis::getTag(TYPE_HEX, order);
-  return BasisFactory::create(tag);
+  return BasisFactory::getNodalBasis(tag);
 }
 void DI_Hexa::computeIntegral() {
     integral_ = TetraVol(pt(0), pt(1), pt(3), pt(4)) + TetraVol(pt(1), pt(4), pt(5), pt(7))
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index c791cfc1fba67c07781382aed74dd52d2e0e4ed4..132405c539969aac7b63a146314f31d1ec55f006 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -49,7 +49,7 @@ class DI_Point
             const std::vector<gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
   virtual ~DI_Point(){}
   virtual const nodalBasis* getFunctionSpace(int o) const
-  { return BasisFactory::create(MSH_PNT); }
+  { return BasisFactory::getNodalBasis(MSH_PNT); }
   virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
   {
     s[0] = 1.;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 481d6a2b8e292c1e0268626dea52378135efbd44..e8f526c1658b04793cf9376732970c4bad5741e8 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -8,7 +8,7 @@
 #include "MVertex.h"
 #include "ParamCoord.h"
 #include "polynomialBasis.h"
-#include "JacobianBasis.h"
+