From f552f5892f68f7b9dc1f984e0e96c53cf29b9e30 Mon Sep 17 00:00:00 2001 From: Sebastien Blaise <sebastien.blaise@uclouvain.be> Date: Mon, 29 Oct 2012 15:28:59 +0000 Subject: [PATCH] gmsh: Added functions to nodalBasis for dg to compile. (I Think) those functions should be there and not in polynomialBasis, as they are general to any base. dg: switch from polynomialBasis to nodalBasis --- Numeric/Numeric.cpp | 2 - Numeric/Numeric.h | 1 + Numeric/nodalBasis.h | 117 +++++++++++++++++++++++++++++++++++- Numeric/polynomialBasis.cpp | 100 ++---------------------------- Numeric/polynomialBasis.h | 17 +----- 5 files changed, 122 insertions(+), 115 deletions(-) diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp index b6ec2dde2a..ae4b8b7c7b 100644 --- a/Numeric/Numeric.cpp +++ b/Numeric/Numeric.cpp @@ -7,8 +7,6 @@ #include "GmshMessage.h" #include "Numeric.h" -#define SQU(a) ((a)*(a)) - double myatan2(double a, double b) { if(a == 0.0 && b == 0) diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h index 6d9b879fc1..85ba4cdabc 100644 --- a/Numeric/Numeric.h +++ b/Numeric/Numeric.h @@ -14,6 +14,7 @@ #define myhypot(a,b) (sqrt((a)*(a)+(b)*(b))) #define sign(x) (((x)>=0)?1:-1) +#define SQU(a) ((a)*(a)) struct mean_plane { diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h index 94791d6bc7..3570b76959 100644 --- a/Numeric/nodalBasis.h +++ b/Numeric/nodalBasis.h @@ -7,6 +7,7 @@ #define BASIS_H #include "fullMatrix.h" +#include "GmshDefines.h" class nodalBasis { @@ -23,10 +24,122 @@ class nodalBasis { // Basis functions gradients evaluation inline virtual void df(double u, double v, double w, double grads[][3]) const {}; - inline void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const; + inline void df(fullMatrix<double> &coord, fullMatrix<double> &dfm) const {}; inline void ddf(double u, double v, double w, double grads[][3][3]) const {}; inline void dddf(double u, double v, double w, double grads[][3][3][3]) const {}; + inline virtual int getNumShapeFunctions() const {Msg::Fatal("Not implemented"); return -1;} + + class closure : public std::vector<int> { + public: + int type; + }; + typedef std::vector<closure> clCont; + // 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 + // considered face becomes the first one in the right orientation; + // For element, like prisms that have different kind of faces, + // fullCLosure[i] rotates the element so that the considered face + // becomes the closureRef[i]-th face (the first tringle or the first + // quad face) + clCont closures, fullClosures; + std::vector<int> closureRef; + inline virtual int getClosureType(int id) const {Msg::Fatal("Not implemented"); return -1;} + inline virtual const std::vector<int> &getClosure(int id) const {Msg::Fatal("Not implemented"); std::vector<int> *ret=NULL; return *ret;} + inline virtual const std::vector<int> &getFullClosure(int id) const {Msg::Fatal("Not implemented"); std::vector<int> *ret=NULL; return *ret;} + inline virtual int getClosureId(int iFace, int iSign=1, int iRot=0) const {Msg::Fatal("Not implemented"); return -1;} + inline void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const {Msg::Fatal("Not implemented");iFace=-1; iSign=-1; iRot=-1;} + + static int getTag(int parentTag, int order, bool serendip = false) + { + 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; + } + } + }; -#endif \ No newline at end of file +#endif diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp index ea67a22634..1965fb9542 100644 --- a/Numeric/polynomialBasis.cpp +++ b/Numeric/polynomialBasis.cpp @@ -43,96 +43,6 @@ static void printClosure(polynomialBasis::clCont &fullClosure, */ - -int polynomialBasis::getTag(int parentTag, int order, bool serendip) -{ - 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; - } -} - fullMatrix<double> generate1DMonomials(int order) { fullMatrix<double> monomials(order + 1, 1); @@ -1024,7 +934,7 @@ static void getFaceClosureTet(int iFace, int iSign, int iRotate, { closure.clear(); closure.resize((order + 1) * (order + 2) / 2); - closure.type = polynomialBasis::getTag(TYPE_TRI, order, false); + closure.type = nodalBasis::getTag(TYPE_TRI, order, false); switch (order){ case 0: @@ -1119,7 +1029,7 @@ static void generate2dEdgeClosureFull(polynomialBasis::clCont &closure, if (serendip) break; } for (int r = 0; r < nNod*2 ; r++) { - closure[r].type = polynomialBasis::getTag(TYPE_LIN, order); + closure[r].type = nodalBasis::getTag(TYPE_LIN, order); closureRef[r] = 0; } } @@ -1322,7 +1232,7 @@ static void generateFaceClosureHex(polynomialBasis::clCont &closure, int order, { closure.clear(); const polynomialBasis &fsFace = *polynomialBases::find - (polynomialBasis::getTag(TYPE_QUA, order, serendip)); + (nodalBasis::getTag(TYPE_QUA, order, serendip)); for (int iRotate = 0; iRotate < 4; iRotate++){ for (int iSign = 1; iSign >= -1; iSign -= 2){ for (int iFace = 0; iFace < 6; iFace++) { @@ -1407,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 = polynomialBasis::getTag(isTriangle ? TYPE_TRI : TYPE_QUA, order); + closure.type = nodalBasis::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]; @@ -1523,7 +1433,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 = polynomialBasis::getTag(TYPE_LIN, order); + closure[j].type = closure[nNod+j].type = nodalBasis::getTag(TYPE_LIN, order); } } diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h index 72e5a01c2d..ad41453e84 100644 --- a/Numeric/polynomialBasis.h +++ b/Numeric/polynomialBasis.h @@ -80,21 +80,7 @@ class polynomialBasis : public nodalBasis // basis, we use the type of the corresponding gmsh element as type //int type, parentType, order, dimension; //bool serendip; - class closure : public std::vector<int> { - public: - int type; - }; - typedef std::vector<closure> clCont; - // 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 - // considered face becomes the first one in the right orientation; - // For element, like prisms that have different kind of faces, - // fullCLosure[i] rotates the element so that the considered face - // becomes the closureRef[i]-th face (the first tringle or the first - // quad face) - clCont closures, fullClosures; - std::vector<int> closureRef; + //fullMatrix<double> points; fullMatrix<double> monomials; fullMatrix<double> coefficients; @@ -456,7 +442,6 @@ class polynomialBasis : public nodalBasis break; } } - static int getTag(int parentTag, int order, bool serendip = false); }; class polynomialBases -- GitLab