Added Gauss quadrature points and weights for pyramids

Numeric/CMakeLists.txt

@@ -22,6 +22,7 @@ set(SRC
     GaussQuadratureTet.cpp
     GaussQuadratureHex.cpp
     GaussQuadraturePri.cpp
+    GaussQuadraturePyr.cpp
     GaussLegendreSimplex.cpp
   GaussJacobi1D.cpp
   robustPredicates.cpp

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.;

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

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];
 }