diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 0d1b8683770c61a595358df320babc1b9071ac65..05ac767ff314e18ae7352b9473f3128e7b03d707 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -54,7 +54,7 @@ GradientBasis::GradientBasis(int tag, int order)
       samplingPoints = gmshGeneratePointsHexahedron(order,false);
       break;
     case TYPE_PYR :
-      samplingPoints = JacobianBasis::generateJacPointsPyramid(order);
+      samplingPoints = gmshGeneratePointsPyramidGeneral(false, order+2, order);
       break;
     default :
       Msg::Error("Unknown Jacobian function space for element tag %d", tag);
@@ -92,10 +92,10 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
 
 JacobianBasis::JacobianBasis(int tag, int jacOrder) :
     _bezier(NULL), _tag(tag), _dim(ElementType::DimensionFromTag(tag)),
-    _jacOrder(jacOrder >= 0 ? jacOrder : JacobianBasis::jacobianOrder(tag))
+    _jacOrder(jacOrder >= 0 ? jacOrder : jacobianOrder(tag))
 {
   const int parentType = ElementType::ParentTypeFromTag(tag);
-  const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
+  const int primJacobianOrder = jacobianOrder(parentType, 1);
 
   fullMatrix<double> lagPoints;                                  // Sampling points
 
@@ -122,7 +122,7 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder) :
       lagPoints = gmshGeneratePointsHexahedron(_jacOrder,false);
       break;
     case TYPE_PYR :
-      lagPoints = generateJacPointsPyramid(_jacOrder);
+      lagPoints = gmshGeneratePointsPyramidGeneral(false, _jacOrder+2, _jacOrder);
       break;
     default :
       Msg::Error("Unknown Jacobian function space for element tag %d", tag);
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 92a25a4116ef8fd0734e6fe60a167b11acb224e9..fc8a48f060d21d8cca7e97cb7e126761e4472ffa 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -252,6 +252,9 @@ std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
     fullMatrix<double> ref, prox;
 
     subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, order+2, order);
+    prox.setAsProxy(subPoints[0], 2, 1);
+    prox.scale(-1);
+    prox.add(order);
     subPoints[0].scale(.5/(order+2));
 
     subPoints[1].copy(subPoints[0]);
@@ -380,10 +383,10 @@ fullMatrix<double> generateBez2LagMatrixPyramid
           * nChoosek(order, exponent(j, 2))
           * pow_int(point(i, 0) / denom, exponent(j, 0))
           * pow_int(point(i, 1) / denom, exponent(j, 1))
-          * pow_int(point(i, 2)        , exponent(j, 2))
+          * pow_int(1. - point(i, 2)   , exponent(j, 2))
           * pow_int(1. - point(i, 0) / denom, order + 2 - exponent(j, 0))
           * pow_int(1. - point(i, 1) / denom, order + 2 - exponent(j, 1))
-          * pow_int(1. - point(i, 2)        , order     - exponent(j, 2));
+          * pow_int(point(i, 2)             , order     - exponent(j, 2));
     }
   }
   return bez2Lag;
@@ -527,7 +530,7 @@ void bezierBasis::_construct(int parentType, int p)
     bezierPoints.resize(_exponents.size1(), _exponents.size2());
     const double ord = _order + 2;
     for (int i = 0; i < bezierPoints.size1(); ++i) {
-      bezierPoints(i, 2) = _exponents(i, 2) / ord;
+      bezierPoints(i, 2) = (_order - _exponents(i, 2)) / ord;
       const double scale = 1. - bezierPoints(i, 2);
       bezierPoints(i, 0) = _exponents(i, 0) / ord * scale;
       bezierPoints(i, 1) = _exponents(i, 1) / ord * scale;
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 485bf985756e970c14ffd47fe1f0186cb123cd86..69ac1fed18b2ec8760ebec39bd71ee2c019ef50b 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -107,9 +107,9 @@ fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, b
   for (int i = 0; i < points.size1(); ++i) {
     int div = pyr ? nk+nij : std::max(nij, nk);
     points(i, 2) = (nk - points(i, 2)) / div;
-    const double duv = -1. + points(i, 2);
-    points(i, 0) = duv + points(i, 0) * 2. / div;
-    points(i, 1) = duv + points(i, 1) * 2. / div;
+    const double scale = 1. - points(i, 2);
+    points(i, 0) = scale * (-1 + points(i, 0) * 2. / div);
+    points(i, 1) = scale * (-1 + points(i, 1) * 2. / div);
   }
   return points;
 }
@@ -833,6 +833,8 @@ fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(
   monomials(3, 1) = nijBase;
   monomials(3, 2) = nk;
 
+  int index = 4;
+
   if (nk > 0) {
     monomials(4, 0) = 0;
     monomials(4, 1) = 0;
@@ -849,9 +851,9 @@ fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(
     monomials(7, 0) = 0;
     monomials(7, 1) = nij;
     monomials(7, 2) = 0;
-  }
 
-  int index = 8;
+    index = 8;
+  }
 
   // Base
   if (nijBase > 1) {
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index c3f88f39abf8da08b1c0de4711541973b6296187..2d581f2f9dde0c3585edea68c3c136fb9845a8b0 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -28,10 +28,6 @@ fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip = false
 fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip = false);
 
 fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip = false);
-
-// If 'serendip' == true, generate points for serendipity pyramid of order 'nk',
-// else if 'pyr' == true, generate points for space {X^i Y^j Z^k | i,j <= k+'nij', k <= 'nk'},
-// else if 'pyr' == false, generate points for space {X^i Y^j Z^k | i,j <= 'nij', k <= 'nk'}
 fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, bool serendip = false);
 
 
@@ -49,9 +45,11 @@ fullMatrix<double> gmshGenerateMonomialsHexaSerendipity(int order);
 fullMatrix<double> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
 fullMatrix<double> gmshGenerateMonomialsPrismSerendipity(int order);
 
+// Generate monomials of pyramidal nodal space {X^i Y^j Z^k | i,j <= k, k <= 'order'},
 fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
-
-// See explanations for 'gmshGeneratePointsPyramidGeneral'
+// If 'serendip' == true, generate monomials of serendipity pyramid at order 'nk',
+// else if 'pyr' == true, generate monomials of space {X^i Y^j Z^k | i,j <= k+'nij', k <= 'nk'},
+// else if 'pyr' == false, generate monomials of space {X^i Y^j Z^k | i,j <= 'nij', k <= 'nk'}
 fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(bool pyr, int nij, int nk, bool forSerendipPoints = false);