Skip to content
Snippets Groups Projects
Commit bdd62177 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

rid of dPsiDx + fullClosures

parent d4b144fa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/polynomialBasis.cpp
+ 289
− 17
View file @ bdd62177
......@@ -5,6 +5,7 @@
//
// Contributor(s):
// Koen Hillewaert
// Jonathan Lambrechts
//
#include "GmshDefines.h"
......@@ -12,6 +13,25 @@
#include "polynomialBasis.h"
#include "Gauss.h"
#include "stdlib.h"
static void printClosure(polynomialBasis::clCont &fullClosure, std::vector<int> &closureRef, polynomialBasis::clCont &closures) {
for (int i = 0; i <closures.size(); i++) {
printf("%3i (%2i): ", i, closureRef[i]);
if(closureRef[i]==-1){
printf("\n");
continue;
}
for (int j = 0; j < fullClosure[i].size(); j++) {
printf("%i ", fullClosure[i][j]);
}
printf (" -- ");
for (int j = 0; j < closures[closureRef[i]].size(); j++) {
printf("%i-%i ", fullClosure[i][closures[closureRef[i]][j]], closures[i][j]);
}
printf("\n");
}
}
static int getTriangleType (int order) {
switch(order) {
case 0 : return MSH_TRI_1;
......@@ -861,6 +881,33 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
return coefficient;
}
static void generate1dVertexClosure(polynomialBasis::clCont &closure)
{
closure.clear();
closure.resize(2);
closure[0].push_back(0);
closure[1].push_back(1);
closure[0].type = MSH_PNT;
closure[1].type = MSH_PNT;
}
static void generate1dVertexClosureFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order)
{
closure.clear();
closure.resize(2);
closure[0].push_back(0);
closure[0].push_back(1);
closure[1].push_back(1);
closure[1].push_back(0);
for (int i = 0; i < order - 1; i++) {
closure[0].push_back(2 + i);
closure[1].push_back(2 + order - 2 - i);
}
closureRef.resize(2);
closureRef[0] = 0;
closureRef[1] = 0;
}
static void getFaceClosureTet(int iFace, int iSign, int iRotate, polynomialBasis::closure &closure,
int order)
{
......@@ -927,10 +974,130 @@ static void getFaceClosureTet(int iFace, int iSign, int iRotate, polynomialBasis
break;
}
}
static void generate2dEdgeClosureFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order, int nNod, bool serendip)
{
closure.clear();
closure.resize(2*nNod);
closureRef.resize(2*nNod);
int shift = 0;
for (int corder = order; corder>=0; corder -=3) {
if (corder == 0) {
for (int r = 0; r < nNod ; r++){
closure[r].push_back(shift);
closure[r+nNod].push_back(shift);
}
break;
}
for (int r = 0; r < nNod ; r++){
for (int j = 0; j < nNod; j++) {
closure[r].push_back(shift + (r + j) % nNod);
closure[r + nNod].push_back(shift + (r - j + 1 + nNod) % nNod);
}
}
shift += nNod;
int n = nNod*(corder-1);
for (int r = 0; r < nNod ; r++){
for (int j = 0; j < n; j++){
closure[r].push_back(shift + (j + (corder - 1) * r) % n);
closure[r + nNod].push_back(shift + (n - j - 1 + (corder - 1) * (r + 1)) % n);
}
}
shift += n;
if (serendip) break;
}
for (int r = 0; r < nNod*2 ; r++) {
closure[r].type = getLineType(order);
closureRef[r] = 0;
}
}
static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
static void generateFaceClosureTetFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order, bool serendip)
{
closure.clear();
//input :
static const short int faces[4][3] = {{0,1,2}, {0,3,1}, {0,2,3}, {3,2,1}};
static const short int edges[6][2] = {{0,1}, {1,2}, {2,0}, {3,0}, {3,2}, {3,1}};
static const int faceOrientation[6] = {0,1,2,5,3,4};
closure.resize(24);
closureRef.resize(24);
for (int i = 0; i < 24; i ++)
closureRef[i] = 0;
if (order == 0) {
for (int i = 0; i < 24; i ++) {
closure[i].push_back(0);
}
return;
}
//Mapping for the p1 nodes
for (int iRotate = 0; iRotate < 3; iRotate ++) {
for (int iFace = 0; iFace < 4; iFace ++) {
for (int iNode = 0; iNode < 3; iNode ++) {
closure[iFace + iRotate * 8 + 0].push_back(faces[iFace][(6 + iNode - iRotate) % 3]);
closure[iFace + iRotate * 8 + 4].push_back(faces[iFace][(6 - iNode - iRotate) % 3]);
}
}
}
short int nodes2Edges [4][4];
for (int i = 0; i < 6; i++) {
nodes2Edges[edges[i][0]][edges[i][1]] = i * 2;
nodes2Edges[edges[i][1]][edges[i][0]] = i * 2 + 1;
}
int nodes2Faces[4][4][4];
for (int i = 0; i < 4; i++) {
for (int iRotate = 0; iRotate < 3; iRotate ++) {
short int n0 = faces[i][(3 - iRotate) % 3];
short int n1 = faces[i][(4 - iRotate) % 3];
short int n2 = faces[i][(5 - iRotate) % 3];
nodes2Faces[n0][n1][n2] = i * 6 + iRotate;
nodes2Faces[n0][n2][n1] = i * 6 + iRotate + 3;
}
}
polynomialBasis::clCont closureTriangles;
std::vector<int> closureTrianglesRef;
if (order >= 3)
generate2dEdgeClosureFull(closureTriangles, closureTrianglesRef, order - 3, 3, false);
for (int iClosure = 0; iClosure < closure.size(); iClosure++) {
std::vector<int> &cl = closure[iClosure];
//last node
cl.push_back(6 - cl[0] - cl[1] - cl[2]);
//edges
for (int iEdge = 0; iEdge < 6; iEdge ++) {
int n0 = cl[edges[iEdge][0]];
int n1 = cl[edges[iEdge][1]];
short int &oEdge = nodes2Edges[n0][n1];
for (int i = 0 ; i < order - 1; i++)
cl.push_back(4 + oEdge/2 * (order - 1) + ((oEdge % 2) ? order - 2 - i : i));
}
//faces
if (order >= 3) {
for (int iFace = 0; iFace < 4; iFace ++) {
int n0 = cl[faces[iFace][0]];
int n1 = cl[faces[iFace][1]];
int n2 = cl[faces[iFace][2]];
short int id = nodes2Faces[n0][n1][n2];
short int iTriClosure = faceOrientation [id % 6];
short int idFace = id/6;
int nOnFace = closureTriangles[iTriClosure].size();
for (int i = 0; i < nOnFace; i++) {
cl.push_back(4 + 6 * (order - 1) + idFace * nOnFace + closureTriangles[iTriClosure][i]);
}
}
}
}
if (order >= 4 && !serendip) {
polynomialBasis::clCont insideClosure;
std::vector<int> fakeClosureRef;
generateFaceClosureTetFull(insideClosure, fakeClosureRef, order - 4, false);
for (int i = 0; i < closure.size(); i ++) {
int shift = closure[i].size();
for (int j = 0; j < insideClosure[i].size(); j++)
closure[i].push_back(insideClosure[i][j] + shift);
}
}
}
static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
{
closure.clear();
for (int iRotate = 0; iRotate < 3; iRotate++){
for (int iSign = 1; iSign >= -1; iSign -= 2){
......@@ -942,6 +1109,7 @@ static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
}
}
}
static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
polynomialBasis::closure &closure, int order)
{
......@@ -977,6 +1145,48 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
closure[nNodes-1] = order2node[iFace][4]; // center
}
}
static void generateFaceClosurePrismFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order) {
closure.clear();
//input :
static const int triRotation[3][6] = {{0,1,2,3,4,5},{1,2,0,4,5,3},{2,0,1,5,3,4}};
static const int triSign[2][6] = {{0,1,2,3,4,5},{0,2,1,3,5,4}};
static const int triFace[5][6] = {{0,2,1,3,5,4},{3,4,5,0,1,2}};
static const int quadFace[3][6] = {{0,1,2,3,4,5},{5,2,4,3,0,1},{1,2,0,4,5,3}};
static const int quadRotation[4][6] = {{0,1,2,3,4,5},{3,0,5,4,1,2},{4,3,5,1,0,2},{1,4,5,0,3,2}};
static const int quadSign[2][6] = {{0,1,2,3,4,5},{4,1,5,3,0,2}};
//static const short int edges[6][2] = {{0,1}, {1,2}, {2,0}, {3,0}, {3,2}, {3,1}};
closure.resize(40);
closureRef.resize(40);
if (order == 0) {
for (int i = 0; i < 40; i ++) {
closure[i].push_back(0);
}
return;
}
for (int i = 0; i < 40; i++)
closure[i].resize(6);
//Mapping for the p1 nodes
//triangular faces
for (int iRotate = 0; iRotate < 3; iRotate ++)
for (int iFace = 0; iFace < 2; iFace ++)
for (int iSign = 0; iSign < 2; iSign ++)
for (int iNode = 0; iNode < 6; iNode ++) {
int i = iFace + iRotate * 10 + iSign * 5;
closureRef[i] = 0;
closure[i][triFace[0][iNode]] =
triFace[iFace][triRotation[iRotate][triSign[iSign][iNode]]];
}
//quad Horiz face
for (int iRotate = 0; iRotate < 4; iRotate ++)
for (int iFace = 0; iFace < 3; iFace ++)
for (int iSign = 0; iSign < 2; iSign ++)
for (int iNode = 0; iNode < 6; iNode ++) {
int i = (2 + iFace) + iRotate * 10 + iSign * 5;
closureRef[i] = (iFace + iRotate + iSign) % 2 ? 3 : 2;
closure[i][quadFace[closureRef[i]-2][quadRotation[0][iNode]]] =
quadFace[iFace][quadRotation[iRotate][quadSign[iSign][iNode]]];
}
}
static void generateFaceClosurePrism(polynomialBasis::clCont &closure, int order)
{
......@@ -993,8 +1203,7 @@ static void generateFaceClosurePrism(polynomialBasis::clCont &closure, int order
}
}
static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order,
int nNod = 3)
static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, int nNod = 3)
{
closure.clear();
closure.resize(2*nNod);
......@@ -1011,16 +1220,6 @@ static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order,
}
}
static void generate1dVertexClosure(polynomialBasis::clCont &closure)
{
closure.clear();
closure.resize(2);
closure[0].push_back(0);
closure[1].push_back(1);
closure[0].type = MSH_PNT;
closure[1].type = MSH_PNT;
}
static void generateClosureOrder0(polynomialBasis::clCont &closure, int nb)
{
closure.clear();
......@@ -1039,7 +1238,6 @@ const polynomialBasis *polynomialBases::find(int tag)
if (it != fs.end()) return &it->second;
polynomialBasis F;
F.numFaces = -1;
switch (tag){
case MSH_PNT:
F.numFaces = 1;
......@@ -1052,6 +1250,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.monomials = generate1DMonomials(0);
F.points = generate1DPoints(0);
generateClosureOrder0(F.closures,2);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 0);
F.parentType = TYPE_LIN;
break;
case MSH_LIN_2 :
......@@ -1059,6 +1258,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.monomials = generate1DMonomials(1);
F.points = generate1DPoints(1);
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 1);
F.parentType = TYPE_LIN;
break;
case MSH_LIN_3 :
......@@ -1066,6 +1266,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.monomials = generate1DMonomials(2);
F.points = generate1DPoints(2);
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 2);
F.parentType = TYPE_LIN;
break;
case MSH_LIN_4:
......@@ -1074,6 +1275,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(3);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 3);
break;
case MSH_LIN_5:
F.numFaces = 2;
......@@ -1081,6 +1283,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(4);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 4);
break;
case MSH_LIN_6:
F.numFaces = 2;
......@@ -1088,6 +1291,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(5);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 5);
break;
case MSH_LIN_7:
F.numFaces = 2;
......@@ -1095,6 +1299,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(6);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 6);
break;
case MSH_LIN_8:
F.numFaces = 2;
......@@ -1102,6 +1307,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(7);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 7);
break;
case MSH_LIN_9:
F.numFaces = 2;
......@@ -1109,6 +1315,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(8);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 8);
break;
case MSH_LIN_10:
F.numFaces = 2;
......@@ -1116,6 +1323,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(9);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 9);
break;
case MSH_LIN_11:
F.numFaces = 2;
......@@ -1123,13 +1331,16 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(10);
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
generate1dVertexClosureFull(F.fullClosures, F.closureRef, 10);
break;
case MSH_TRI_1 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(0);
F.points = gmshGeneratePointsTriangle(0, false);
F.parentType = TYPE_TRI;
generateClosureOrder0(F.closures,6);
generateClosureOrder0(F.closures, 6);
generateClosureOrder0(F.fullClosures, 6);
F.closureRef.resize(6, 0);
break;
case MSH_TRI_3 :
F.numFaces = 3;
......@@ -1137,6 +1348,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(1, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 1);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 1, 3, false);
break;
case MSH_TRI_6 :
F.numFaces = 3;
......@@ -1144,6 +1356,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(2, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 2);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 2, 3, false);
break;
case MSH_TRI_9 :
F.numFaces = 3;
......@@ -1151,6 +1364,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(3, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 3);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 3, true);
break;
case MSH_TRI_10 :
F.numFaces = 3;
......@@ -1158,6 +1372,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(3, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 3);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 3, false);
break;
case MSH_TRI_12 :
F.numFaces = 3;
......@@ -1165,6 +1380,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(4, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 3, true);
break;
case MSH_TRI_15 :
F.numFaces = 3;
......@@ -1172,6 +1388,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(4, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 3, false);
break;
case MSH_TRI_15I :
F.numFaces = 3;
......@@ -1179,6 +1396,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(5, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 5);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 3, true);
break;
case MSH_TRI_21 :
F.numFaces = 3;
......@@ -1186,6 +1404,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(5, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 5);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 3, false);
break;
case MSH_TRI_28 :
F.numFaces = 3;
......@@ -1193,6 +1412,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(6, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 6);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 3, false);
break;
case MSH_TRI_36 :
F.numFaces = 3;
......@@ -1200,6 +1420,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(7, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 7);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 3, false);
break;
case MSH_TRI_45 :
F.numFaces = 3;
......@@ -1207,6 +1428,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(8, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 8);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 3, false);
break;
case MSH_TRI_55 :
F.numFaces = 3;
......@@ -1214,6 +1436,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(9, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 9);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 3, false);
break;
case MSH_TRI_66 :
F.numFaces = 3;
......@@ -1221,6 +1444,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(10, false);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 10);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 3, false);
break;
case MSH_TRI_18 :
F.numFaces = 3;
......@@ -1228,6 +1452,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(6, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 6);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 3, true);
break;
case MSH_TRI_21I :
F.numFaces = 3;
......@@ -1235,6 +1460,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(7, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 7);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 3, true);
break;
case MSH_TRI_24 :
F.numFaces = 3;
......@@ -1242,6 +1468,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(8, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 8);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 3, true);
break;
case MSH_TRI_27 :
F.numFaces = 3;
......@@ -1249,6 +1476,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(9, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 9);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 3, true);
break;
case MSH_TRI_30 :
F.numFaces = 3;
......@@ -1256,6 +1484,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTriangle(10, true);
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 10);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 3, true);
break;
case MSH_TET_1 :
F.numFaces = 4;
......@@ -1263,6 +1492,8 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(0, false);
F.parentType = TYPE_TET;
generateClosureOrder0(F.closures,24);
generateClosureOrder0(F.fullClosures, 24);
F.closureRef.resize(24, 0);
break;
case MSH_TET_4 :
F.numFaces = 4;
......@@ -1270,6 +1501,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(1, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 1);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 1, false);
break;
case MSH_TET_10 :
F.numFaces = 4;
......@@ -1277,6 +1509,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(2, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 2);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 2, false);
break;
case MSH_TET_20 :
F.numFaces = 4;
......@@ -1284,6 +1517,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(3, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 3);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 3, false);
break;
case MSH_TET_35 :
F.numFaces = 4;
......@@ -1291,6 +1525,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(4, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 4);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 4, false);
break;
case MSH_TET_34 :
F.numFaces = 4;
......@@ -1298,6 +1533,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(4, true);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 4);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 4, true);
break;
case MSH_TET_52 :
F.numFaces = 4;
......@@ -1305,6 +1541,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(5, true);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 5);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 5, true);
break;
case MSH_TET_56 :
F.numFaces = 4;
......@@ -1312,6 +1549,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(5, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 5);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 5, false);
break;
case MSH_TET_74 :
F.numFaces = 4;
......@@ -1319,6 +1557,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(6, true);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 6);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 6, true);
break;
case MSH_TET_84 :
F.numFaces = 4;
......@@ -1326,6 +1565,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(6, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 6);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 6, false);
break;
case MSH_TET_100 :
F.numFaces = 4;
......@@ -1333,6 +1573,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(7, true);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 7);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 7, true);
break;
case MSH_TET_120 :
F.numFaces = 4;
......@@ -1340,6 +1581,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(7, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 7);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 7, false);
break;
case MSH_TET_130 :
F.numFaces = 4;
......@@ -1347,6 +1589,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(8, true);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 8);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 8, true);
break;
case MSH_TET_164 :
F.numFaces = 4;
......@@ -1354,6 +1597,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(9, true);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 9);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 9, true);
break;
case MSH_TET_165 :
F.numFaces = 4;
......@@ -1361,6 +1605,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(8, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 8);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 8, false);
break;
case MSH_TET_202 :
F.numFaces = 4;
......@@ -1368,6 +1613,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(10, true);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 10);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 10, true);
break;
case MSH_TET_220 :
F.numFaces = 4;
......@@ -1375,6 +1621,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(9, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 9);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 9, false);
break;
case MSH_TET_286 :
F.numFaces = 4;
......@@ -1382,6 +1629,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(10, false);
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 10);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, 10, false);
break;
case MSH_QUA_1 :
F.numFaces = 4;
......@@ -1389,6 +1637,8 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(0,false);
F.parentType = TYPE_QUA;
generateClosureOrder0(F.closures,8);
generateClosureOrder0(F.fullClosures,8);
F.closureRef.resize(8, 0);
break;
case MSH_QUA_4 :
F.numFaces = 4;
......@@ -1396,6 +1646,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(1,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 1, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 1, 4, false);
break;
case MSH_QUA_9 :
F.numFaces = 4;
......@@ -1403,6 +1654,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(2,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 2, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 2, 4, false);
break;
case MSH_QUA_16 :
F.numFaces = 4;
......@@ -1410,6 +1662,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(3,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 3, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 4, false);
break;
case MSH_QUA_25 :
F.numFaces = 4;
......@@ -1417,6 +1670,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(4,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 4, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 4, false);
break;
case MSH_QUA_36 :
F.numFaces = 4;
......@@ -1424,6 +1678,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(5,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 5, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 4, false);
break;
case MSH_QUA_49 :
F.numFaces = 4;
......@@ -1431,6 +1686,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(6,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 6, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 4, false);
break;
case MSH_QUA_64 :
F.numFaces = 4;
......@@ -1438,6 +1694,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(7,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 7, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 4, false);
break;
case MSH_QUA_81 :
F.numFaces = 4;
......@@ -1445,6 +1702,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(8,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 8, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 4, false);
break;
case MSH_QUA_100 :
F.numFaces = 4;
......@@ -1452,6 +1710,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(9,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 9, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 4, false);
break;
case MSH_QUA_121 :
F.numFaces = 4;
......@@ -1459,6 +1718,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(10,false);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 10, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 4, false);
break;
case MSH_QUA_8 :
F.numFaces = 4;
......@@ -1466,6 +1726,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(2,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 2, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 2, 4, true);
break;
case MSH_QUA_12 :
F.numFaces = 4;
......@@ -1473,6 +1734,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(3,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 3, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 4, true);
break;
case MSH_QUA_16I :
F.numFaces = 4;
......@@ -1480,6 +1742,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(4,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 4, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 4, true);
break;
case MSH_QUA_20 :
F.numFaces = 4;
......@@ -1487,6 +1750,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(5,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 5, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 4, true);
break;
case MSH_QUA_24 :
F.numFaces = 4;
......@@ -1494,6 +1758,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(6,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 6, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 4, true);
break;
case MSH_QUA_28 :
F.numFaces = 4;
......@@ -1501,6 +1766,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(7,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 7, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 4, true);
break;
case MSH_QUA_32 :
F.numFaces = 4;
......@@ -1508,6 +1774,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(8,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 8, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 4, true);
break;
case MSH_QUA_36I :
F.numFaces = 4;
......@@ -1515,6 +1782,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(9,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 9, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 4, true);
break;
case MSH_QUA_40 :
F.numFaces = 4;
......@@ -1522,6 +1790,7 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsQuad(10,true);
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 10, 4);
generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 4, true);
break;
case MSH_PRI_1 :
F.numFaces = 5;
......@@ -1529,6 +1798,8 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsPrism(0, false);
F.parentType = TYPE_PRI;
generateClosureOrder0(F.closures,48);
generateClosureOrder0(F.fullClosures,48);
F.closureRef.resize(48, 0);
break;
case MSH_PRI_6 :
F.numFaces = 5;
......@@ -1536,14 +1807,15 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsPrism(1, false);
F.parentType = TYPE_PRI;
generateFaceClosurePrism(F.closures, 1);
generateFaceClosurePrismFull(F.fullClosures, F.closureRef, 1);
break;
case MSH_PRI_18 :
/*case MSH_PRI_18 :
F.numFaces = 5;
F.monomials = generatePascalPrism(2);
F.points = gmshGeneratePointsPrism(2, false);
F.parentType = TYPE_PRI;
generateFaceClosurePrism(F.closures, 2);
break;
break;*/
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.numFaces = 4;
......
Numeric/polynomialBasis.h
+ 7
− 1
View file @ bdd62177
......@@ -80,7 +80,13 @@ class polynomialBasis
int type;
};
typedef std::vector<closure> clCont;
clCont closures;
// 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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment