Select Git revision
BergotBasis.cpp
BergotBasis.cpp 5.74 KiB
// Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
#include <cmath>
#include "BergotBasis.h"
#include "MElement.h"
#include "orthogonalBasis.h"
BergotBasis::BergotBasis(int p, bool incpl) : order(p), incomplete(incpl)
{
if(incomplete && order > 2) {
Msg::Error("Incomplete pyramids of order %i not yet implemented", order);
}
}
BergotBasis::~BergotBasis() {}
bool BergotBasis::validIJ(int i, int j) const
{
if(!incomplete) return (i <= order) && (j <= order);
if(i + j <= order) return true;
if(i + j == order + 1) return i == 1 || j == 1;
return false;
}
// 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
{
// Compute Legendre polynomials at uhat
double uhat = (w == 1.) ? 0. : u / (1. - w);
std::vector<double> uFcts(order + 1);
LegendrePolynomials::f(order, uhat, &(uFcts[0]));
// Compute Legendre polynomials at vhat
double vhat = (w == 1.) ? 0. : v / (1. - w);
std::vector<double> vFcts(order + 1);
LegendrePolynomials::f(order, 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++) {
int kMax = order - mIJ;
std::vector<double> &wf = wFcts[mIJ];
wf.resize(kMax + 1);
JacobiPolynomials::f(kMax, 2. * mIJ + 2., 0., what, &(wf[0]));
}
// Recombine to find shape function values
int index = 0;
for(int i = 0; i <= order; i++) {
for(int j = 0; j <= order; j++) {
if(validIJ(i, 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;
}
}
}
}
// Derivatives of Bergot basis functions for coordinates (u,v,w) in the unit