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

Fixed indeterminate form in Legendre and Jacobi polynomials. Fixed pyramidal basis functions.

parent 35074cd4
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 113 additions and 100 deletions
Geo/MPyramid.h
+ 0
− 44
View file @ 5e4793b6
......@@ -124,50 +124,6 @@ class MPyramid : public MElement {
MVertex *tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
}
virtual int getVolumeSign();
/*virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
double r = (w != 1.) ? (u * v * w / (1. - w)) : 0.;
s[0] = 0.25 * ((1. - u) * (1. - v) - w + r);
s[1] = 0.25 * ((1. + u) * (1. - v) - w - r);
s[2] = 0.25 * ((1. + u) * (1. + v) - w + r);
s[3] = 0.25 * ((1. - u) * (1. + v) - w - r);
s[4] = w;
}
*/
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
{
if(w == 1.) {
s[0][0] = -0.25 ;
s[0][1] = -0.25 ;
s[0][2] = -0.25 ;
s[1][0] = 0.25 ;
s[1][1] = -0.25 ;
s[1][2] = -0.25 ;
s[2][0] = 0.25 ;
s[2][1] = 0.25 ;
s[2][2] = -0.25 ;
s[3][0] = -0.25 ;
s[3][1] = 0.25 ;
s[3][2] = -0.25 ;
}
else{
s[0][0] = 0.25 * ( -(1. - v) + v * w / (1. - w)) ;
s[0][1] = 0.25 * ( -(1. - u) + u * w / (1. - w)) ;
s[0][2] = 0.25 * ( -1. + u * v / (1. - w) + u * v * w / (1. - w) / (1. - w)) ;
s[1][0] = 0.25 * ( (1. - v) - v * w / (1. - w)) ;
s[1][1] = 0.25 * ( -(1. + u) - u * w / (1. - w)) ;
s[1][2] = 0.25 * ( -1. - u * v / (1. - w) - u * v * w / (1. - w) / (1. - w)) ;
s[2][0] = 0.25 * ( (1. + v) + v * w / (1. - w)) ;
s[2][1] = 0.25 * ( (1. + u) + u * w / (1. - w)) ;
s[2][2] = 0.25 * ( -1. + u * v / (1. - w) + u * v * w / (1. - w) / (1. - w)) ;
s[3][0] = 0.25 * ( -(1. + v) - v * w / (1. - w)) ;
s[3][1] = 0.25 * ( (1. - u) - u * w / (1. - w)) ;
s[3][2] = 0.25 * ( -1. - u * v / (1. - w) - u * v * w / (1. - w) / (1. - w)) ;
}
s[4][0] = 0.;
s[4][1] = 0.;
s[4][2] = 1.;
}
virtual void getNode(int num, double &u, double &v, double &w)
{
switch(num) {
......
Numeric/BergotBasis.cpp
+ 6
− 36
View file @ 5e4793b6
......@@ -3,36 +3,9 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include <cmath>
#include "BergotBasis.h"
/*BergotBasis::BergotBasis() {
}*/
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]);
}
}
}
#include "MElement.h"
......@@ -125,13 +98,10 @@ 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.) 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]));
}
}
// Recombine to find the shape function gradients
int index = 0;
......
Numeric/JacobianBasis.cpp
+ 22
− 0
View file @ 5e4793b6
......@@ -1025,6 +1025,26 @@ static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &poi
static void generateGradShapesPyramid(JacobianBasis &jfs, const fullMatrix<double> &points, const pyramidalBasis *F)
{
fullMatrix<double> allDPsi;
F->df(points, allDPsi);
const int NBez = points.size1(), NLag = allDPsi.size1();
jfs.gradShapeMatX.resize(NBez,NLag);
jfs.gradShapeMatY.resize(NBez,NLag);
jfs.gradShapeMatZ.resize(NBez,NLag);
for (int i=0; i<NBez; i++)
for (int j=0; j<NLag; j++) {
jfs.gradShapeMatX(i,j) = allDPsi(j,3*i);
jfs.gradShapeMatY(i,j) = allDPsi(j,3*i+1);
jfs.gradShapeMatZ(i,j) = allDPsi(j,3*i+2);
}
}
std::map<int, bezierBasis> bezierBasis::_bbs;
const bezierBasis *bezierBasis::find(int tag)
{
......@@ -1142,6 +1162,8 @@ const JacobianBasis *JacobianBasis::find(int tag)
J.bezier = bezierBasis::find(MSH_PYR_14);
break;
}
pyramidalBasis *F = (pyramidalBasis*)BasisFactory::create(tag);
generateGradShapesPyramid(J, J.bezier->points, F);
}
else {
const int parentType = MElement::ParentTypeFromTag(tag), order = MElement::OrderFromTag(tag);
......
Numeric/jacobiPolynomials.cpp
+ 41
− 1
View file @ 5e4793b6
......@@ -4,6 +4,10 @@
// bugs and problems to <gmsh@geuz.org>.
#include "jacobiPolynomials.h"
#include <cmath>
#include <iostream>
inline double Pochhammer(double x,int n)
{
......@@ -12,11 +16,17 @@ inline double Pochhammer(double x,int n)
return result;
}
JacobiPolynomials::JacobiPolynomials(double a, double b, int o):
alpha(a),beta(b),n(o),alphaPlusBeta(a+b),a2MinusB2(a*a-b*b) {}
JacobiPolynomials::~JacobiPolynomials() {;}
void JacobiPolynomials::f(double u, double *val) const
{
val[0] = 1.;
......@@ -37,21 +47,51 @@ void JacobiPolynomials::f(double u, double *val) const
}
}
void JacobiPolynomials::df(double u, double *val) const
{
// Indeterminate form for u == -1 and u == 1
// TODO: Extend to non-integer alpha & beta?
if ((u == 1.) || (u == -1.)) {
// alpha or beta in formula, depending on u
int coeff = (u == 1.) ? (int)alpha : (int)beta;
// Compute factorial
const int fMax = std::max(n,1)+coeff;
std::vector<double> fact(fMax+1);
fact[0] = 1.;
for (int i=1;i<=fMax;i++) fact[i] = i*fact[i-1];
// Compute formula (with appropriate sign at even orders for u == -1)
val[0] = 0.;
for (int k=1;k<=n;k++) val[k] = 0.5*(k+alphaPlusBeta+1)*fact[k+coeff]/(fact[coeff+1]*fact[k-1]);
if ((u == -1.) && (n >= 2)) for (int k=2;k<=n;k+=2) val[k] = -val[k];
return;
}
// Now general case
// Values of the Jacobi polynomials
std::vector<double> tmp(n+1);
f(u,&(tmp[0]));
// First 2 values of the derivatives
val[0] = 0;
if (n>=1) val[1] = 0.5*(alphaPlusBeta + 2.);
// Values of the derivative for orders > 1 computed from the values of the polynomials
for (int i=2;i<=n;i++) {
double ii = (double) i;
double aa = (2.*ii + alphaPlusBeta);
double g2 = aa*(1.-u*u);
double g1 = ii*(alpha - beta - aa*u);
double g0 = 2.*(ii+alpha)*(ii+beta);
val[i] = (g1 * tmp[i] + g0 * tmp[i-1])/g2;
}
}
Numeric/legendrePolynomials.cpp
+ 26
− 3
View file @ 5e4793b6
......@@ -5,11 +5,16 @@
#include "legendrePolynomials.h"
LegendrePolynomials::LegendrePolynomials(int o): n(o) {
}
LegendrePolynomials::LegendrePolynomials(int o): n(o) {}
LegendrePolynomials::~LegendrePolynomials() {}
void LegendrePolynomials::f(double u, double *val) const
{
val[0] = 1;
......@@ -26,18 +31,36 @@ void LegendrePolynomials::f(double u, double *val) const
}
}
void LegendrePolynomials::df(double u, double *val) const
{
// Indeterminate form for u == -1 and u == 1
if ((u == 1.) || (u == -1.)) {
for (int k=0;k<=n;k++) val[k] = 0.5*k*(k+1);
if ((u == -1.) && (n >= 2)) for (int k=2;k<=n;k+=2) val[k] = -val[k];
return;
}
// Now general case
// Values of the Legendre polynomials
std::vector<double> tmp(n+1);
f(u,&(tmp[0]));
// First value of the derivative
val[0] = 0;
double g2 = (1.-u*u);
// Values of the derivative for orders > 1 computed from the values of the polynomials
for (int i=1;i<=n;i++) {
double g1 = - u*i;
double g0 = (double) i;
val[i] = (g1 * tmp[i] + g0 * tmp[i-1])/g2;
}
}
Numeric/nodalBasis.h
+ 2
− 2
View file @ 5e4793b6
......@@ -21,11 +21,11 @@ class nodalBasis {
// Basis functions evaluation
virtual void f(double u, double v, double w, double *sf) const {Msg::Fatal("Not implemented");};
virtual void f(fullMatrix<double> &coord, fullMatrix<double> &sf) const {Msg::Fatal("Not implemented");};
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const {Msg::Fatal("Not implemented");};
// Basis functions gradients evaluation
virtual void df(double u, double v, double w, double grads[][3]) const {Msg::Fatal("Not implemented");};
virtual void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const {Msg::Fatal("Not implemented");};
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const {Msg::Fatal("Not implemented");};
virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");};
virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");};
......
Numeric/polynomialBasis.cpp
+ 2
− 2
View file @ 5e4793b6
......@@ -404,7 +404,7 @@ void polynomialBasis::f(double u, double v, double w, double *sf) const
void polynomialBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf) const
void polynomialBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const
{
double p[1256];
sf.resize (coord.size1(), coefficients.size1());
......@@ -418,7 +418,7 @@ void polynomialBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf) const
void polynomialBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
void polynomialBasis::df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const
{
double dfv[1256][3];
dfm.resize (coefficients.size1(), coord.size1() * 3, false);
......
Numeric/polynomialBasis.h
+ 2
− 2
View file @ 5e4793b6
......@@ -89,8 +89,8 @@ class polynomialBasis : public nodalBasis
~polynomialBasis();
virtual void f(double u, double v, double w, double *sf) const;
virtual void f(fullMatrix<double> &coord, fullMatrix<double> &sf) const;
virtual void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void ddf(double u, double v, double w, double hess[][3][3]) const;
virtual void dddf(double u, double v, double w, double third[][3][3][3]) const;
......
Numeric/pyramidalBasis.cpp
+ 2
− 2
View file @ 5e4793b6
......@@ -62,7 +62,7 @@ void pyramidalBasis::f(double u, double v, double w, double *val) const
void pyramidalBasis::f(fullMatrix<double> &coord, fullMatrix<double> &sf)
void pyramidalBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf)
{
const int N = bergot->size(), NPts = coord.size1();
......@@ -107,7 +107,7 @@ void pyramidalBasis::df(double u, double v, double w, double grads[][3]) const
void pyramidalBasis::df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const
void pyramidalBasis::df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const
{
const int N = bergot->size(), NPts = coord.size1();
......
Numeric/pyramidalBasis.h
+ 2
− 2
View file @ 5e4793b6
......@@ -27,9 +27,9 @@ class pyramidalBasis: public nodalBasis
~pyramidalBasis();
virtual void f(double u, double v, double w, double *val) const;
virtual void f(fullMatrix<double> &coord, fullMatrix<double> &sf);
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf);
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
};
......
benchmarks/misc/testPyramidalBasis.py
+ 8
− 6
View file @ 5e4793b6
......@@ -6,7 +6,7 @@ import math
# Parameters
Nr = 100
p = 8
p = 2
......@@ -144,12 +144,14 @@ polyEval(nodalCoeff,p,poly,basis.points)
#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)
# Generate Nr random points in pyramid (well, Nr-5 random plus vertices 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(0,0,-1); uvwr.set(0,1,-1); uvwr.set(0,2,0)
uvwr.set(1,0,-1); uvwr.set(1,1,1); uvwr.set(1,2,0)
uvwr.set(2,0,1); uvwr.set(2,1,1); uvwr.set(2,2,0)
uvwr.set(3,0,1); uvwr.set(3,1,-1); uvwr.set(3,2,0)
uvwr.set(4,0,0); uvwr.set(4,1,0); uvwr.set(4,2,1)
for i in range(5,Nr):
uvwr.set(i,0,1000)
uvwr.set(i,1,1000)
uvwr.set(i,2,1000)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment