From 0675a29bc1c2fccfc41455cf768e7dba42a59361 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Wed, 21 Nov 2012 10:05:52 +0000
Subject: [PATCH] Fixed non-standard C++ code in nodalBasis and pyramidalBasis
---
Numeric/BergotBasis.cpp | 28 +++++++++++++++-------------
Numeric/nodalBasis.cpp | 2 +-
Numeric/pyramidalBasis.cpp | 12 ++++++------
3 files changed, 22 insertions(+), 20 deletions(-)
diff --git a/Numeric/BergotBasis.cpp b/Numeric/BergotBasis.cpp
index b465082c78..ad9a2c0433 100644
--- a/Numeric/BergotBasis.cpp
+++ b/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";
}
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index e11af47cad..ada63ee201 100644
--- a/Numeric/nodalBasis.cpp
+++ b/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);
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index ae04985de2..4a67e1f0f7 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/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;
+
}
--
GitLab