From 18cf20b19dfa20d558b01d9322e985a10b21f914 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Tue, 11 Jun 2013 08:58:22 +0000
Subject: [PATCH] Validate new monomials/points generation algorithms
---
Common/Gmsh.cpp | 37 ++++++++++++++++++++++++++
Common/GmshDefines.h | 1 +
Geo/MElement.cpp | 22 +++++++++++++++-
Numeric/nodalBasis.cpp | 39 ++++++++++++++++++++++++++--
Numeric/nodalBasis.h | 18 +++++++------
Numeric/pointsGenerators.cpp | 14 ++++++----
Numeric/pointsGenerators.h | 3 +++
Numeric/polynomialBasis.cpp | 50 ++++++++++++++++++++++++++++++++++--
Numeric/polynomialBasis.h | 12 +++++----
9 files changed, 173 insertions(+), 23 deletions(-)
diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index c88d471b28..45e19235bc 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -22,6 +22,9 @@
#include "PView.h"
#endif
+//test new algo generation points
+#include "BasisFactory.h"
+
#if defined(HAVE_ONELAB)
#include "gmshLocalNetworkClient.h"
#endif
@@ -241,6 +244,40 @@ int GmshBatch()
solver_batch_cb(0, (void*)CTX::instance()->launchSolverAtStartup);
#endif
+ if (false) {
+ // 11/06/13 : This will be removed later !
+ for (int i = 1; i < MSH_NUM_TYPE+1; ++i) {
+ if (i == 76 || i == 77 || i == 78)
+ continue;
+ if (i == 67 || i == 68 || i == 70) {
+ Msg::Warning("ignoring unknown %d", i);
+ continue;
+ }
+ if (MElement::ParentTypeFromTag(i) == TYPE_PRI) {
+ Msg::Info("ignoring prism %d", i);
+ continue;
+ }
+ if (MElement::ParentTypeFromTag(i) == TYPE_PYR) {
+ Msg::Info("ignoring pyramid %d", i);
+ continue;
+ }
+ if (MElement::ParentTypeFromTag(i) == TYPE_POLYG) {
+ Msg::Info("ignoring polygone %d", i);
+ continue;
+ }
+ if (MElement::ParentTypeFromTag(i) == TYPE_POLYH) {
+ Msg::Info("ignoring polyhŹdre %d", i);
+ continue;
+ }
+ if (MElement::ParentTypeFromTag(i) == TYPE_XFEM) {
+ Msg::Info("ignoring xfem %d", i);
+ continue;
+ }
+ const nodalBasis *fs = BasisFactory::getNodalBasis(i);
+ if (fs) fs->compareNewAlgoPointsWithOld();
+ }
+ }
+
time_t now;
time(&now);
std::string currtime = ctime(&now);
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 66f24c25f5..a63733cce8 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -59,6 +59,7 @@
#define TYPE_HEX 8
#define TYPE_POLYG 9
#define TYPE_POLYH 10
+#define TYPE_XFEM 11
// Element types in .msh file format (numbers should not be changed)
#define MSH_LIN_2 1
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 9fa590258b..5f12e1aac5 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1167,6 +1167,7 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
{
switch(typeMSH){
case MSH_PNT : if(name) *name = "Point"; return 1;
+ case MSH_LIN_1 : if(name) *name = "Line 1"; return 1;
case MSH_LIN_2 : if(name) *name = "Line 2"; return 2;
case MSH_LIN_3 : if(name) *name = "Line 3"; return 2 + 1;
case MSH_LIN_4 : if(name) *name = "Line 4"; return 2 + 2;
@@ -1179,6 +1180,7 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_LIN_11 : if(name) *name = "Line 11"; return 2 + 9;
case MSH_LIN_B : if(name) *name = "Line Border"; return 2;
case MSH_LIN_C : if(name) *name = "Line Child"; return 2;
+ case MSH_TRI_1 : if(name) *name = "Triangle 1"; return 1;
case MSH_TRI_3 : if(name) *name = "Triangle 3"; return 3;
case MSH_TRI_6 : if(name) *name = "Triangle 6"; return 3 + 3;
case MSH_TRI_9 : if(name) *name = "Triangle 9"; return 3 + 6;
@@ -1198,6 +1200,7 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_TRI_27 : if(name) *name = "Triangle 27"; return 3 + 24;
case MSH_TRI_30 : if(name) *name = "Triangle 30"; return 3 + 27;
case MSH_TRI_B : if(name) *name = "Triangle Border"; return 3;
+ case MSH_QUA_1 : if(name) *name = "Quadrilateral 1"; return 1;
case MSH_QUA_4 : if(name) *name = "Quadrilateral 4"; return 4;
case MSH_QUA_8 : if(name) *name = "Quadrilateral 8"; return 4 + 4;
case MSH_QUA_9 : if(name) *name = "Quadrilateral 9"; return 9;
@@ -1212,8 +1215,14 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_QUA_12 : if(name) *name = "Quadrilateral 12"; return 12;
case MSH_QUA_16I : if(name) *name = "Quadrilateral 16I";return 16;
case MSH_QUA_20 : if(name) *name = "Quadrilateral 20"; return 20;
+ case MSH_QUA_24 : if(name) *name = "Quadrilateral 24"; return 24;
+ case MSH_QUA_28 : if(name) *name = "Quadrilateral 28"; return 28;
+ case MSH_QUA_32 : if(name) *name = "Quadrilateral 32"; return 32;
+ case MSH_QUA_36I : if(name) *name = "Quadrilateral 36I";return 36;
+ case MSH_QUA_40 : if(name) *name = "Quadrilateral 40"; return 40;
case MSH_POLYG_ : if(name) *name = "Polygon"; return 0;
case MSH_POLYG_B : if(name) *name = "Polygon Border"; return 0;
+ case MSH_TET_1 : if(name) *name = "Tetrahedron 1"; return 1;
case MSH_TET_4 : if(name) *name = "Tetrahedron 4"; return 4;
case MSH_TET_10 : if(name) *name = "Tetrahedron 10"; return 4 + 6;
case MSH_TET_20 : if(name) *name = "Tetrahedron 20"; return 4 + 12 + 4;
@@ -1226,6 +1235,12 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_TET_165 : if(name) *name = "Tetrahedron 165"; return (9*10*11)/6;
case MSH_TET_220 : if(name) *name = "Tetrahedron 220"; return (10*11*12)/6;
case MSH_TET_286 : if(name) *name = "Tetrahedron 286"; return (11*12*13)/6;
+ case MSH_TET_74 : if(name) *name = "Tetrahedron 74"; return 74;
+ case MSH_TET_100 : if(name) *name = "Tetrahedron 100"; return 100;
+ case MSH_TET_130 : if(name) *name = "Tetrahedron 130"; return 130;
+ case MSH_TET_164 : if(name) *name = "Tetrahedron 164"; return 164;
+ case MSH_TET_202 : if(name) *name = "Tetrahedron 202"; return 202;
+ case MSH_HEX_1 : if(name) *name = "Hexahedron 1"; return 1;
case MSH_HEX_8 : if(name) *name = "Hexahedron 8"; return 8;
case MSH_HEX_20 : if(name) *name = "Hexahedron 20"; return 8 + 12;
case MSH_HEX_27 : if(name) *name = "Hexahedron 27"; return 8 + 12 + 6 + 1;
@@ -1243,9 +1258,11 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_HEX_296 : if(name) *name = "Hexahedron 296"; return 296;
case MSH_HEX_386 : if(name) *name = "Hexahedron 386"; return 386;
case MSH_HEX_488 : if(name) *name = "Hexahedron 488"; return 488;
+ case MSH_PRI_1 : if(name) *name = "Prism 1"; return 1;
case MSH_PRI_6 : if(name) *name = "Prism 6"; return 6;
case MSH_PRI_15 : if(name) *name = "Prism 15"; return 6 + 9;
case MSH_PRI_18 : if(name) *name = "Prism 18"; return 6 + 9 + 3;
+ case MSH_PYR_1 : if(name) *name = "Pyramid 1"; return 1;
case MSH_PYR_5 : if(name) *name = "Pyramid 5"; return 5;
case MSH_PYR_13 : if(name) *name = "Pyramid 13"; return 5 + 8;
case MSH_PYR_14 : if(name) *name = "Pyramid 14"; return 5 + 8 + 1;
@@ -1424,10 +1441,13 @@ int MElement::ParentTypeFromTag(int tag)
case(MSH_HEX_218): case(MSH_HEX_296):
case(MSH_HEX_386): case(MSH_HEX_488):
return TYPE_HEX;
- case(MSH_POLYG_): case(MSH_POLYG_B):
+ case(MSH_POLYG_): case(MSH_POLYG_B):
return TYPE_POLYG;
case(MSH_POLYH_):
return TYPE_POLYH;
+ case(MSH_PNT_SUB): case(MSH_LIN_SUB):
+ case(MSH_TRI_SUB): case(MSH_TET_SUB):
+ return TYPE_XFEM;
default:
Msg::Error("Unknown type %i, assuming tetrahedron.", tag);
return TYPE_TET;
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 8a2c1b3c65..bd3411536c 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -4,11 +4,46 @@
// bugs and problems to the public mailing list <gmsh@geuz.org>.
#include <limits>
+#include <cmath>
#include "nodalBasis.h"
#include "BasisFactory.h"
//#include "pointsGenerators.h"
+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);
+ return 1;
+ }
+ if (points_newAlgo.size1() != points.size1()) {
+ Msg::Error("%d: pas meme taille (%d, %d, %d) %s", type, parentType, order, serendip, *name);
+ return 2;
+ }
+ for (int i = 0; i < points.size1(); ++i) {
+ for (int j = 0; j < points.size2(); ++j) {
+ //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);
+ 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)
+ Msg::Info("(%d,%d) : points %f / %f newPoints (mon %d)", i, j, points(i, j), points_newAlgo(i, j), monomials_newAlgo(i, j));
+ else
+ Msg::Error("(%d,%d) : points %f / %f newPoints (mon %d)", i, j, points(i, j), points_newAlgo(i, j), monomials_newAlgo(i, j));
+ }
+ Msg::Info(" ");
+ }
+ return 3;
+ }
+ }
+ }
+ Msg::Info("%d: ok (%d, %d, %d) %s", type, parentType, order, serendip, *name);
+ return 0;
+}
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
@@ -1373,7 +1408,7 @@ static void generateClosureOrder0(nodalBasis::clCont &closure, int nb)
nodalBasis::nodalBasis(int tag)
{
type = tag;
-
+
switch (tag) {
case MSH_PNT : parentType = TYPE_PNT; order = 0; serendip = false; break;
case MSH_LIN_1 : parentType = TYPE_LIN; order = 0; serendip = false; break;
@@ -1472,7 +1507,7 @@ nodalBasis::nodalBasis(int tag)
case MSH_HEX_512 : parentType = TYPE_HEX; order = 7; serendip = false; break;
case MSH_HEX_729 : parentType = TYPE_HEX; order = 8; serendip = false; break;
case MSH_HEX_1000: parentType = TYPE_HEX; order = 9; serendip = false; break;
- case MSH_HEX_20 : parentType = TYPE_HEX; order = 2; serendip = false; break;
+ case MSH_HEX_20 : parentType = TYPE_HEX; order = 2; serendip = true; break;
case MSH_HEX_56 : parentType = TYPE_HEX; order = 3; serendip = true; break;
case MSH_HEX_98 : parentType = TYPE_HEX; order = 4; serendip = true; break;
case MSH_HEX_152 : parentType = TYPE_HEX; order = 5; serendip = true; break;
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index eb545bbeaa..1957fdc3bb 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -14,12 +14,14 @@ class nodalBasis {
int type, parentType, order, dimension, numFaces;
bool serendip;
fullMatrix<double> points;
-
+ fullMatrix<double> points_newAlgo;
+ fullMatrix<int> monomials_newAlgo;
+
nodalBasis(int tag);
virtual ~nodalBasis() {}
-
+
virtual int getNumShapeFunctions() const = 0;
-
+
// Basis functions & gradients evaluation
virtual void f(double u, double v, double w, double *sf) const = 0;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const = 0;
@@ -27,9 +29,9 @@ class nodalBasis {
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const = 0;
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");};
-
-
-
+
+ int compareNewAlgoPointsWithOld() 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
// nodes of the element that rotates the element so that the
@@ -45,7 +47,7 @@ class nodalBasis {
typedef std::vector<closure> clCont;
clCont closures, fullClosures;
std::vector<int> closureRef;
-
+
// for a given face/edge, with both a sign and a rotation, give an
// ordered list of nodes on this face/edge
virtual int getClosureType(int id) const { return closures[id].type; }
@@ -53,7 +55,7 @@ class nodalBasis {
virtual const std::vector<int> &getFullClosure(int id) const { return fullClosures[id]; }
inline int getClosureId(int iFace, int iSign=1, int iRot=0) const;
inline void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const;
-
+
static inline int getTag(int parentTag, int order, bool serendip = false);
};
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 5b4e4221ac..bb50681b91 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -784,9 +784,6 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
return points;
}
-
-
-
fullMatrix<int> gmshGenerateMonomialsLine(int order)
{
fullMatrix<int> monomials(order + 1, 1);
@@ -797,6 +794,7 @@ fullMatrix<int> gmshGenerateMonomialsLine(int order)
}
return monomials;
}
+
fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
{
int nbMonomials = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
@@ -826,13 +824,14 @@ fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
monomials(index, 1) = monomials(i0, 1) + u_1 * i;
}
}
- fullMatrix<int> inner = gmshGenerateMonomialsQuadrangle(order-2);
+ fullMatrix<int> inner = gmshGenerateMonomialsTriangle(order-3);
inner.add(1);
monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
}
}
return monomials;
}
+
fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order)
{
int nbMonomials = (order+1)*(order+1);
@@ -872,6 +871,7 @@ fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order)
}
return monomials;
}
+
//KH : caveat : node coordinates are not yet coherent with node numbering associated
// to numbering of principal vertices of face !!!!
fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
@@ -922,6 +922,7 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
if (order > 2) {
fullMatrix<int> dudv = gmshGenerateMonomialsTriangle(order - 3);
+ dudv.add(1);
for (int iface = 0; iface < 4; ++iface) {
int i0 = MTetrahedron::faces_tetra(iface, 0);
@@ -954,6 +955,7 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
}
return monomials;
}
+
/*
fullMatrix<int> gmshGenerateMonomialsPrism(int order)
{
@@ -1003,6 +1005,7 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order)
}
*/
+
fullMatrix<int> gmshGenerateMonomialsHexahedron(int order)
{
int nbMonomials = (order+1)*(order+1)*(order+1);
@@ -1058,7 +1061,8 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order)
}
}
- fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 3);
+ fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 2);
+ dudv.add(1);
for (int iface = 0; iface < 6; ++iface) {
int i0 = MHexahedron::faces_hexa(iface, 0);
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index 5aa57b2c7f..8199e874f0 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -17,6 +17,7 @@
*/
// Points
+
fullMatrix<double> gmshGeneratePointsLine(int order);
fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip);
@@ -28,7 +29,9 @@ fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip);
fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip);
+
// Monomial exponents
+
fullMatrix<int> gmshGenerateMonomialsLine(int order);
fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip = false);
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index d5c06f6c06..bdecd24afd 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -13,6 +13,7 @@
#include "GmshMessage.h"
#include "polynomialBasis.h"
#include "GaussIntegration.h"
+#include "pointsGenerators.h"
#include <limits>
/*
@@ -342,7 +343,15 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
return coefficient;
}
-
+static void copy(const fullMatrix<int> &orig, fullMatrix<double> &b)
+{
+ b.resize(orig.size1(), orig.size2());
+ for (int i = 0; i < orig.size1(); ++i) {
+ for (int j = 0; j < orig.size2(); ++j) {
+ b(i, j) = static_cast<double>(orig(i, j));
+ }
+ }
+}
polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
{
@@ -376,6 +385,44 @@ polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
coefficients = generateLagrangeMonomialCoefficients(monomials, points);
+
+ // TEST NEW ALGO POINTS / MONOMIALS
+
+ bool centered = false;
+ switch (parentType) {
+ case TYPE_PNT :
+ monomials_newAlgo = gmshGenerateMonomialsLine(0);
+ break;
+ case TYPE_LIN :
+ monomials_newAlgo = gmshGenerateMonomialsLine(order);
+ centered = true;
+ break;
+ case TYPE_TRI :
+ monomials_newAlgo = gmshGenerateMonomialsTriangle(order, serendip);
+ break;
+ case TYPE_QUA :
+ if (serendip) return;
+ monomials_newAlgo = gmshGenerateMonomialsQuadrangle(order);
+ centered = true;
+ break;
+ case TYPE_TET :
+ monomials_newAlgo = gmshGenerateMonomialsTetrahedron(order, serendip);
+ break;
+ case TYPE_HEX :
+ if (serendip) return;
+ monomials_newAlgo = gmshGenerateMonomialsHexahedron(order);
+ centered = true;
+ break;
+ }
+ copy(monomials_newAlgo, points_newAlgo);
+ if (order == 0) return;
+ if (centered) {
+ points_newAlgo.scale(2./order);
+ points_newAlgo.add(-1.);
+ }
+ else {
+ points_newAlgo.scale(1./order);
+ }
}
@@ -386,7 +433,6 @@ polynomialBasis::~polynomialBasis()
-int polynomialBasis::getNumShapeFunctions() const { return coefficients.size1(); }
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 395f55ada3..fc04fc62d5 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -71,19 +71,21 @@ class polynomialBasis : public nodalBasis
// basis, we use the type of the corresponding gmsh element as type
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
-
+
polynomialBasis(int tag);
~polynomialBasis();
-
- virtual int getNumShapeFunctions() const;
-
+
+ int compareNewAlgoPointsWithOld() const;
+
+ virtual inline int getNumShapeFunctions() const {return coefficients.size1();}
+
virtual void f(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;
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;
-
+
inline void evaluateMonomials(double u, double v, double w, double p[]) const {
for (int j = 0; j < monomials.size1(); j++) {
p[j] = pow_int(u, (int)monomials(j, 0));
--
GitLab