From ecab3ae465c4da774efd08d69d1ac44bd82acafd Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Thu, 28 Apr 2011 07:51:08 +0000
Subject: [PATCH] move bezierBasis outside of jacobianBasis so that it can be
used in other contexts
---
Common/gmshpy.i | 8 +-
Geo/MLine.cpp | 20 ++---
Geo/MPoint.h | 2 +-
Geo/MPrism.cpp | 8 +-
Geo/MQuadrangle.cpp | 40 +++++-----
Geo/MTetrahedron.cpp | 30 ++++----
Geo/MTriangle.cpp | 42 +++++-----
Mesh/qualityMeasures.cpp | 4 +-
Numeric/JacobianBasis.cpp | 57 ++++++++------
Numeric/JacobianBasis.h | 32 ++++----
Numeric/polynomialBasis.cpp | 143 ++++++++++++++++++++++-------------
Numeric/polynomialBasis.h | 1 +
Plugin/AnalyseCurvedMesh.cpp | 35 +++++----
Plugin/AnalyseCurvedMesh.h | 2 +-
14 files changed, 238 insertions(+), 186 deletions(-)
diff --git a/Common/gmshpy.i b/Common/gmshpy.i
index 670ea9d373..f4f9effa5b 100644
--- a/Common/gmshpy.i
+++ b/Common/gmshpy.i
@@ -1,5 +1,9 @@
%feature("autodoc", "1");
%module gmshpy
+%include std_string.i
+%include std_vector.i
+%include std_list.i
+
%{
#include "GmshConfig.h"
#include "GModel.h"
@@ -61,10 +65,6 @@
};
%}
-%include std_string.i
-%include std_vector.i
-%include std_list.i
-
%init %{
errorHandler *eH = new errorHandler;
Msg::SetCallback(eH);
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 96184bdce3..90efd2e595 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -35,16 +35,16 @@ const JacobianBasis* MLine::getJacobianFuncSpace(int o) const
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 1: return JacobianBases::find(MSH_LIN_2);
- case 2: return JacobianBases::find(MSH_LIN_3);
- case 3: return JacobianBases::find(MSH_LIN_4);
- case 4: return JacobianBases::find(MSH_LIN_5);
- case 5: return JacobianBases::find(MSH_LIN_6);
- case 6: return JacobianBases::find(MSH_LIN_7);
- case 7: return JacobianBases::find(MSH_LIN_8);
- case 8: return JacobianBases::find(MSH_LIN_9);
- case 9: return JacobianBases::find(MSH_LIN_10);
- case 10: return JacobianBases::find(MSH_LIN_11);
+ case 1: return JacobianBasis::find(MSH_LIN_2);
+ case 2: return JacobianBasis::find(MSH_LIN_3);
+ case 3: return JacobianBasis::find(MSH_LIN_4);
+ case 4: return JacobianBasis::find(MSH_LIN_5);
+ case 5: return JacobianBasis::find(MSH_LIN_6);
+ case 6: return JacobianBasis::find(MSH_LIN_7);
+ case 7: return JacobianBasis::find(MSH_LIN_8);
+ case 8: return JacobianBasis::find(MSH_LIN_9);
+ case 9: return JacobianBasis::find(MSH_LIN_10);
+ case 10: return JacobianBasis::find(MSH_LIN_11);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 01402a5553..efbf2844d4 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -60,7 +60,7 @@ class MPoint : public MElement {
}
virtual const JacobianBasis* getJacobianFuncSpace(int o) const
{
- return JacobianBases::find(MSH_PNT);
+ return JacobianBasis::find(MSH_PNT);
}
virtual bool isInside(double u, double v, double w)
{
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index e87e3948a0..adbfdb176e 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -63,15 +63,15 @@ const JacobianBasis* MPrism::getJacobianFuncSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return JacobianBases::find(MSH_PRI_6);
- case 2: return JacobianBases::find(MSH_PRI_18);
+ case 1: return JacobianBasis::find(MSH_PRI_6);
+ case 2: return JacobianBasis::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return JacobianBases::find(MSH_PRI_6);
- case 2: return JacobianBases::find(MSH_PRI_18);
+ case 1: return JacobianBasis::find(MSH_PRI_6);
+ case 2: return JacobianBasis::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 8a41e387bc..fe21cb7cd9 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -62,29 +62,29 @@ const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 1: return JacobianBases::find(MSH_QUA_4);
- case 2: return JacobianBases::find(MSH_QUA_8);
- case 3: return JacobianBases::find(MSH_QUA_12);
- case 4: return JacobianBases::find(MSH_QUA_16I);
- case 5: return JacobianBases::find(MSH_QUA_20);
- case 6: return JacobianBases::find(MSH_QUA_24);
- case 7: return JacobianBases::find(MSH_QUA_28);
- case 8: return JacobianBases::find(MSH_QUA_32);
- case 9: return JacobianBases::find(MSH_QUA_36I);
- case 10: return JacobianBases::find(MSH_QUA_40);
+ case 1: return JacobianBasis::find(MSH_QUA_4);
+ case 2: return JacobianBasis::find(MSH_QUA_8);
+ case 3: return JacobianBasis::find(MSH_QUA_12);
+ case 4: return JacobianBasis::find(MSH_QUA_16I);
+ case 5: return JacobianBasis::find(MSH_QUA_20);
+ case 6: return JacobianBasis::find(MSH_QUA_24);
+ case 7: return JacobianBasis::find(MSH_QUA_28);
+ case 8: return JacobianBasis::find(MSH_QUA_32);
+ case 9: return JacobianBasis::find(MSH_QUA_36I);
+ case 10: return JacobianBasis::find(MSH_QUA_40);
}
}
switch (order) {
- case 1: return JacobianBases::find(MSH_QUA_4);
- case 2: return JacobianBases::find(MSH_QUA_9);
- case 3: return JacobianBases::find(MSH_QUA_16);
- case 4: return JacobianBases::find(MSH_QUA_25);
- case 5: return JacobianBases::find(MSH_QUA_36);
- case 6: return JacobianBases::find(MSH_QUA_49);
- case 7: return JacobianBases::find(MSH_QUA_64);
- case 8: return JacobianBases::find(MSH_QUA_81);
- case 9: return JacobianBases::find(MSH_QUA_100);
- case 10: return JacobianBases::find(MSH_QUA_121);
+ case 1: return JacobianBasis::find(MSH_QUA_4);
+ case 2: return JacobianBasis::find(MSH_QUA_9);
+ case 3: return JacobianBasis::find(MSH_QUA_16);
+ case 4: return JacobianBasis::find(MSH_QUA_25);
+ case 5: return JacobianBasis::find(MSH_QUA_36);
+ case 6: return JacobianBasis::find(MSH_QUA_49);
+ case 7: return JacobianBasis::find(MSH_QUA_64);
+ case 8: return JacobianBasis::find(MSH_QUA_81);
+ case 9: return JacobianBasis::find(MSH_QUA_100);
+ case 10: return JacobianBasis::find(MSH_QUA_121);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
return 0;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 9b4efb8d83..11632f250e 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -166,26 +166,26 @@ const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return JacobianBases::find(MSH_TET_4);
- case 2: return JacobianBases::find(MSH_TET_10);
- case 3: return JacobianBases::find(MSH_TET_20);
- case 4: return JacobianBases::find(MSH_TET_34);
- case 5: return JacobianBases::find(MSH_TET_52);
+ case 1: return JacobianBasis::find(MSH_TET_4);
+ case 2: return JacobianBasis::find(MSH_TET_10);
+ case 3: return JacobianBasis::find(MSH_TET_20);
+ case 4: return JacobianBasis::find(MSH_TET_34);
+ case 5: return JacobianBasis::find(MSH_TET_52);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return JacobianBases::find(MSH_TET_4);
- case 2: return JacobianBases::find(MSH_TET_10);
- case 3: return JacobianBases::find(MSH_TET_20);
- case 4: return JacobianBases::find(MSH_TET_35);
- case 5: return JacobianBases::find(MSH_TET_56);
- case 6: return JacobianBases::find(MSH_TET_84);
- case 7: return JacobianBases::find(MSH_TET_120);
- case 8: return JacobianBases::find(MSH_TET_165);
- case 9: return JacobianBases::find(MSH_TET_220);
- case 10: return JacobianBases::find(MSH_TET_286);
+ case 1: return JacobianBasis::find(MSH_TET_4);
+ case 2: return JacobianBasis::find(MSH_TET_10);
+ case 3: return JacobianBasis::find(MSH_TET_20);
+ case 4: return JacobianBasis::find(MSH_TET_35);
+ case 5: return JacobianBasis::find(MSH_TET_56);
+ case 6: return JacobianBasis::find(MSH_TET_84);
+ case 7: return JacobianBasis::find(MSH_TET_120);
+ case 8: return JacobianBasis::find(MSH_TET_165);
+ case 9: return JacobianBasis::find(MSH_TET_220);
+ case 10: return JacobianBasis::find(MSH_TET_286);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index c1633e3ae6..2adf4d9787 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -160,32 +160,32 @@ const JacobianBasis* MTriangle::getJacobianFuncSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 0: return JacobianBases::find(MSH_TRI_1);
- case 1: return JacobianBases::find(MSH_TRI_3);
- case 2: return JacobianBases::find(MSH_TRI_6);
- case 3: return JacobianBases::find(MSH_TRI_9);
- case 4: return JacobianBases::find(MSH_TRI_12);
- case 5: return JacobianBases::find(MSH_TRI_15I);
- case 6: return JacobianBases::find(MSH_TRI_18);
- case 7: return JacobianBases::find(MSH_TRI_21I);
- case 8: return JacobianBases::find(MSH_TRI_24);
- case 9: return JacobianBases::find(MSH_TRI_27);
- case 10: return JacobianBases::find(MSH_TRI_30);
+ case 0: return JacobianBasis::find(MSH_TRI_1);
+ case 1: return JacobianBasis::find(MSH_TRI_3);
+ case 2: return JacobianBasis::find(MSH_TRI_6);
+ case 3: return JacobianBasis::find(MSH_TRI_9);
+ case 4: return JacobianBasis::find(MSH_TRI_12);
+ case 5: return JacobianBasis::find(MSH_TRI_15I);
+ case 6: return JacobianBasis::find(MSH_TRI_18);
+ case 7: return JacobianBasis::find(MSH_TRI_21I);
+ case 8: return JacobianBasis::find(MSH_TRI_24);
+ case 9: return JacobianBasis::find(MSH_TRI_27);
+ case 10: return JacobianBasis::find(MSH_TRI_30);
default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return JacobianBases::find(MSH_TRI_3);
- case 2: return JacobianBases::find(MSH_TRI_6);
- case 3: return JacobianBases::find(MSH_TRI_10);
- case 4: return JacobianBases::find(MSH_TRI_15);
- case 5: return JacobianBases::find(MSH_TRI_21);
- case 6: return JacobianBases::find(MSH_TRI_28);
- case 7: return JacobianBases::find(MSH_TRI_36);
- case 8: return JacobianBases::find(MSH_TRI_45);
- case 9: return JacobianBases::find(MSH_TRI_55);
- case 10: return JacobianBases::find(MSH_TRI_66);
+ case 1: return JacobianBasis::find(MSH_TRI_3);
+ case 2: return JacobianBasis::find(MSH_TRI_6);
+ case 3: return JacobianBasis::find(MSH_TRI_10);
+ case 4: return JacobianBasis::find(MSH_TRI_15);
+ case 5: return JacobianBasis::find(MSH_TRI_21);
+ case 6: return JacobianBasis::find(MSH_TRI_28);
+ case 7: return JacobianBasis::find(MSH_TRI_36);
+ case 8: return JacobianBasis::find(MSH_TRI_45);
+ case 9: return JacobianBasis::find(MSH_TRI_55);
+ case 10: return JacobianBasis::find(MSH_TRI_66);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index b17f53269c..c311c180b0 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -304,7 +304,7 @@ double mesh_functional_distorsion_p2_exact(MTriangle *t)
double mesh_functional_distorsion_pN(MElement *t)
{
- const JacobianBasis *jac = t->getJacobianFuncSpace();
+ const bezierBasis *jac = t->getJacobianFuncSpace()->bezier;
// jac->monomials.print("coucou");
fullVector<double>Ji(jac->points.size1());
for (int i=0;i<jac->points.size1();i++){
@@ -386,7 +386,7 @@ static double mesh_functional_distorsion(MTetrahedron *t, double u, double v, do
double qmDistorsionOfMapping(MTetrahedron *t)
{
- const JacobianBasis *jac = t->getJacobianFuncSpace();
+ const bezierBasis *jac = t->getJacobianFuncSpace()->bezier;
fullVector<double>Ji(jac->points.size1());
for (int i=0;i<jac->points.size1();i++){
const double u = jac->points(i,0);
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 76077e344f..0666b34c71 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -1187,19 +1187,17 @@ static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &poi
}
return;
}
+std::map<int, bezierBasis> bezierBasis::_bbs;
-std::map<int, JacobianBasis> JacobianBases::fs;
-
-const JacobianBasis *JacobianBases::find(int tag)
+const bezierBasis *bezierBasis::find(int tag)
{
- std::map<int, JacobianBasis>::const_iterator it = fs.find(tag);
- if (it != fs.end()) return &it->second;
-
- JacobianBasis B;
- B.numLagPts = -1;
+ std::map<int, bezierBasis>::const_iterator it = _bbs.find(tag);
+ if (it != _bbs.end())
+ return &it->second;
+ bezierBasis &B = _bbs[tag];
const polynomialBasis *F = polynomialBases::find(tag);
- int order = F->order;
+ int o = F->order;
switch (F->parentType) {
case TYPE_PNT :
B.numLagPts = 1;
@@ -1208,7 +1206,6 @@ const JacobianBasis *JacobianBases::find(int tag)
B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, 0, 0);
break;
case TYPE_LIN : {
- int o = order - 1;
B.numLagPts = 2;
B.exposants = generate1DExposants(o);
B.points = generate1DPoints(o);
@@ -1216,68 +1213,82 @@ const JacobianBasis *JacobianBases::find(int tag)
std::vector< fullMatrix<double> > subPoints = generateSubPointsLine(0);
B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 0);
const polynomialBasis *F = polynomialBases::find(tag);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
break;
}
case TYPE_TRI : {
- int o = 2 * (order - 1);
B.numLagPts = 3;
B.exposants = generateExposantsTriangle(o);
B.points = gmshGeneratePointsTriangle(o,false);
B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 2);
std::vector< fullMatrix<double> > subPoints = generateSubPointsTriangle(o);
B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 2);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
break;
}
case TYPE_TET : {
- int o = 3 * (order - 1);
B.numLagPts = 4;
B.exposants = generateExposantsTetrahedron(o);
B.points = gmshGeneratePointsTetrahedron(o,false);
B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 3);
std::vector< fullMatrix<double> > subPoints = generateSubPointsTetrahedron(o);
B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 3);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
break;
}
case TYPE_QUA : {
- int o = 2 * order - 1;
B.numLagPts = 4;
B.exposants = generateExposantsQuad(o);
B.points = generatePointsQuad(o,false);
B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 0);
std::vector< fullMatrix<double> > subPoints = generateSubPointsQuad(o);
B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 0);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
break;
}
case TYPE_PRI : {
- int o = order * 3 - 1; // o=3*ord-2 on triangle base and =3*ord-1 for third dimension
B.numLagPts = 4;
B.exposants = generateExposantsPrism(o);
B.points = generatePointsPrism(o,false);
B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, o, 2);
std::vector< fullMatrix<double> > subPoints = generateSubPointsPrism(o);
B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, o, 2);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
break;
}
default : {
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
B.numLagPts = 4;
B.exposants = generateExposantsTetrahedron(0);
- B.points = gmshGeneratePointsTetrahedron(0,false);
+ B.points = gmshGeneratePointsTetrahedron(0, false);
B.matrixLag2Bez = generateLag2BezMatrix(B.exposants, B.points, 0, 3);
std::vector< fullMatrix<double> > subPoints = generateSubPointsTetrahedron(0);
B.divisor = generateDivisor(B.exposants, subPoints, B.matrixLag2Bez, 0, 3);
F = polynomialBases::find(MSH_TET_4);
- generateGradShapes(B, B.points, F->monomials, F->coefficients);
}
}
B.numDivisions = (int) pow(2., (int) B.points.size2());
+ return &B;
+}
+
+std::map<int, JacobianBasis> JacobianBasis::fs;
- fs.insert(std::make_pair(tag, B));
- return &fs[tag];
+const JacobianBasis *JacobianBasis::find(int tag)
+{
+ std::map<int, JacobianBasis>::const_iterator it = fs.find(tag);
+ if (it != fs.end())
+ return &it->second;
+ JacobianBasis &B = fs[tag];
+ const polynomialBasis *F = polynomialBases::find(tag);
+ int jacobianOrder;
+ switch (F->parentType) {
+ case TYPE_PNT : jacobianOrder = 0; break;
+ case TYPE_LIN : jacobianOrder = F->order - 1; break;
+ case TYPE_TRI : jacobianOrder = 2 * (F->order - 1); break;
+ case TYPE_TET : jacobianOrder = 3 * (F->order - 1); break;
+ case TYPE_QUA : jacobianOrder = 2 * F->order - 1; break;
+ case TYPE_PRI : jacobianOrder = F->order * 3 - 1; break; // o=3*ord-2 on triangle base and =3*ord-1 for third dimension
+ default :
+ Msg::Error("Unknown function space %d: reverting to TET_4", tag);
+ jacobianOrder = 0;
+ }
+ B.bezier = bezierBasis::find(polynomialBasis::getTag(F->parentType, jacobianOrder, false));
+ generateGradShapes(B, B.bezier->points, F->monomials, F->coefficients);
+ return &B;
}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 645c512c75..bcd6c248cf 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -9,23 +9,25 @@
#include <map>
#include <vector>
#include "fullMatrix.h"
-
-class JacobianBasis
-{
- public:
- int numLagPts;
- int numDivisions;
- fullMatrix<double> exposants; //exposants of Bezier FF
- fullMatrix<double> points; //sampling point
- fullMatrix<double> matrixLag2Bez;
- fullMatrix<double> gradShapeMatX;
- fullMatrix<double> gradShapeMatY;
- fullMatrix<double> gradShapeMatZ;
- fullMatrix<double> divisor;
+class bezierBasis {
+ private :
+ static std::map<int, bezierBasis> _bbs;
+ public :
+ int numLagPts;
+ int numDivisions;
+ fullMatrix<double> exposants; //exposants of Bezier FF
+ fullMatrix<double> points; //sampling point
+ fullMatrix<double> matrixLag2Bez;
+ fullMatrix<double> divisor;
+ static const bezierBasis *find(int);
};
-class JacobianBases
-{
+class JacobianBasis {
+ public :
+ const bezierBasis *bezier;
+ fullMatrix<double> gradShapeMatX;
+ fullMatrix<double> gradShapeMatY;
+ fullMatrix<double> gradShapeMatZ;
private:
static std::map<int, JacobianBasis> fs;
public:
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index aef8ae5a0f..e2bb340ccd 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -38,57 +38,92 @@ static void printClosure(polynomialBasis::clCont &fullClosure, std::vector<int>
}
}
-static int getTriangleType (int order)
+int polynomialBasis::getTag(int parentTag, int order, bool serendip)
{
- switch(order) {
- case 0 : return MSH_TRI_1;
- case 1 : return MSH_TRI_3;
- case 2 : return MSH_TRI_6;
- case 3 : return MSH_TRI_10;
- case 4 : return MSH_TRI_15;
- case 5 : return MSH_TRI_21;
- case 6 : return MSH_TRI_28;
- case 7 : return MSH_TRI_36;
- case 8 : return MSH_TRI_45;
- case 9 : return MSH_TRI_55;
- case 10 : return MSH_TRI_66;
- default : Msg::Error("triangle order %i unknown", order); return 0;
- }
-}
-
-static int getQuadType (int order)
-{
- switch(order) {
- case 0 : return MSH_QUA_1;
- case 1 : return MSH_QUA_4;
- case 2 : return MSH_QUA_9;
- case 3 : return MSH_QUA_16;
- case 4 : return MSH_QUA_25;
- case 5 : return MSH_QUA_36;
- case 6 : return MSH_QUA_49;
- case 7 : return MSH_QUA_64;
- case 8 : return MSH_QUA_81;
- case 9 : return MSH_QUA_100;
- case 10 : return MSH_QUA_121;
- default : Msg::Error("quad order %i unknown", order); return 0;
- }
-}
-
-static int getLineType (int order)
-{
- switch(order) {
- case 0 : return MSH_LIN_1;
- case 1 : return MSH_LIN_2;
- case 2 : return MSH_LIN_3;
- case 3 : return MSH_LIN_4;
- case 4 : return MSH_LIN_5;
- case 5 : return MSH_LIN_6;
- case 6 : return MSH_LIN_7;
- case 7 : return MSH_LIN_8;
- case 8 : return MSH_LIN_9;
- case 9 : return MSH_LIN_10;
- case 10 : return MSH_LIN_11;
- default : Msg::Error("line order %i unknown", order); return 0;
+ switch (parentTag) {
+ case TYPE_PNT :
+ return MSH_PNT;
+ case TYPE_LIN :
+ switch(order) {
+ case 0 : return MSH_LIN_1;
+ case 1 : return MSH_LIN_2;
+ case 2 : return MSH_LIN_3;
+ case 3 : return MSH_LIN_4;
+ case 4 : return MSH_LIN_5;
+ case 5 : return MSH_LIN_6;
+ case 6 : return MSH_LIN_7;
+ case 7 : return MSH_LIN_8;
+ case 8 : return MSH_LIN_9;
+ case 9 : return MSH_LIN_10;
+ case 10: return MSH_LIN_11;
+ default : Msg::Error("line order %i unknown", order); return 0;
+ }
+ case TYPE_TRI :
+ switch(order) {
+ case 0 : return MSH_TRI_1;
+ case 1 : return MSH_TRI_3;
+ case 2 : return MSH_TRI_6;
+ case 3 : return serendip ? MSH_TRI_9 : MSH_TRI_10;
+ case 4 : return serendip ? MSH_TRI_12 : MSH_TRI_15;
+ case 5 : return serendip ? MSH_TRI_15I: MSH_TRI_21;
+ case 6 : return serendip ? MSH_TRI_18 : MSH_TRI_28;
+ case 7 : return serendip ? MSH_TRI_21I: MSH_TRI_36;
+ case 8 : return serendip ? MSH_TRI_24 : MSH_TRI_45;
+ case 9 : return serendip ? MSH_TRI_27 : MSH_TRI_55;
+ case 10: return serendip ? MSH_TRI_30 : MSH_TRI_66;
+ default : Msg::Error("triangle order %i unknown", order); return 0;
+ }
+ case TYPE_QUA :
+ switch(order) {
+ case 0 : return MSH_QUA_1;
+ case 1 : return MSH_QUA_4;
+ case 2 : return serendip ? MSH_QUA_8 : MSH_QUA_9;
+ case 3 : return serendip ? MSH_QUA_12 : MSH_QUA_16;
+ case 4 : return serendip ? MSH_QUA_16I: MSH_QUA_25;
+ case 5 : return serendip ? MSH_QUA_20 : MSH_QUA_36;
+ case 6 : return serendip ? MSH_QUA_24 : MSH_QUA_49;
+ case 7 : return serendip ? MSH_QUA_28 : MSH_QUA_64;
+ case 8 : return serendip ? MSH_QUA_32 : MSH_QUA_81;
+ case 9 : return serendip ? MSH_QUA_36I: MSH_QUA_100;
+ case 10: return serendip ? MSH_QUA_40 : MSH_QUA_121;
+ default : Msg::Error("quad order %i unknown", order); return 0;
+ }
+ case TYPE_TET :
+ switch(order) {
+ case 0 : return MSH_TET_1;
+ case 1 : return MSH_TET_4;
+ case 2 : return MSH_TET_10;
+ case 3 : return MSH_TET_20;
+ case 4 : return serendip ? MSH_TET_34 : MSH_TET_35;
+ case 5 : return serendip ? MSH_TET_52 : MSH_TET_56;
+ case 6 : return serendip ? MSH_TET_74 : MSH_TET_84;
+ case 7 : return serendip ? MSH_TET_100: MSH_TET_120;
+ case 8 : return serendip ? MSH_TET_130: MSH_TET_165;
+ case 9 : return serendip ? MSH_TET_164: MSH_TET_220;
+ case 10: return serendip ? MSH_TET_202: MSH_TET_286;
+ default : Msg::Error("terahedron order %i unknown", order); return 0;
+ }
+ case TYPE_HEX :
+ switch(order) {
+ case 1 : return MSH_HEX_8;
+ case 2 : return serendip ? MSH_HEX_20 : MSH_HEX_27;
+ case 3 : return serendip ? MSH_HEX_56 : MSH_HEX_64;
+ case 4 : return serendip ? MSH_HEX_98 : MSH_HEX_125;
+ case 5 : return serendip ? MSH_HEX_152: MSH_HEX_216;
+ case 6 : return serendip ? MSH_HEX_222: MSH_HEX_343;
+ case 7 : return serendip ? MSH_HEX_296: MSH_HEX_512;
+ case 8 : return serendip ? MSH_HEX_386: MSH_HEX_729;
+ case 9 : return serendip ? MSH_HEX_488: MSH_HEX_1000;
+ default : Msg::Error("hexahedron order %i unknown", order); return 0;
+ }
+ case TYPE_PRI :
+ switch(order) {
+ case 0 : return MSH_PRI_1;
+ case 1 : return MSH_PRI_6;
+ case 2 : return MSH_PRI_18;
+ default : Msg::Error("prism order %i unknown", order); return 0;
+ }
+ default : Msg::Error("unknown element type %i", parentTag); return 0;
}
}
@@ -1045,7 +1080,7 @@ static void getFaceClosureTet(int iFace, int iSign, int iRotate, polynomialBasis
closure.clear();
closure.resize((order + 1) * (order + 2) / 2);
- closure.type = getTriangleType(order);
+ closure.type = polynomialBasis::getTag(TYPE_TRI, order, false);
switch (order){
case 0:
@@ -1138,7 +1173,7 @@ static void generate2dEdgeClosureFull(polynomialBasis::clCont &closure, std::vec
if (serendip) break;
}
for (int r = 0; r < nNod*2 ; r++) {
- closure[r].type = getLineType(order);
+ closure[r].type = polynomialBasis::getTag(TYPE_LIN, order);
closureRef[r] = 0;
}
}
@@ -1282,7 +1317,7 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
// int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8},
// {8, 13, 11, 7}, {9, 11, 14, 10}};
int nVertex = isTriangle ? 3 : 4;
- closure.type = isTriangle ? getTriangleType(order) : getQuadType(order);
+ closure.type = polynomialBasis::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order);
for (int i = 0; i < nVertex; ++i){
int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
closure[i] = order1node[iFace][k];
@@ -1394,7 +1429,7 @@ static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, i
closure[j].push_back( nNod + (order-1)*j + i );
closure[nNod+j].push_back(nNod + (order-1)*(j+1) -i -1);
}
- closure[j].type = closure[nNod+j].type = getLineType(order);
+ closure[j].type = closure[nNod+j].type = polynomialBasis::getTag(TYPE_LIN, order);
}
}
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 54131bbc0d..f26edd968f 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -279,6 +279,7 @@ class polynomialBasis
}
}
const fullMatrix<double> &getGradientAtFaceIntegrationPoints(int integrationOrder, int closureId) const;
+ static int getTag(int parentTag, int order, bool serendip = false);
};
class polynomialBases
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 5312f26f91..50097b259e 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -517,22 +517,23 @@ int GMSH_AnalyseCurvedMeshPlugin::method_1_1(MElement *el, int depth)
Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
return 0;
}
+ const bezierBasis *bb = jfs->bezier;
- int numSamplingPt = jfs->points.size1();
- int dim = jfs->points.size2();
+ int numSamplingPt = bb->points.size1();
+ int dim = bb->points.size2();
fullVector<double> jacobian(numSamplingPt);
for (int i = 0; i < numSamplingPt; i++) {
double gsf[256][3];
switch (dim) {
case 1 :
- fs->df(jfs->points(i,0),0,0, gsf);
+ fs->df(bb->points(i,0),0,0, gsf);
break;
case 2 :
- fs->df(jfs->points(i,0),jfs->points(i,1),0, gsf);
+ fs->df(bb->points(i,0),bb->points(i,1),0, gsf);
break;
case 3 :
- fs->df(jfs->points(i,0),jfs->points(i,1),jfs->points(i,2), gsf);
+ fs->df(bb->points(i,0),bb->points(i,1),bb->points(i,2), gsf);
break;
default :
Msg::Error("Can't get the gradient for %dD elements.", dim);
@@ -548,7 +549,7 @@ int GMSH_AnalyseCurvedMeshPlugin::method_1_1(MElement *el, int depth)
fullVector<double> jacBez(jacobian.size());
- jfs->matrixLag2Bez.mult(jacobian, jacBez);
+ bb->matrixLag2Bez.mult(jacobian, jacBez);
bool allPtPositive = true;
for (int i = 0; i < jacBez.size(); i++) {
@@ -561,7 +562,7 @@ int GMSH_AnalyseCurvedMeshPlugin::method_1_1(MElement *el, int depth)
return 0;
}
else{
- int tag = division(jfs, jacBez, depth-1);
+ int tag = division(bb, jacBez, depth-1);
if (tag < 0)
return tag - 1;
if (tag > 0)
@@ -583,9 +584,10 @@ int GMSH_AnalyseCurvedMeshPlugin::method_1_2(MElement *el, int depth)
Msg::Error("Jacobian function space not implemented for type of element %d", el->getNum());
return 0;
}
+ const bezierBasis *bb = jfs->bezier;
- int numSamplingPt = jfs->points.size1();
- int dim = jfs->points.size2();
+ int numSamplingPt = bb->points.size1();
+ int dim = bb->points.size2();
fullVector<double> jacobian(numSamplingPt);
@@ -606,7 +608,7 @@ int GMSH_AnalyseCurvedMeshPlugin::method_1_2(MElement *el, int depth)
fullVector<double> jacBez(jacobian.size());
- jfs->matrixLag2Bez.mult(jacobian, jacBez);
+ bb->matrixLag2Bez.mult(jacobian, jacBez);
bool allPtPositive = true;
for (int i = 0; i < jacBez.size(); i++) {
@@ -619,7 +621,7 @@ int GMSH_AnalyseCurvedMeshPlugin::method_1_2(MElement *el, int depth)
return 0;
}
else{
- int tag = division(jfs, jacBez, depth-1);
+ int tag = division(bb, jacBez, depth-1);
if (tag < 0)
return tag - 1;
if (tag > 0)
@@ -642,9 +644,10 @@ void GMSH_AnalyseCurvedMeshPlugin::method_2_2
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
+ const bezierBasis *bb = jfs->bezier;
- int numSamplingPt = jfs->points.size1();
- int dim = jfs->points.size2();
+ int numSamplingPt = bb->points.size1();
+ int dim = bb->points.size2();
int numEl = tags.size();
fullMatrix<double> jacobian(numSamplingPt, numEl);
@@ -671,7 +674,7 @@ void GMSH_AnalyseCurvedMeshPlugin::method_2_2
double critere = .2; // à définir suivant le temps de chaque opération
if ((double) numBad / (double) numEl < critere) {// il vaut mieux calculer les points de controle de chaque élément
jacBez.resize(numSamplingPt, numEl);
- jfs->matrixLag2Bez.mult(jacobian, jacBez);
+ bb->matrixLag2Bez.mult(jacobian, jacBez);
correspondingEl.resize(numEl);
for (int i = 0; i < numEl; i++) correspondingEl[i] = i;
}
@@ -690,7 +693,7 @@ void GMSH_AnalyseCurvedMeshPlugin::method_2_2
}
}
jacBez.resize(numSamplingPt, numEl-numBad);
- jfs->matrixLag2Bez.mult(jacobian, jacBez);
+ bb->matrixLag2Bez.mult(jacobian, jacBez);
}
@@ -728,7 +731,7 @@ void GMSH_AnalyseCurvedMeshPlugin::method_2_2
}
int GMSH_AnalyseCurvedMeshPlugin::division
- (const JacobianBasis *jfs, const fullVector<double> &jacobian, int depth)
+ (const bezierBasis *jfs, const fullVector<double> &jacobian, int depth)
{
if (jfs->divisor.size2() != jacobian.size()) {
Msg::Error("Wrong sizes in division : [%d,%d] * [%d]",
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index f00a97a87e..63de47ec14 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -38,7 +38,7 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
//int checkJacobian(MElement *, int depth);
//int *checkJacobian2(MElement *const *, int numEl, int depth);
int *checkJacobian(MElement *const *, int numEl, int depth, int method);
- int division(const JacobianBasis *, const fullVector<double> &, int depth);
+ int division(const bezierBasis *, const fullVector<double> &, int depth);
};
#endif
--
GitLab