diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 36bbca031974586a96d7e87ba858bf76c62b6101..836bff3b3595ff38e601f70904edd91f7e8efeb0 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -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) {
diff --git a/Numeric/BergotBasis.cpp b/Numeric/BergotBasis.cpp
index 8611c1e0a7835513ec8b7fff6e2f7d55500b6f91..52ac4fa1d93ef9ccf2fbaec1840edf36f158a75b 100644
--- a/Numeric/BergotBasis.cpp
+++ b/Numeric/BergotBasis.cpp
@@ -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,12 +98,9 @@ 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]));
- }
+ JacobiPolynomials jacobi(2.*mIJ+2.,0.,kMax);
+ jacobi.f(what,&(wf[0]));
+ jacobi.df(what,&(wg[0]));
}
// Recombine to find the shape function gradients
@@ -177,7 +147,7 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const
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]);
+ + uFcts[i]*vFcts[j]*powM1*(2.*oMW*wg[k]-mIJ*wf[k]);
}
}
}
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index e21b2a1736e189c95ad3d29dffa43c3823ee5184..ab0b938ec0bdf6a8d56668d65f325fa823a3f852 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -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);
diff --git a/Numeric/jacobiPolynomials.cpp b/Numeric/jacobiPolynomials.cpp
index 7a03c9c9b1f7d37ab1f1ac8bc8b05aa0312d5edf..a8693fe5cedad3b9545ccd370e7217cbe9d494e6 100644
--- a/Numeric/jacobiPolynomials.cpp
+++ b/Numeric/jacobiPolynomials.cpp
@@ -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;
}
+
}
diff --git a/Numeric/legendrePolynomials.cpp b/Numeric/legendrePolynomials.cpp
index b0cc7281c14b5b3ed8bbab47c1925d395bfee216..e5cc0cf48be4c63da9d6d5cfd04dbe02ffe63b77 100644
--- a/Numeric/legendrePolynomials.cpp
+++ b/Numeric/legendrePolynomials.cpp
@@ -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;
}
+
}
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index da1bc0d723db5e373bfaf0d334dec90da7b45d28..4584bba32e9d8f7565361f0fabd249b974635898 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -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");};
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index bda418798d7236708948fa4af0c03466bbb250ff..9334222f53aa690a80db01cf9e42737905645849 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -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);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 4444efbb899053ac4368494879d8fc25e20f0f17..0222972faa4c76414f74bfd7a60ed8f813ebed9d 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -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;
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 4a67e1f0f7f8fc86376c85c812b1adf21ba6cb68..c8b9f8e7537c652bca5e503e35a3dda1bb7e5d6d 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -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();
diff --git a/Numeric/pyramidalBasis.h b/Numeric/pyramidalBasis.h
index 5e7e33283695b283580bec06941c8c60accfbc89..a4cf260705392c2e6f6c47816caf0a5362fec0eb 100644
--- a/Numeric/pyramidalBasis.h
+++ b/Numeric/pyramidalBasis.h
@@ -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;
};
diff --git a/benchmarks/misc/testPyramidalBasis.py b/benchmarks/misc/testPyramidalBasis.py
index b2c3b1c45e2abee29f1dc12cdba17e02da1f822b..417a5ee2d70a311e1a546e93db884ad0456812a8 100644
--- a/benchmarks/misc/testPyramidalBasis.py
+++ b/benchmarks/misc/testPyramidalBasis.py
@@ -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)