diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index 251674c8d4581d33c9e005a0cfdb1d1c5d58cee9..c17b913e65caf5c6d276fdaf934bb1c5c4a642d1 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -275,8 +275,12 @@ int GmshBatch()
Msg::Warning("ignoring xfem %d", i);
continue;
}
+
const nodalBasis *fs = BasisFactory::getNodalBasis(i);
- if (fs) fs->compareNewAlgoPointsWithOld();
+
+ if (fs && !fs->compareNewAlgoPointsWithOld() && !fs->serendip) {
+ fs->compareNewAlgoBaseFunctionsWithOld();
+ }
}
}
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 8fb9be5916ba7239fdb1e3ace21f95734f152546..5856ae2a099b5a7b357546aeac7c8534edd43ed1 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -15,11 +15,11 @@ int nodalBasis::compareNewAlgoPointsWithOld() const
const char **name = new const char*[1];
MElement::getInfoMSH(type, name);
if (points_newAlgo.size1() == 0) {
- Msg::Warning("%d: pas de points (%d, %d, %d) %s", type, parentType, order, serendip, *name);
+ Msg::Warning("%d: pas de points (%d, %d, %d) %s", type, parentType, serendip, order, *name);
return 1;
}
if (points_newAlgo.size1() != points.size1()) {
- Msg::Error("%d: pas meme taille (%d, %d, %d) %s", type, parentType, order, serendip, *name);
+ Msg::Error("%d: pas meme taille (%d, %d, %d) %s", type, parentType, serendip, order, *name);
Msg::Error(" |--> points[%d, %d] vs newPoints[%d, %d]", points.size1(), points.size2(), points_newAlgo.size1(), points_newAlgo.size2());
return 2;
}
@@ -28,7 +28,7 @@ int nodalBasis::compareNewAlgoPointsWithOld() const
//Msg::Info("(i, j) = (%d, %d)", i, j);
//Msg::Info(" ");
if (std::abs(points(i, j) - points_newAlgo(i, j)) > 1e-15) {
- Msg::Error("%d: correspond pas (%d, %d, %d) %s", type, parentType, order, serendip, *name);
+ Msg::Error("%d: correspond pas (%d, %d, %d) %s", type, parentType, serendip, order, *name);
for (int i = 0; i < points.size1(); ++i) {
for (int j = 0; j < points.size2(); ++j) {
if(std::abs(points(i, j) - points_newAlgo(i, j)) <= 1e-15)
@@ -42,7 +42,73 @@ int nodalBasis::compareNewAlgoPointsWithOld() const
}
}
}
- Msg::Info("%d: ok (%d, %d, %d) %s", type, parentType, order, serendip, *name);
+ Msg::Info("%d: ok (%d, %d, %d) %s", type, parentType, serendip, order, *name);
+ return 0;
+}
+
+int nodalBasis::compareNewAlgoBaseFunctionsWithOld() const
+{
+ const char **name = new const char*[1];
+ MElement::getInfoMSH(type, name);
+ int ndof = points.size1();
+ int P[3];
+ double oldVal[ndof], newVal[ndof];
+ for (int i = 0; i < 10; ++i) {
+ if (i == 0) {
+ P[0] = 0;
+ P[1] = 0;
+ P[2] = 0;
+ }
+ else if (i == 1) {
+ P[0] = 0;
+ P[1] = 0;
+ P[2] = 1000;
+ }
+ else if (i == 2) {
+ P[0] = 1000;
+ P[1] = 0;
+ P[2] = 0;
+ }
+ else {
+ P[0] = std::rand() % 1001;
+ P[1] = std::rand() % (1001 - P[0]);
+ P[2] = std::rand() % (1001 - P[0] - P[1]);
+ }
+
+ f(P[0] / 1000., P[1] / 1000., P[2] / 1000., oldVal);
+ fnew(P[0] / 1000., P[1] / 1000., P[2] / 1000., newVal);
+
+ double sumError = .0;
+ double sumOne[2];
+ sumOne[0] = .0;
+ sumOne[1] = .0;
+ for (int j = 0; j < ndof; ++j) {
+ double diff = std::abs(oldVal[j] - newVal[j]);
+ //double mean = std::sqrt(.5 * (oldVal[j]*oldVal[j] + newVal[j]*newVal[j]));
+ //double error = diff/(mean + 1e-15);
+ double error = diff;
+ sumError += error*error;
+ sumOne[0] += oldVal[j];
+ sumOne[1] += newVal[j];
+ }
+ sumError = std::sqrt(sumError/ndof);
+ /*if (sumError > 1e-4) {
+ Msg::Error("(%.2f, %.2f, %.2f) -> fold=%.2e / fnew=%.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumOne[0], sumOne[1]);
+ //for (int j = 0; j < ndof; ++j) {
+ // Msg::Info("old %.3e vs %.3e new", oldVal[j], newVal[j]);
+ //}
+ }*/
+ if (sumError > 1e-13) {
+ if (type < 10) Msg::Warning(" f(%.2f, %.2f, %.2f) bad precision diff %.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
+ else if (type < 100) Msg::Warning(" f(%.2f, %.2f, %.2f) bad precision diff %.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
+ else Msg::Warning(" f(%.2f, %.2f, %.2f) bad precision diff %.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
+ return 1;
+ }
+ /*else if (type < 10) Msg::Info(" f(%.2f, %.2f, %.2f) ok", P[0] / 1000., P[1] / 1000., P[2] / 1000.);
+ else if (type < 100) Msg::Info(" f(%.2f, %.2f, %.2f) ok", P[0] / 1000., P[1] / 1000., P[2] / 1000.);
+ else Msg::Info(" f(%.2f, %.2f, %.2f) ok", P[0] / 1000., P[1] / 1000., P[2] / 1000.);
+ */
+ }
return 0;
}
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index 1957fdc3bb9211dc39d3a23206c30b74851c5553..c7339acca791daeb6f1f5bcdfeadc80bb98d0009 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -16,6 +16,7 @@ class nodalBasis {
fullMatrix<double> points;
fullMatrix<double> points_newAlgo;
fullMatrix<int> monomials_newAlgo;
+ fullMatrix<double> coefficients_newAlgo;
nodalBasis(int tag);
virtual ~nodalBasis() {}
@@ -24,6 +25,7 @@ class nodalBasis {
// Basis functions & gradients evaluation
virtual void f(double u, double v, double w, double *sf) const = 0;
+ virtual void fnew(double u, double v, double w, double *sf) const = 0;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const = 0;
virtual void df(double u, double v, double w, double grads[][3]) const = 0;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const = 0;
@@ -31,6 +33,7 @@ class nodalBasis {
virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");};
int compareNewAlgoPointsWithOld() const;
+ int compareNewAlgoBaseFunctionsWithOld() const;
// closures is the list of the nodes of each face, in the order of
// the polynomialBasis of the face; fullClosures is mapping of the
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 4d12d058af3f39c868b4564399e5abd6ef369840..2b6365b1273740875c05df2e13783e783449b4ea 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -826,9 +826,11 @@ fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
monomials(index, 1) = monomials(i0, 1) + u_1 * i;
}
}
- fullMatrix<int> inner = gmshGenerateMonomialsTriangle(order-3);
- inner.add(1);
- monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
+ if (!serendip) {
+ fullMatrix<int> inner = gmshGenerateMonomialsTriangle(order-3);
+ inner.add(1);
+ monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
+ }
}
}
return monomials;
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index d4cef22282809c3f54f24709ec9ba66c3f8f83af..27717851935b57674a330ea10ddfe4d78770e484 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -9,6 +9,7 @@
//
#include <stdlib.h>
+#include <cmath>
#include "GmshDefines.h"
#include "GmshMessage.h"
#include "polynomialBasis.h"
@@ -137,6 +138,87 @@ static fullMatrix<double> generatePascalQuad(int order)
static fullMatrix<double> generatePascalHex(int order, bool serendip)
{
+ if (false && serendip) {
+ fullMatrix<double> monomials( 8+(order-1)*12, 3);
+ monomials(0,0)=0.;
+ monomials(0,1)=0.;
+ monomials(0,2)=0.;
+ monomials(1,0)=1.;
+ monomials(1,1)=0.;
+ monomials(1,2)=0.;
+ monomials(2,0)=0.;
+ monomials(2,1)=1.;
+ monomials(2,2)=0.;
+ monomials(3,0)=1.;
+ monomials(3,1)=1.;
+ monomials(3,2)=0.;
+
+ monomials(4,0)=0.;
+ monomials(4,1)=0.;
+ monomials(4,2)=1.;
+ monomials(5,0)=1.;
+ monomials(5,1)=0.;
+ monomials(5,2)=1.;
+ monomials(6,0)=0.;
+ monomials(6,1)=1.;
+ monomials(6,2)=1.;
+ monomials(7,0)=1.;
+ monomials(7,1)=1.;
+ monomials(7,2)=1.;
+ int index = 8;
+ for (int p = 2; p <= order; p++) {
+ monomials(index, 0) = p;
+ monomials(index, 1) = 0;
+ monomials(index, 2) = 0;
+ index++;
+ monomials(index, 0) = p;
+ monomials(index, 1) = 1;
+ monomials(index, 2) = 0;
+ index++;
+ monomials(index, 0) = p;
+ monomials(index, 1) = 1;
+ monomials(index, 2) = 1;
+ index++;
+ monomials(index, 0) = p;
+ monomials(index, 1) = 0;
+ monomials(index, 2) = 1;
+ index++;
+ monomials(index, 0) = 0;
+ monomials(index, 1) = p;
+ monomials(index, 2) = 0;
+ index++;
+ monomials(index, 0) = 1;
+ monomials(index, 1) = p;
+ monomials(index, 2) = 0;
+ index++;
+ monomials(index, 0) = 1;
+ monomials(index, 1) = p;
+ monomials(index, 2) = 1;
+ index++;
+ monomials(index, 0) = 0;
+ monomials(index, 1) = p;
+ monomials(index, 2) = 1;
+ index++;
+ monomials(index, 0) = 0;
+ monomials(index, 1) = 0;
+ monomials(index, 2) = p;
+ index++;
+ monomials(index, 0) = 1;
+ monomials(index, 1) = 0;
+ monomials(index, 2) = p;
+ index++;
+ monomials(index, 0) = 1;
+ monomials(index, 1) = 1;
+ monomials(index, 2) = p;
+ index++;
+ monomials(index, 0) = 0;
+ monomials(index, 1) = 1;
+ monomials(index, 2) = p;
+ index++;
+ }
+ return monomials;
+ }
+
int siz = (order+1)*(order+1)*(order+1);
if (serendip) siz -= (order-1)*(order-1)*(order-1);
fullMatrix<double> monomials( siz, 3);
@@ -340,6 +422,19 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
fullMatrix<double> coefficient(ndofs, ndofs);
Vandermonde.invert(coefficient);
+
+ fullMatrix<double> unity(ndofs, ndofs);
+ Vandermonde.mult(coefficient, unity);
+ double max = .0;
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ if (i == j) unity(i, j) -= 1.;
+ //Msg::Info(" unity(%d, %d) = %.3e", i, j, unity(i, j));
+ max = std::max(max, std::abs(unity(i, j)));
+ }
+ }
+ if (max > 1e-10) Msg::Info(" max unity = %.3e", max);
+
return coefficient;
}
@@ -417,7 +512,7 @@ polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
break;
}
copy(monomials_newAlgo, points_newAlgo);
- if (order == 0) return;
+ //if (order == 0) return;
switch (rescale) {
case 0 :
points_newAlgo.scale(1./order);
@@ -440,6 +535,23 @@ polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
default :
break;
}
+
+ fullMatrix<double> monDouble;
+ switch (parentType) {
+ case TYPE_PNT :
+ case TYPE_LIN :
+ case TYPE_TRI :
+ case TYPE_TET :
+ copy(monomials_newAlgo, monDouble);
+ break;
+ case TYPE_QUA :
+ case TYPE_HEX :
+ case TYPE_PRI :
+ //if (serendip) monDouble = monomials;
+ /*else*/ copy(monomials_newAlgo, monDouble);
+ break;
+ }
+ coefficients_newAlgo = generateLagrangeMonomialCoefficients(monDouble, points_newAlgo);
}
@@ -451,6 +563,23 @@ polynomialBasis::~polynomialBasis()
+inline void polynomialBasis::evaluateMonomialsNew(double u, double v, double w, double p[]) const
+{
+ for (int j = 0; j < monomials_newAlgo.size1(); ++j) {
+ p[j] = 1.;
+ switch (dimension) {
+ case 3 :
+ p[j] *= pow_int(w, monomials_newAlgo(j, 2));
+ case 2 :
+ p[j] *= pow_int(v, monomials_newAlgo(j, 1));
+ case 1 :
+ p[j] *= pow_int(u, monomials_newAlgo(j, 0));
+ default :
+ break;
+ }
+ }
+}
+
void polynomialBasis::f(double u, double v, double w, double *sf) const
@@ -464,6 +593,17 @@ void polynomialBasis::f(double u, double v, double w, double *sf) const
}
}
}
+void polynomialBasis::fnew(double u, double v, double w, double *sf) const
+{
+ double p[1256];
+ evaluateMonomialsNew(u, v, w, p);
+ for (int i = 0; i < coefficients_newAlgo.size1(); i++) {
+ sf[i] = 0.0;
+ for (int j = 0; j < coefficients_newAlgo.size2(); j++) {
+ sf[i] += coefficients_newAlgo(i, j) * p[j];
+ }
+ }
+}
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 5c67c3580b3ef1f9548cbfed20be1bbefea358a1..da48cc29cfe9901a17ab42f1809e37ed6f49cf9f 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -78,6 +78,7 @@ class polynomialBasis : public nodalBasis
virtual inline int getNumShapeFunctions() const {return coefficients.size1();}
virtual void f(double u, double v, double w, double *sf) const;
+ virtual void fnew(double u, double v, double w, double *sf) 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;
@@ -91,6 +92,7 @@ class polynomialBasis : public nodalBasis
if (monomials.size2() > 2) p[j] *= pow_int(w, (int)monomials(j, 2));
}
}
+ inline void evaluateMonomialsNew(double u, double v, double w, double p[]) const;
};
#endif
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 13c030b81b6c5d5c9e432d4f3fe855f9c05324e7..888752d749e907a55c55fd94f7d54bd4cf7ddc2d 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -5,6 +5,7 @@
#include "pyramidalBasis.h"
#include "pointsGenerators.h"
+#include <cmath>
static void copy(const fullMatrix<int> &orig, fullMatrix<double> &b)
@@ -17,6 +18,55 @@ static void copy(const fullMatrix<int> &orig, fullMatrix<double> &b)
}
}
+static fullMatrix<double> generateLagrangeMonomialCoefficientsPyr
+ (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
+{
+ if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
+ Msg::Fatal("Wrong sizes for Lagrange coefficients generation %d %d -- %d %d",
+ monomial.size1(),point.size1(),
+ monomial.size2(),point.size2() );
+ return fullMatrix<double>(1, 1);
+ }
+
+ int ndofs = monomial.size1();
+ int dim = monomial.size2();
+ fullMatrix<double> ppoint(point.size1(), point.size2());
+
+ for (int i = 0; i < ndofs; ++i) {
+ ppoint(i, 2) = 1. - point(i, 2);
+ if (ppoint(i, 2) != .0) {
+ ppoint(i, 0) = point(i, 0) / ppoint(i, 2);
+ ppoint(i, 1) = point(i, 1) / ppoint(i, 2);
+ }
+ }
+
+ fullMatrix<double> Vandermonde(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double dd = 1.;
+ for (int k = 0; k < dim; k++) dd *= pow(ppoint(j, k), monomial(i, k));
+ Vandermonde(i, j) = dd;
+ }
+ }
+
+ fullMatrix<double> coefficient(ndofs, ndofs);
+ Vandermonde.invert(coefficient);
+
+ fullMatrix<double> unity(ndofs, ndofs);
+ Vandermonde.mult(coefficient, unity);
+ double max = .0;
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ if (i == j) unity(i, j) -= 1.;
+ //Msg::Info(" unity(%d, %d) = %.3e", i, j, unity(i, j));
+ max = std::max(max, std::abs(unity(i, j)));
+ }
+ }
+ if (max > 1e-10) Msg::Info("mon max unity = %.3e", max);
+
+ return coefficient;
+}
+
pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
{
@@ -38,6 +88,19 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
delete[] fval;
+ fullMatrix<double> unity(num_points, num_points);
+ VDM.mult(VDMinv, unity);
+ double max = .0;
+ for (int i = 0; i < num_points; i++) {
+ for (int j = 0; j < num_points; j++) {
+ if (i == j) unity(i, j) -= 1.;
+ //Msg::Info(" unity(%d, %d) = %.3e", i, j, unity(i, j));
+ max = std::max(max, std::abs(unity(i, j)));
+ }
+ }
+ if (max > 1e-10) Msg::Info("leg max unity = %.3e", max);
+
+
// TEST NEW ALGO POINTS / MONOMIAL
monomials_newAlgo = gmshGenerateMonomialsPyramid(order, serendip);
@@ -52,6 +115,10 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
points_newAlgo(i, 1) = duv + points_newAlgo(i, 1) * 2. / order;
}
}
+
+ fullMatrix<double> monDouble;
+ copy(monomials_newAlgo, monDouble);
+ coefficients_newAlgo = generateLagrangeMonomialCoefficientsPyr(monDouble, points_newAlgo);
}
@@ -83,6 +150,17 @@ void pyramidalBasis::f(double u, double v, double w, double *val) const
delete[] fval;
}
+void pyramidalBasis::fnew(double u, double v, double w, double *sf) const
+{
+ double p[1256];
+ evaluateMonomialsNew(u, v, w, p);
+ for (int i = 0; i < coefficients_newAlgo.size1(); i++) {
+ sf[i] = 0.0;
+ for (int j = 0; j < coefficients_newAlgo.size2(); j++) {
+ sf[i] += coefficients_newAlgo(i, j) * p[j];
+ }
+ }
+}
diff --git a/Numeric/pyramidalBasis.h b/Numeric/pyramidalBasis.h
index b570a9b1e1d311269de84021c117341ef230d1f9..998f49d4d757affa77b395ac0ae5fdc13f2b9f4a 100644
--- a/Numeric/pyramidalBasis.h
+++ b/Numeric/pyramidalBasis.h
@@ -27,12 +27,25 @@ class pyramidalBasis: public nodalBasis
~pyramidalBasis();
virtual void f(double u, double v, double w, double *val) const;
+ virtual void fnew(double u, double v, double w, double *val) const;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const;
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual int getNumShapeFunctions() const;
+ inline void evaluateMonomialsNew(double u, double v, double w, double p[]) const {
+ w = 1-w;
+ if (w != .0) {
+ u = u/w;
+ v = v/w;
+ }
+ for (int j = 0; j < monomials_newAlgo.size1(); j++) {
+ p[j] = pow_int(u, monomials_newAlgo(j, 0));
+ if (dimension > 1) p[j] *= pow_int(v, monomials_newAlgo(j, 1));
+ if (dimension > 2) p[j] *= pow_int(w, monomials_newAlgo(j, 2));
+ }
+ }
};