diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 383a213684f25dd41036d695349ef2ba589e3562..86fce5eec93587866eef0dc19042151dc44279eb 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -22,6 +22,7 @@ set(SRC
     GaussQuadratureTet.cpp
     GaussQuadratureHex.cpp
     GaussQuadraturePri.cpp
+    GaussQuadraturePyr.cpp
     GaussLegendreSimplex.cpp
   GaussJacobi1D.cpp
   robustPredicates.cpp
diff --git a/Numeric/GaussIntegration.cpp b/Numeric/GaussIntegration.cpp
index 756ca222d2f73a9c3796cd4674feffbebd1a6034..69be54233b4b6b81ebbb854114430873a046e796 100644
--- a/Numeric/GaussIntegration.cpp
+++ b/Numeric/GaussIntegration.cpp
@@ -55,6 +55,12 @@ void gaussIntegration::getPrism(int order, fullMatrix<double> &pts,
   pts2fullMatrix(getNGQPriPts(order),getGQPriPts(order),pts,weights);
 }
 
+void gaussIntegration::getPyramid(int order, fullMatrix<double> &pts,
+                                fullVector<double> &weights)
+{
+  pts2fullMatrix(getNGQPyrPts(order),getGQPyrPts(order),pts,weights);
+}
+
 void gaussIntegration::get(int elementType, int order, fullMatrix<double> &pts,
                            fullVector<double> &weights)
 {
@@ -65,6 +71,7 @@ void gaussIntegration::get(int elementType, int order, fullMatrix<double> &pts,
   case TYPE_TET : pts2fullMatrix(getNGQTetPts(order), getGQTetPts(order), pts, weights); break;
   case TYPE_HEX : pts2fullMatrix(getNGQHPts(order), getGQHPts(order), pts, weights); break;
   case TYPE_PRI : pts2fullMatrix(getNGQPriPts(order), getGQPriPts(order), pts, weights); break;
+  case TYPE_PYR : pts2fullMatrix(getNGQPyrPts(order), getGQPyrPts(order), pts, weights); break;
   case TYPE_PNT :
     weights.resize(1,1);
     weights(0) = 1.;
diff --git a/Numeric/GaussIntegration.h b/Numeric/GaussIntegration.h
index 9c9079965a42016e81752c577968e3a2a4df9f2a..54f6f06bfa77959e957ef341d38bd231008c8268 100644
--- a/Numeric/GaussIntegration.h
+++ b/Numeric/GaussIntegration.h
@@ -31,6 +31,9 @@ IntPt *getGQTetPts(int order);
 int getNGQPriPts(int order);
 IntPt *getGQPriPts(int order);
 
+int getNGQPyrPts(int order);
+IntPt *getGQPyrPts(int order);
+
 int getNGQHPts(int order);
 IntPt *getGQHPts(int order);
 
@@ -50,6 +53,8 @@ class gaussIntegration {
                             fullVector<double> &weights);
   static void getPrism(int order, fullMatrix<double> &pts,
                        fullVector<double> &weights);
+  static void getPyramid(int order, fullMatrix<double> &pts,
+                       fullVector<double> &weights);
 };
 
 #endif
diff --git a/Numeric/GaussQuadraturePyr.cpp b/Numeric/GaussQuadraturePyr.cpp
index 35b2d7b83cb066d5f2b8f08460e8f61d26f1cfab..30c4bda94b476496f230bc4a600b329844c96189 100644
--- a/Numeric/GaussQuadraturePyr.cpp
+++ b/Numeric/GaussQuadraturePyr.cpp
@@ -6,37 +6,65 @@
 #include "GmshMessage.h"
 #include "GaussIntegration.h"
 #include "GaussLegendre1D.h"
+#include "GaussJacobi1D.h"
 
 IntPt *getGQPyrPts(int order);
 int getNGQPyrPts(int order);
 
-IntPt * GQP[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+IntPt * GQPyr[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 
 IntPt *getGQPyrPts(int order)
 {
-  int nLin = (order+3)/2;
-  int nTri = getNGQTPts(order);
-  int n = nLin*nTri;
+
   int index = order;
-  if (order >= (int)(sizeof(GQP) / sizeof(IntPt*)))
-    Msg::Fatal("Increase size of GQP in gauss quadrature prism");
-  if(!GQP[index]){
+
+  if(!GQPyr[index]) {
+
+    int nbPtUV = order/2 + 1;
+    int nbPtW  = order/2 + 1;
+    int nbPtUV2 = nbPtUV * nbPtUV;
+
     double *linPt,*linWt;
-    IntPt *triPts = getGQTPts(order);
-    gmshGaussLegendre1D(nLin,&linPt,&linWt);
-    GQP[index] = new IntPt[n];
+    gmshGaussLegendre1D(nbPtUV,&linPt,&linWt);
+
+    double *GJ20Pt, *GJ20Wt;
+    getGaussJacobiQuadrature(2, 0, nbPtW, &GJ20Pt, &GJ20Wt);
+
+    GQPyr[index] = new IntPt[getNGQPyrPts(order)];
+    if (order >= (int)(sizeof(GQPyr) / sizeof(IntPt*)))
+      Msg::Fatal("Increase size of GQPyr in gauss quadrature prism");
+
     int l = 0;
-    for(int i=0; i < nTri; i++) {
-      for(int j=0; j < nLin; j++) {
-        GQP[index][l].pt[0] = triPts[i].pt[0];
-        GQP[index][l].pt[1] = triPts[i].pt[1];
-        GQP[index][l].pt[2] = linPt[j];
-        GQP[index][l++].weight = triPts[i].weight*linWt[j];
-      }
+    for (int i = 0; i < getNGQPyrPts(order); i++) {
+
+      // compose an integration rule for (1-w)^2 f(u,v,w) on the standard hexahedron
+
+      int iW = i / (nbPtUV2);
+      int iU = (i - iW*nbPtUV2)/nbPtUV;
+      int iV = (i - iW*nbPtUV2 - iU*nbPtUV);
+
+      // std::cout << "Points " << iU << " " << iV << " " << iW << std::endl;
+
+      int wt = linWt[iU]*linWt[iV]*GJ20Wt[iW];
+
+      double up = linPt[iU];
+      double vp = linPt[iV];
+      double wp = GJ20Pt[iW];
+      
+      // now incorporate the Duffy transformation from pyramid to hexahedron
+
+      GQPyr[index][l].pt[0] = 0.5*(1-wp)*up;
+      GQPyr[index][l].pt[1] = 0.5*(1-wp)*vp;
+      GQPyr[index][l].pt[2] = 0.5*(1+wp);
+
+      wt *= 0.125;
+      GQPyr[index][l++].weight = wt;
+
     }
+
   }
-  return GQP[index];
+  return GQPyr[index];
 }