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

Fixed bugs in shape functions for pyramidal elements (still buggy for w = 1).

parent 2ae1c81f
No related branches found
No related tags found
No related merge requests found
......@@ -14,52 +14,23 @@
BergotBasis::BergotBasis(int p): order(p)
{
// allocate function information and fill
iOrder = new int[size()];
jOrder = new int[size()];
kOrder = new int[size()];
LegendrePolynomials lp(p);
legendre = lp;
int index = 0;
for (int i=0;i<=order;i++) {
for (int j=0;j<=order;j++) {
int mIJ = std::max(i,j);
for (int k=0; k<=order-mIJ; k++,index++) {
iOrder[index] = i;
jOrder[index] = j;
kOrder[index] = k;
if (jacobi.find(mIJ) == jacobi.end()) {
JacobiPolynomials jp(2*mIJ+2,0,order);
jacobi[mIJ] = jp;
}
}
}
}
}
BergotBasis::~BergotBasis()
{
delete [] iOrder;
delete [] jOrder;
delete [] kOrder;
}
void BergotBasis::f(double u, double v, double w, double* val) const
{
double uhat = (w == 1.) ? 1. : u/(1.-w);
LegendrePolynomials legendre(order);
double uhat = (w == 1.) ? 1. : u/(1.-w);
std::vector<double> uFcts(order+1);
//double uFcts[order+1];
legendre.f(uhat,&(uFcts[0]));
double vhat = (w == 1.) ? 1. : v/(1.-w);
......@@ -67,14 +38,13 @@ void BergotBasis::f(double u, double v, double w, double* val) const
legendre.f(vhat,&(vFcts[0]));
double what = 2.*w - 1.;
std::map<int,double* > wFcts;
std::map<int,JacobiPolynomials>::const_iterator jIter=jacobi.begin();
for (;jIter!=jacobi.end();jIter++) {
int a = jIter->first;
wFcts[a] = new double[order + 1];
double* wf = wFcts[a];
jIter->second.f(what,wf);
std::vector<std::vector<double> > wFcts(order+1), wGrads(order+1);
for (int mIJ=0; mIJ<=order; mIJ++) {
int kMax = order-mIJ;
std::vector<double> &wf = wFcts[mIJ];
wf.resize(kMax+1);
JacobiPolynomials jacobi(2.*mIJ+2.,0.,kMax);
jacobi.f(what,&(wf[0]));
}
// recombine to find shape function values
......@@ -83,50 +53,60 @@ void BergotBasis::f(double u, double v, double w, double* val) const
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);
double* wf = (wFcts.find(mIJ))->second;
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;
}
}
}
jIter=jacobi.begin();
for (;jIter!=jacobi.end();jIter++) {
int a = jIter->first;
delete [] wFcts[a];
}
}
void BergotBasis::df(double u, double v, double w, double grads[][3]) const
{
std::vector<double> uFcts(order+1);
legendre.f(u,&(uFcts[0]));
std::vector<double> uGrads(order+1);
legendre.df(u,&(uGrads[0]));
LegendrePolynomials legendre(order);
std::vector<double> vFcts(order+1);
legendre.f(v,&(vFcts[0]));
// std::cout << "DBGTT: u = " << u << ", v = " << v << ", w = " << w << "\n";
std::vector<double> vGrads(order+1);
legendre.df(v,&(vGrads[0]));
double uhat = (w == 1.) ? 1. : u/(1.-w);
std::vector<double> uFcts(order+1), uGrads(order+1);
legendre.f(uhat,&(uFcts[0]));
legendre.df(uhat,&(uGrads[0]));
std::map<int,double* > wFcts;
std::map<int,double* > wGrads;
std::map<int,JacobiPolynomials>::const_iterator jIter=jacobi.begin();
double vhat = (w == 1.) ? 1. : v/(1.-w);
std::vector<double> vFcts(order+1), vGrads(order+1);
legendre.f(vhat,&(vFcts[0]));
legendre.df(vhat,&(vGrads[0]));
for (;jIter!=jacobi.end();jIter++) {
int a = jIter->first;
wFcts[a] = new double[order+1];
double* wf = wFcts[a];
jIter->second.f(w,wf);
wGrads[a] = new double[order+1];
double* wg = wGrads[a];
jIter->second.df(w,wg);
double what = 2.*w - 1.;
std::vector<std::vector<double> > wFcts(order+1), wGrads(order+1);
for (int mIJ=0; mIJ<=order; mIJ++) {
int kMax = order-mIJ;
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 {
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
......@@ -134,38 +114,46 @@ 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);
double* wf = (wFcts .find(mIJ))->second;
double* wg = (wGrads.find(mIJ))->second;
double wPowM2 = pow(1.-w, mIJ-2);
double wPowM1 = wPowM2*(1.-w);
double wPowM0 = wPowM1*(1.-w);
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++) {
grads[index][0] = uGrads[i] * vFcts[j] * wf[k] * wPowM1;
grads[index][1] = uFcts[i] * vGrads[j] * wf[k] * wPowM1;
grads[index][2] = 0;
grads[index][2] += u * uGrads[i] * vFcts[j] * wf[k] * wPowM2;
grads[index][2] += v * uFcts[i] * vGrads[j] * wf[k] * wPowM2;
grads[index][2] -= mIJ * uFcts[i] * vFcts[j] * wf[k] * wPowM1;
grads[index][2] += 2 * uFcts[i] * vFcts[j] * wg[k] * wPowM0;
if (mIJ == 0) {
grads[index][0] = 0.;
grads[index][1] = 0.;
grads[index][2] = 2.*wg[k];
}
else if (mIJ == 1) {
if (i == 0) {
grads[index][0] = 0.;
grads[index][1] = wf[k];
grads[index][2] = 2.*v*wg[k];
}
else if (j == 0) {
grads[index][0] = wf[k];
grads[index][1] = 0.;
grads[index][2] = 2.*u*wg[k];
}
else {
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 {
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";
}
}
jIter=jacobi.begin();
for (;jIter!=jacobi.end();jIter++) {
int a = jIter->first;
delete [] wFcts[a];
delete [] wGrads[a];
}
}
......@@ -26,17 +26,8 @@ class BergotBasis {
void initialize() {};
private:
int order; //!< maximal order of surrounding functional spaces (on triangle / quad)
int *iOrder; //!< order of \f$\hat \xi \f$ polynomial
int *jOrder; //!< order of \f$\hat \eta \f$ polynomial
int *kOrder; //!< order of \f$\hat \zeta \f$ polynomial
//! list of Legendre polynomials up to order p
LegendrePolynomials legendre;
//! list of Jacobi polynomials up to order p in function of index i (\f$ \alpha = 2*i + 2\f$)
std::map<int,JacobiPolynomials> jacobi;
int order; //!< maximal order of surrounding functional spaces (on triangle / quad)
};
......
......@@ -19,10 +19,10 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
num_points -= (order-2)*((order-2)+1)*(2*(order-2)+1)/6;
}
VDMinv = new fullMatrix<double>(num_points, num_points);
VDMinv.resize(num_points, num_points);
// Invert the Vandermonde matrix
fullMatrix<double> VDM = fullMatrix<double>(num_points, num_points);
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);
......@@ -30,7 +30,7 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
VDM(i,j) = f[i];
}
}
VDM.invert(*VDMinv);
VDM.invert(VDMinv);
}
......@@ -48,14 +48,16 @@ void pyramidalBasis::f(double u, double v, double w, double *val) const
const int N = bergot->size();
double f[N];
bergot->f(u, v, w, f);
double *fval = new double[N];
bergot->f(u, v, w, fval);
for (int i = 0; i < N; i++) {
val[i] = 0.;
for (int j = 0; j < N; j++) val[i] += (*VDMinv)(i,j)*f[j];
for (int j = 0; j < N; j++) val[i] += VDMinv(i,j)*fval[j];
}
delete[] fval;
}
......@@ -66,17 +68,18 @@ void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
const int N = bergot->size(), NPts = coord.size1();
sf.resize(NPts, N);
double f[N];
fullVector<double> fVect(f,N);
double *fval = new double[N];
for (int iPt=0; iPt<NPts; iPt++) {
bergot->f(coord(iPt,0), coord(iPt,0), coord(iPt,0), f);
bergot->f(coord(iPt,0), coord(iPt,1), coord(iPt,2), fval);
for (int i = 0; i < N; i++) {
sf(iPt,i) = 0.;
for (int j = 0; j < N; j++) sf(iPt,i) += (*VDMinv)(i,j)*f[j];
for (int j = 0; j < N; j++) sf(iPt,i) += VDMinv(i,j)*fval[j];
}
}
delete[] fval;
}
......@@ -86,18 +89,20 @@ void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
const int N = bergot->size();
double df[N][3];
bergot->df(u, v, w, df);
double (*dfval)[3] = new double[N][3];
bergot->df(u, v, w, dfval);
for (int i = 0; i < N; i++) {
grads[i][0] = 0.; grads[i][1] = 0.; grads[i][2] = 0.;
for (int j = 0; j < N; j++) {
grads[i][0] += (*VDMinv)(i,j)*df[j][0];
grads[i][1] += (*VDMinv)(i,j)*df[j][1];
grads[i][2] += (*VDMinv)(i,j)*df[j][2];
grads[i][0] += VDMinv(i,j)*dfval[j][0];
grads[i][1] += VDMinv(i,j)*dfval[j][1];
grads[i][2] += VDMinv(i,j)*dfval[j][2];
}
}
delete[] dfval;
}
......@@ -107,7 +112,7 @@ void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) cons
const int N = bergot->size(), NPts = coord.size1();
double dfv[N][3];
double (*dfv)[3] = new double[N][3];
dfm.resize (N, 3*NPts, false);
for (int iPt=0; iPt<NPts; iPt++) {
......@@ -119,4 +124,6 @@ void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) cons
}
}
delete[] dfv;
}
......@@ -17,7 +17,7 @@ class pyramidalBasis: public nodalBasis
{
private:
// Inverse of the Vandermonde matrix
fullMatrix<double>* VDMinv;
fullMatrix<double> VDMinv;
// Orthogonal basis for the pyramid
BergotBasis *bergot;
......
......@@ -12,6 +12,7 @@
#include "fullMatrix.h"
#include "nodalBasis.h"
#include "polynomialBasis.h"
#include "pyramidalBasis.h"
%}
%include "GaussIntegration.h"
......@@ -21,3 +22,4 @@
%template(fullVectorDouble) fullVector<double>;
%include "nodalBasis.h"
%include "polynomialBasis.h"
%include "pyramidalBasis.h"
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment