A

Mesh/HighOrder.cpp

@@ -516,7 +516,7 @@ static int getNewFacePointsInVolume(MElement *incomplete, int nPts,
   case TYPE_PYR:
     switch (nPts){
     case 0:
-    case 1: return -1;
+    case 1: points = BasisFactory::getNodalBasis(MSH_PYR_14)->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;
@@ -1224,10 +1224,10 @@ static MPyramid *setHighOrder(MPyramid *p, GRegion *gr,
 {
   std::vector<MVertex*> ve, vf, vr;
   getEdgeVertices(gr, p, ve, newHOVert, edgeVertices, linear, nPts);
-//  MPyramidN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
-//                  p->getVertex(3), p->getVertex(4), ve, nPts + 1, 0, p->getPartition());
-//  getFaceVertices(gr, &incpl, p, vf, faceVertices, edgeVertices, linear, nPts);
-  getFaceVertices(gr, 0, p, vf, newHOVert, faceVertices, edgeVertices, linear, nPts);
+  MPyramidN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+                  p->getVertex(3), p->getVertex(4), ve, nPts + 1, 0, p->getPartition());
+  getFaceVertices(gr, &incpl, p, vf, newHOVert, faceVertices, edgeVertices, linear, nPts);
+  // getFaceVertices(gr, 0, p, vf, newHOVert, faceVertices, edgeVertices, linear, nPts);
   ve.insert(ve.end(), vf.begin(), vf.end());
   getVolumeVerticesPyramid(gr, p, ve, vr, newHOVert, linear, nPts);
   ve.insert(ve.end(), vr.begin(), vr.end());

Numeric/BergotBasis.cpp

@@ -9,8 +9,12 @@
 
 
 
-BergotBasis::BergotBasis(int p): order(p)
+BergotBasis::BergotBasis(int p,bool incpl): order(p),incomplete(incpl)
 {
+  
+  if (incomplete && order > 2) {
+    Msg::Error("Incomplete pyramids of order %i not yet implemented",order);
+  }
 }
 
 
@@ -20,6 +24,13 @@ BergotBasis::~BergotBasis()
 }
 
 
+bool BergotBasis::validIJ(int i,int j) const {
+
+  if (!incomplete)    return (i<=order) && (j<=order);
+  if (i+j <= order)   return true;
+  if (i+j == order+1) return i==1 || j==1;
+  return false;
+}
 
 // Values of Bergot basis functions for coordinates (u,v,w) in the unit pyramid:
 // f = L_i(uhat)*L_j(vhat)*(1-w)^max(i,j)*P_k^{2*max(i,j)+2,0}(what)
@@ -55,13 +66,14 @@ void BergotBasis::f(double u, double v, double w, double* val) const
   int index = 0;
   for (int i=0;i<=order;i++) {
     for (int j=0;j<=order;j++) {
-      int mIJ = std::max(i,j);
-      double fact = pow(1.-w, mIJ);
-      std::vector<double> &wf = wFcts[mIJ];
-      for (int k=0; k<=order-mIJ; k++,index++) val[index] = uFcts[i] * vFcts[j] * wf[k] * fact;
+      if (validIJ(i,j)) {
+        int mIJ = std::max(i,j);
+        double fact = pow(1.-w, mIJ);
+        std::vector<double> &wf = wFcts[mIJ];
+        for (int k=0; k<=order-mIJ; k++,index++) val[index] = uFcts[i] * vFcts[j] * wf[k] * fact;
+      }
     }
   }
-
 }
 
 
@@ -107,50 +119,52 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const
   int index = 0;
   for (int i=0;i<=order;i++) {
     for (int j=0;j<=order;j++) {
-      int mIJ = std::max(i,j);
-      std::vector<double> &wf = wFcts[mIJ], &wg = wGrads[mIJ];
-      if (mIJ == 0) {                                                     // Indeterminate form for mIJ = 0
-        for (int k=0; k<=order-mIJ; k++,index++) {
-          grads[index][0] = 0.;
-          grads[index][1] = 0.;
-          grads[index][2] = 2.*wg[k];
-        }
-      }
-      else if (mIJ == 1) {                                                // Indeterminate form for mIJ = 1
-        if (i == 0) {
+      if (validIJ(i,j)) {
+        int mIJ = std::max(i,j);
+        std::vector<double> &wf = wFcts[mIJ], &wg = wGrads[mIJ];
+        if (mIJ == 0) {                                                     // Indeterminate form for mIJ = 0
           for (int k=0; k<=order-mIJ; k++,index++) {
             grads[index][0] = 0.;
-            grads[index][1] = wf[k];
-            grads[index][2] = 2.*v*wg[k];
+            grads[index][1] = 0.;
+            grads[index][2] = 2.*wg[k];
           }
         }
-        else if (j == 0) {
-          for (int k=0; k<=order-mIJ; k++,index++) {
-            grads[index][0] = wf[k];
-            grads[index][1] = 0.;
-            grads[index][2] = 2.*u*wg[k];
+        else if (mIJ == 1) {                                                // Indeterminate form for mIJ = 1
+          if (i == 0) {
+            for (int k=0; k<=order-mIJ; k++,index++) {
+              grads[index][0] = 0.;
+              grads[index][1] = wf[k];
+              grads[index][2] = 2.*v*wg[k];
+            }
+          }
+          else if (j == 0) {
+            for (int k=0; k<=order-mIJ; k++,index++) {
+              grads[index][0] = wf[k];
+              grads[index][1] = 0.;
+              grads[index][2] = 2.*u*wg[k];
+            }
+          }
+          else {
+            for (int k=0; k<=order-mIJ; k++,index++) {
+              grads[index][0] = vhat*wf[k];
+              grads[index][1] = uhat*wf[k];
+              grads[index][2] = uhat*vhat*wf[k]+2.*uhat*v*wg[k];
+            }
           }
         }
-        else {
+        else {                                                              // General formula
+          double oMW = 1.-w;
+          double powM2 = pow(oMW, mIJ-2);
+          double powM1 = powM2*oMW;
           for (int k=0; k<=order-mIJ; k++,index++) {
-            grads[index][0] = vhat*wf[k];
-            grads[index][1] = uhat*wf[k];
-            grads[index][2] = uhat*vhat*wf[k]+2.*uhat*v*wg[k];
+            grads[index][0] = uGrads[i] * vFcts[j]  * wf[k] * powM1;
+            grads[index][1] = uFcts[i]  * vGrads[j] * wf[k] * powM1;
+            grads[index][2] = wf[k]*powM2*(u*uGrads[i]*vFcts[j]+v*uFcts[i]*vGrads[j])
+              + uFcts[i]*vFcts[j]*powM1*(2.*oMW*wg[k]-mIJ*wf[k]);
           }
         }
       }
-      else {                                                              // General formula
-        double oMW = 1.-w;
-        double powM2 = pow(oMW, mIJ-2);
-        double powM1 = powM2*oMW;
-        for (int k=0; k<=order-mIJ; k++,index++) {
-          grads[index][0] = uGrads[i] * vFcts[j]  * wf[k] * powM1;
-          grads[index][1] = uFcts[i]  * vGrads[j] * wf[k] * powM1;
-          grads[index][2] = wf[k]*powM2*(u*uGrads[i]*vFcts[j]+v*uFcts[i]*vGrads[j])
-                                    + uFcts[i]*vFcts[j]*powM1*(2.*oMW*wg[k]-mIJ*wf[k]);
-        }
-      }
     }
-  }
-
+  }  
 }
+  

Numeric/BergotBasis.h

@@ -20,7 +20,7 @@
 class BergotBasis {
  public:
 
-  BergotBasis(int p);
+  BergotBasis(int p,bool incpl=false);
   virtual ~BergotBasis();
 
   int size() const { const int n = order+1; return n*(n+1)*(2*n+1)/6; }
@@ -31,9 +31,12 @@ class BergotBasis {
 
   void initialize() {};
 
+  bool validIJ(int i,int j) const;
+
  private:
 
   int order; //!< maximal order of surrounding functional spaces (on triangle / quad)
+  bool incomplete; //!< serendipity interpolation
 
 };
 

Numeric/GaussQuadraturePyr.cpp

@@ -46,8 +46,6 @@ IntPt *getGQPyrPts(int order)
       int iU = (i - iW*nbPtUV2)/nbPtUV;
       int iV = (i - iW*nbPtUV2 - iU*nbPtUV);
 
-      // std::cout << "Points " << iU << " " << iV << " " << iW << std::endl;
-
       double wt = linWt[iU]*linWt[iV]*GJ20Wt[iW];
 
       double up = linPt[iU];

Numeric/pointsGenerators.cpp

@@ -137,6 +137,7 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
     points(i, 0) = duv + points(i, 0) * 2. / order;
     points(i, 1) = duv + points(i, 1) * 2. / order;
   }
+
   return points;
 }
 

Numeric/pyramidalBasis.cpp

@@ -10,12 +10,12 @@
 
 pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag), bergot(0)
 {
-  if (serendip) {
-    Msg::Warning("Serendipity pyramid not implemented yet");
+  if (serendip && order > 2) {
+    Msg::Warning("Serendipity pyramid for order %i not yet implemented",order);
     return;
   }
-
-  bergot = new BergotBasis(order);
+  
+  bergot = new BergotBasis(order,serendip);
 
   int num_points = points.size1();
 
@@ -36,17 +36,19 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag), bergot(0)
   int idx = 0;
   for (int i=0;i<=order;i++) {
     for (int j=0;j<=order;j++) {
-      for (int k=0;k<=order-std::max(i,j);k++,idx++) {
-        monomials(idx,0) = i;
-        monomials(idx,1) = j;
-        monomials(idx,2) = k;
-
-        for (int l=0;l<num_points;l++) {
-          double oneMinW = std::max(1e-14,1-points(l,2));
-          VDM(idx,l)  = std::pow(points(l,0),i);
-          VDM(idx,l) *= std::pow(points(l,1),j);
-          VDM(idx,l) *= std::pow(points(l,2),k);
-          VDM(idx,l) *= std::pow(oneMinW,std::max(i,j)-i-j);
+      if (bergot->validIJ(i,j)) {
+        for (int k=0;k<=order-std::max(i,j);k++,idx++) {
+          monomials(idx,0) = i;
+          monomials(idx,1) = j;
+          monomials(idx,2) = k;
+          
+          for (int l=0;l<num_points;l++) {
+            double oneMinW = std::max(1e-14,1-points(l,2));
+            VDM(idx,l)  = std::pow(points(l,0),i);
+            VDM(idx,l) *= std::pow(points(l,1),j);
+            VDM(idx,l) *= std::pow(points(l,2),k);
+              VDM(idx,l) *= std::pow(oneMinW,std::max(i,j)-i-j);
+          }
         }
       }
     }