From 704d10d77cce4e0c2b291e532b59e5db8c247e7c Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Wed, 21 Nov 2012 14:32:44 +0000
Subject: [PATCH] Fixed and cleaned pyramidal basis functions. Added test
script.
---
Numeric/BergotBasis.cpp | 109 +++++++++------
benchmarks/misc/testPyramidalBasis.py | 194 ++++++++++++++++++++++++++
2 files changed, 261 insertions(+), 42 deletions(-)
create mode 100644 benchmarks/misc/testPyramidalBasis.py
diff --git a/Numeric/BergotBasis.cpp b/Numeric/BergotBasis.cpp
index ad9a2c0433..8611c1e0a7 100644
--- a/Numeric/BergotBasis.cpp
+++ b/Numeric/BergotBasis.cpp
@@ -12,6 +12,30 @@
+namespace {
+
+// Private function to compute value and gradients of Jacobi polynomial at 1 (the general
+// implementation in class jacobiPolynomials does not handle the indefinite form)
+inline void jacobiDiffOne(int alpha, int beta, int n, double *wf, double *wg)
+{
+
+ const int fMax = std::max(n,1)+alpha;
+ std::vector<double> fact(fMax+1);
+ 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<=n;k++) {
+ wf[k] = fact[k+alpha]/(fact[alpha]*fact[k]);
+ wg[k] = 0.5*(k+alpha+beta+1)*fact[k+alpha]/(fact[alpha+1]*fact[k-1]);
+ }
+
+}
+
+}
+
+
+
+
BergotBasis::BergotBasis(int p): order(p)
{
}
@@ -24,19 +48,26 @@ BergotBasis::~BergotBasis()
+// 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)
+// with i,j = 0...order and k = 0...order-max(i,j)
+// and (uhat,vhat,what) = (u/(1-w),v/(1-w),2*w-1) reduced coordinates in the unit cube [-1,1]^3
void BergotBasis::f(double u, double v, double w, double* val) const
{
LegendrePolynomials legendre(order);
- double uhat = (w == 1.) ? 1. : u/(1.-w);
+ // Compute Legendre polynomials at uhat
+ double uhat = (w == 1.) ? 0. : u/(1.-w);
std::vector<double> uFcts(order+1);
legendre.f(uhat,&(uFcts[0]));
- double vhat = (w == 1.) ? 1. : v/(1.-w);
+ // Compute Legendre polynomials at vhat
+ double vhat = (w == 1.) ? 0. : v/(1.-w);
std::vector<double> vFcts(order+1);
legendre.f(vhat,&(vFcts[0]));
+ // Compute Jacobi polynomials at what
double what = 2.*w - 1.;
std::vector<std::vector<double> > wFcts(order+1), wGrads(order+1);
for (int mIJ=0; mIJ<=order; mIJ++) {
@@ -47,17 +78,14 @@ void BergotBasis::f(double u, double v, double w, double* val) const
jacobi.f(what,&(wf[0]));
}
- // recombine to find shape function values
-
+ // Recombine to find shape function values
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;
- }
+ for (int k=0; k<=order-mIJ; k++,index++) val[index] = uFcts[i] * vFcts[j] * wf[k] * fact;
}
}
@@ -65,23 +93,31 @@ void BergotBasis::f(double u, double v, double w, double* val) const
+// Derivatives of Bergot basis functions for coordinates (u,v,w) in the unit pyramid:
+// dfdu = L'_i(uhat)*L_j(vhat)*(1-w)^(max(i,j)-1)*P_k^{2*max(i,j)+2,0}(what)
+// dfdv = L_i(uhat)*L'_j(vhat)*(1-w)^(max(i,j)-1)*P_k^{2*max(i,j)+2,0}(what)
+// dfdw = (1-w)^(max(i,j)-2)*P_k^{2*max(i,j)+2,0}(what)*(u*L'_i(uhat)*L_j(vhat)+v*L_i(uhat)*L'_j(vhat))
+// + u*v*(1-w)^(max(i,j)-1)*(2*(1-w)*P'_k^{2*max(i,j)+2,0}(what)-max(i,j)*P_k^{2*max(i,j)+2,0}(what))
+// with i,j = 0...order and k = 0...order-max(i,j)
+// and (uhat,vhat,what) = (u/(1-w),v/(1-w),2*w-1) reduced coordinates in the unit cube [-1,1]^3
void BergotBasis::df(double u, double v, double w, double grads[][3]) const
{
LegendrePolynomials legendre(order);
-// std::cout << "DBGTT: u = " << u << ", v = " << v << ", w = " << w << "\n";
-
- double uhat = (w == 1.) ? 1. : u/(1.-w);
+ // Compute Legendre polynomials and derivatives at uhat
+ double uhat = (w == 1.) ? 0. : u/(1.-w);
std::vector<double> uFcts(order+1), uGrads(order+1);
legendre.f(uhat,&(uFcts[0]));
legendre.df(uhat,&(uGrads[0]));
- double vhat = (w == 1.) ? 1. : v/(1.-w);
+ // Compute Legendre polynomials and derivatives at vhat
+ double vhat = (w == 1.) ? 0. : v/(1.-w);
std::vector<double> vFcts(order+1), vGrads(order+1);
legendre.f(vhat,&(vFcts[0]));
legendre.df(vhat,&(vGrads[0]));
+ // Compute Jacobi polynomials and derivatives at what
double what = 2.*w - 1.;
std::vector<std::vector<double> > wFcts(order+1), wGrads(order+1);
for (int mIJ=0; mIJ<=order; mIJ++) {
@@ -89,71 +125,60 @@ 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 = 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";
- }
- }
+ if (what == 1.) jacobiDiffOne(2*mIJ+2,0,kMax,&(wf[0]),&(wg[0]));
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";
}
- // now recombine to find the shape function gradients
-
+ // Recombine to find the shape function gradients
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];
- double oMW = 1.-w;
- double powM2 = ((w == 1.) && (mIJ < 2)) ? 0. : pow(oMW, mIJ-2);
- double powM1 = powM2*oMW;
- for (int k=0; k<=order-mIJ; k++,index++) {
- if (mIJ == 0) {
+ 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) {
- if (i == 0) {
+ }
+ 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) {
+ }
+ 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 {
+ }
+ 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] = 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]);
}
-//// std::cout << " -> i = " << i << ", j = " << j << ", k = " << k << " -> grad = ("
-//// << grads[index][0] << ", " << grads[index][1] << ", "<< grads[index][2] << ")\n";
-// std::cout << " -> i = " << i << ", j = " << j << ", k = " << k
-// << " -> wf[k] = " << wf[k] << ", wg[k] = " << wg[k] << "\n";
}
}
}
diff --git a/benchmarks/misc/testPyramidalBasis.py b/benchmarks/misc/testPyramidalBasis.py
new file mode 100644
index 0000000000..b2c3b1c45e
--- /dev/null
+++ b/benchmarks/misc/testPyramidalBasis.py
@@ -0,0 +1,194 @@
+from dgpy import *
+import random
+import math
+
+
+
+# Parameters
+Nr = 100
+p = 8
+
+
+
+# Evaluate polynomial
+def polyEval(FCT, p, coeff, uvw):
+ for n in range(uvw.size1()):
+ u = uvw(n,0)
+ v = uvw(n,1)
+ w = uvw(n,2)
+ FCT.set(n,0)
+ ind = 0
+ for i in range(p+1):
+ upi = math.pow(u,i)
+ for j in range(p+1-i):
+ vpj = math.pow(v,j)
+ for k in range(p+1-i-j):
+ FCT.set(n,FCT(n)+coeff(ind)*upi*vpj*math.pow(w,k))
+ ind += 1
+
+
+
+# Evaluate gradient of polynomial
+def gradPolyEval(FCT, p, coeff, uvw):
+ for n in range(uvw.size1()):
+ u = uvw(n,0)
+ v = uvw(n,1)
+ w = uvw(n,2)
+ FCT.set(3*n,0)
+ FCT.set(3*n+1,0)
+ FCT.set(3*n+2,0)
+ ind = 0
+ for i in range(p+1):
+ upi = math.pow(u,i)
+ if (i == 0):
+ dupi = 0
+ else:
+ dupi = i*math.pow(u,i-1)
+ for j in range(p+1-i):
+ vpj = math.pow(v,j)
+ if (j == 0):
+ dvpj = 0
+ else:
+ dvpj = j*math.pow(v,j-1)
+ for k in range(p+1-i-j):
+ wpk = math.pow(w,k)
+ if (k == 0):
+ dwpk = 0
+ else:
+ dwpk = k*math.pow(w,k-1)
+ FCT.set(3*n,FCT(3*n)+coeff(ind)*dupi*vpj*wpk)
+ FCT.set(3*n+1,FCT(3*n+1)+coeff(ind)*upi*dvpj*wpk)
+ FCT.set(3*n+2,FCT(3*n+2)+coeff(ind)*upi*vpj*dwpk)
+ ind += 1
+
+
+
+# Evaluate approx. on nodal basis
+def interp(FCT, basis, coeff, uvw):
+ sf = fullMatrixDouble(uvw.size1(),coeff.size())
+ basis.f(uvw,sf)
+ for i in range(uvw.size1()):
+ FCT.set(i,0)
+ for j in range(coeff.size()):
+ FCT.set(i,FCT(i)+coeff(j)*sf(i,j))
+
+
+
+# Evaluate gradient of approx. on nodal basis
+def gradInterp(FCT, basis, coeff, uvw):
+ gsf = fullMatrixDouble(coeff.size(),3*uvw.size1())
+ basis.df(uvw,gsf)
+ for i in range(3*uvw.size1()):
+ FCT.set(i,0)
+ for j in range(coeff.size()):
+ FCT.set(i,FCT(i)+coeff(j)*gsf(j,i))
+
+
+
+# Test if point is in unit pyramid
+def isInPyr(u,v,w):
+ if ((w >= 0) and (w < 1)):
+ uh = u/(1-w)
+ vh = v/(1-w)
+ return ((uh >= -1) and (uh <= 1) and (vh >= -1) and (vh <= 1))
+ elif (w == 1):
+ return ((u == 0) and (v == 0))
+ else:
+ return False
+
+
+
+# Generate polynomial of order p with random coefficients
+Np = (p+1)*(p+2)*(p+3)/6
+poly = fullVectorDouble(Np)
+for i in range(Np):
+ poly.set(i,random.random())
+
+# Get appropriate basis
+if (p == 1):
+ basis = pyramidalBasis(MSH_PYR_5)
+elif (p == 2):
+ basis = pyramidalBasis(MSH_PYR_14)
+elif (p == 3):
+ basis = pyramidalBasis(MSH_PYR_30)
+elif (p == 4):
+ basis = pyramidalBasis(MSH_PYR_55)
+elif (p == 5):
+ basis = pyramidalBasis(MSH_PYR_91)
+elif (p == 6):
+ basis = pyramidalBasis(MSH_PYR_140)
+elif (p == 7):
+ basis = pyramidalBasis(MSH_PYR_204)
+elif (p == 8):
+ basis = pyramidalBasis(MSH_PYR_285)
+elif (p == 9):
+ basis = pyramidalBasis(MSH_PYR_385)
+
+print("Basis:\n -> type = %i\n -> parentType = %i\n -> order = %i\n -> points = %i\n -> dimension = %i\n -> numFaces = %i\n -> serendip = %s" %
+ (basis.type, basis.parentType, basis.order, basis.points.size1(), basis.dimension, basis.numFaces, basis.serendip))
+
+#print("%i points:" % basis.points.size1())
+#for i in range(basis.points.size1()) :
+# print(" -> (%g, %g, %g)" % (basis.points(i,0), basis.points(i,1), basis.points(i,2)))
+
+#sf = fullMatrixDouble(basis.points.size1(),5)
+#basis.f(basis.points,sf)
+#print("SF:")
+#for i in range(5) :
+# print(" -> %g %g %g %g %g" % (sf(i,0), sf(i,1), sf(i,2), sf(i,3), sf(i,4)))
+
+# nodal coefficients of function to be evaluated
+nodalCoeff = fullVectorDouble(basis.points.size1())
+polyEval(nodalCoeff,p,poly,basis.points)
+#print("%i nodalCoeffs:" % nodalCoeff.size())
+#for i in range(nodalCoeff.size()) :
+# print(" -> %g" % nodalCoeff(i))
+
+# Generate Nr random points in pyramid (well, Nr-1 random plus (0,0,1) to test indeterminate form)
+uvwr = fullMatrixDouble(Nr,3)
+uvwr.set(0,0,0)
+uvwr.set(0,1,0)
+uvwr.set(0,2,1)
+for i in range(1,Nr):
+ uvwr.set(i,0,1000)
+ uvwr.set(i,1,1000)
+ uvwr.set(i,2,1000)
+ while (not isInPyr(uvwr(i,0),uvwr(i,1),uvwr(i,2))):
+ uvwr.set(i,0,random.random())
+ uvwr.set(i,1,random.random())
+ uvwr.set(i,2,random.random())
+
+# Exact value and gradient of function at random points
+valExact = fullVectorDouble(Nr)
+polyEval(valExact,p,poly,uvwr)
+gradExact = fullVectorDouble(3*Nr)
+gradPolyEval(gradExact,p,poly,uvwr)
+
+# Approx. value and gradient of function at random points
+valApp = fullVectorDouble(Nr)
+interp(valApp,basis,nodalCoeff,uvwr)
+gradApp = fullVectorDouble(3*Nr)
+gradInterp(gradApp,basis,nodalCoeff,uvwr)
+
+# Compute error
+avgErrVal = 0
+avgErrGradX = 0
+avgErrGradY = 0
+avgErrGradZ = 0
+for i in range(Nr):
+ avgErrVal += (valApp(i)-valExact(i))*(valApp(i)-valExact(i))
+ avgErrGradX += (gradApp(3*i)-gradExact(3*i))*(gradApp(3*i)-gradExact(3*i))
+ avgErrGradY += (gradApp(3*i+1)-gradExact(3*i+1))*(gradApp(3*i+1)-gradExact(3*i+1))
+ avgErrGradZ += (gradApp(3*i+2)-gradExact(3*i+2))*(gradApp(3*i+2)-gradExact(3*i+2))
+avgErrVal = math.sqrt(avgErrVal)/Nr
+avgErrGradX = math.sqrt(avgErrGradX)/Nr
+avgErrGradY = math.sqrt(avgErrGradY)/Nr
+avgErrGradZ = math.sqrt(avgErrGradZ)/Nr
+
+# Print error
+print("Check on %i random points:" % Nr)
+print(" -> average error value = %g" % avgErrVal)
+print(" -> average error X gradient = %g" % avgErrGradX)
+print(" -> average error Y gradient = %g" % avgErrGradY)
+print(" -> average error Z gradient = %g" % avgErrGradZ)
+
--
GitLab