diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 7a3623ffabb9f4da3a97a31af631d30dd08f43a9..7bfdc27caf9f9e51539bca5eccc4bd2775c7c0aa 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -165,7 +165,6 @@ class MPyramid : public MElement {
   }
   static int faces_pyramid(const int face, const int vert)
   {
-    // only triangular faces
     static const int f[5][4] = {
       {0, 1, 4, -1},
       {3, 0, 4, -1},
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index af719e909c81e907efc9edebac1338bdd4e648e7..20af53c78f1f85203916b7a58c3c8eb028e2d4a6 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -27,8 +27,8 @@ JacobianBasis::JacobianBasis(int tag)
 {
   const int parentType = ElementType::ParentTypeFromTag(tag);
   const int order = ElementType::OrderFromTag(tag);
-  const int jacobianOrder =  JacobianBasis::jacobianOrder(parentType, order);
-  const int primJacobianOrder =  JacobianBasis::jacobianOrder(parentType, 1);
+  const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+  const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
 
   fullMatrix<double> lagPoints;                                  // Sampling points
 
@@ -728,7 +728,8 @@ fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
   return points;
 }
 
-int JacobianBasis::jacobianOrder(int parentType, int order) {
+int JacobianBasis::jacobianOrder(int parentType, int order)
+{
   switch (parentType) {
     case TYPE_PNT : return 0;
     case TYPE_LIN : return order - 1;
@@ -745,3 +746,14 @@ int JacobianBasis::jacobianOrder(int parentType, int order) {
       return 0;
   }
 }
+
+
+void JacobianBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
+                                          fullMatrix<double> *dxyzdX,
+                                          fullMatrix<double> *dxyzdY,
+                                          fullMatrix<double> *dxyzdZ) const
+{
+  if (dxyzdX) gradShapeMatX.mult(nodes, *dxyzdX);
+  if (dxyzdY) gradShapeMatY.mult(nodes, *dxyzdY);
+  if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index f1655c2e16f9520c5260bcb17e1048c487be3c57..8a08f18e69ef38e11b2d73acaca3d1b8cfff4b95 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -14,8 +14,6 @@
 
 class JacobianBasis {
  private:
-  const bezierBasis *bezier;
-
   fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
   fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
   fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
@@ -41,6 +39,8 @@ class JacobianBasis {
                                        fullMatrix<double> &JDJ) const;
 
  public :
+  const bezierBasis *bezier;
+
   JacobianBasis(int tag);
 
   // Get methods
@@ -106,6 +106,12 @@ class JacobianBasis {
   }
 
   // Jacobian basis order and pyramidal basis
+  void getGradientsFromNodes(const fullMatrix<double> &nodes,
+                             fullMatrix<double> *dxyzdX,
+                             fullMatrix<double> *dxyzdY,
+                             fullMatrix<double> *dxyzdZ) const;
+
+  //
   static int jacobianOrder(int parentType, int order);
   static fullMatrix<double> generateJacMonomialsPyramid(int order);
   static fullMatrix<double> generateJacPointsPyramid(int order);
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index ae93f045880cef9b80962b95cea025b13800f907..054fdaa958ddff38f2f7265284e44e1dd11ae5fe 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -247,6 +247,10 @@ namespace {
       return subPoints;
     }
 
+    Msg::Error("Subpoints for pyramid p>1 not implemented !");
+    return std::vector< fullMatrix<double> >();
+
+    /*
     std::vector< fullMatrix<double> > subPoints(8);
     fullMatrix<double> ref, prox;
 
@@ -296,6 +300,7 @@ namespace {
     }
 
     return subPoints;
+    */
   }
 
   // Matrices generation
@@ -378,18 +383,18 @@ namespace {
       for (int j = 0; j < ndofs; j++) {
         bez2Lag(i, j) =
             nChoosek(order, exponent(j, 2))
-            * pow(point(i, 2), exponent(j, 2)) //FIXME use pow_int
-            * pow(1. - point(i, 2), order - exponent(j, 2));
-
+            * pow_int(point(i, 2), exponent(j, 2))
+            * pow_int(1. - point(i, 2), order - exponent(j, 2));
+  
         double p = order + 2 - exponent(j, 2);
         double denom = 1. - point(i, 2);
         bez2Lag(i, j) *=
               nChoosek(p, exponent(j, 0))
             * nChoosek(p, exponent(j, 1))
-            * pow(point(i, 0) / denom, exponent(j, 0))
-            * pow(point(i, 1) / denom, exponent(j, 1))
-            * pow(1. - point(i, 0) / denom, p - exponent(j, 0))
-            * pow(1. - point(i, 1) / denom, p - exponent(j, 1));
+            * pow_int(point(i, 0) / denom, exponent(j, 0))
+            * pow_int(point(i, 1) / denom, exponent(j, 1))
+            * pow_int(1. - point(i, 0) / denom, p - exponent(j, 0))
+            * pow_int(1. - point(i, 1) / denom, p - exponent(j, 1));
       }
     }
     return bez2Lag;
@@ -458,16 +463,23 @@ void bezierBasis::_construct(int parentType, int p)
     dim = 3;
     numLagCoeff = 8;
     exponents = JacobianBasis::generateJacMonomialsPyramid(order);
-    subPoints = generateSubPointsPyr(order);
-
-    numDivisions = static_cast<int>(subPoints.size());
+    if (order < 2) {
+      subPoints = generateSubPointsPyr(order);
+      numDivisions = static_cast<int>(subPoints.size());
+    }
+    else {
+      numDivisions = 8;
+      static int a = 0;
+      if (++a == 1) Msg::Warning("Subdivision not available for pyramid p>1");
+    }
 
     fullMatrix<double> bezierPoints = exponents;
     bezierPoints.scale(1. / (order + 2));
 
     matrixBez2Lag = generateBez2LagMatrixPyramid(exponents, bezierPoints, order);
     matrixBez2Lag.invert(matrixLag2Bez);
-    subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
+    if (order < 2)
+      subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
     return;
   }