Skip to content
Snippets Groups Projects
Commit 704d10d7 authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Fixed and cleaned pyramidal basis functions. Added test script.

parent 85ad61fc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/BergotBasis.cpp
+ 67
− 42
View file @ 704d10d7
......@@ -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";
}
}
}
......
benchmarks/misc/testPyramidalBasis.py 0 → 100644
+ 194
− 0
View file @ 704d10d7
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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment