Fixed non-standard C++ code in nodalBasis and pyramidalBasis

Numeric/BergotBasis.cpp

@@ -89,22 +89,24 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const
     std::vector<double> &wf = wFcts[mIJ], &wg = wGrads[mIJ];
     wf.resize(kMax+1);
     wg.resize(kMax+1);
-//    if (what == 1.) {
-//      const int alpha = 2*mIJ+2, fMax = kMax+alpha;
-//      std::vector<double> fact(fMax);
-//      fact[0] = 0.;
-//      for (int i=1;i<=fMax;i++) fact[i] = i*fact[i-1];
-//      wf[0] = 1.; wg[0] = 0.;
-//      for (int k=1;k<=kMax;k++) {
-//        wf[k] = fact[k+alpha]/(fact[alpha]*fact[k]);
-//        wg[k] = 0.5*(k+alpha+2)*fact[k+alpha]/(fact[alpha+1]*fact[k-1]);
-//      }
-//    }
-//    else {
+    if (what == 1.) {
+      const int alpha = 2*mIJ+2, fMax = std::max(kMax,1)+alpha;
+      std::vector<double> fact(fMax);
+      fact[0] = 1.;
+      for (int i=1;i<=fMax;i++) fact[i] = i*fact[i-1];
+      wf[0] = 1.; wg[0] = 0.;
+      for (int k=1;k<=kMax;k++) {
+//        std::cout << "DBGTT: calc. fMax = " << fMax << ", alpha = " << alpha << ", k = " << k << "\n";
+        wf[k] = fact[k+alpha]/(fact[alpha]*fact[k]);
+        wg[k] = 0.5*(k+alpha+2)*fact[k+alpha]/(fact[alpha+1]*fact[k-1]);
+//        std::cout << "DBGTT:     -> wf[k] = " << wf[k] << ", wg[k] = " << wg[k] << "\n";
+      }
+    }
+    else {
       JacobiPolynomials jacobi(2.*mIJ+2.,0.,kMax);
       jacobi.f(what,&(wf[0]));
       jacobi.df(what,&(wg[0]));
-//    }
+    }
 //    for (int k=0; k<=order-mIJ; k++) std::cout << " -> mIJ = " << mIJ << ", k = " << k
 //        << " -> wf[k] = " << wf[k] << ", wg[k] = " << wg[k] << "\n";
   }

Numeric/nodalBasis.cpp

@@ -765,7 +765,7 @@ static fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
       }
 
       // Volume
-      if (!serendip and order > 2) {
+      if ((!serendip) && (order > 2)) {
         fullMatrix<double> volume_points = gmshGeneratePointsPyramid(order - 3, false);
         // scale to order-3/order
         fullMatrix<double> T(3,3);

Numeric/pyramidalBasis.cpp

@@ -15,23 +15,23 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
 
   int n = order+1;
   int num_points = n*(n+1)*(2*n+1)/6;
-  if (serendip and order > 2) {
+  if (serendip && (order > 2)) {
     num_points -= (order-2)*((order-2)+1)*(2*(order-2)+1)/6;
   }
 
   VDMinv.resize(num_points, num_points);
+  double *fval = new double[num_points];
 
   // Invert the Vandermonde matrix
   fullMatrix<double> VDM(num_points, num_points);
   for (int j = 0; j < num_points; j++) {
-    double f[num_points];
-    bergot->f(points(j,0), points(j,1), points(j, 2), f);
-    for (int i = 0; i < num_points; i++) {
-      VDM(i,j) = f[i];
-    }
+    bergot->f(points(j,0), points(j,1), points(j, 2), fval);
+    for (int i = 0; i < num_points; i++) VDM(i,j) = fval[i];
   }
   VDM.invert(VDMinv);
 
+  delete[] fval;
+
 }