Skip to content
Snippets Groups Projects
Commit b604d39b authored by Karim Slaoui's avatar Karim Slaoui
Browse files

Added P0 elements

parent 2341211b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Common/GmshDefines.h
+ 6
− 0
View file @ b604d39b
......@@ -142,6 +142,12 @@
#define MSH_TET_164 82
#define MSH_TET_202 83
#define MSH_NUM_TYPE 84
#define MSH_LIN_1 85
#define MSH_TRI_1 86
#define MSH_QUA_1 87
#define MSH_TET_1 88
#define MSH_HEX_1 89
#define MSH_PRI_1 90
// Geometric entities
#define ENT_NONE 0
......
Geo/MLine.cpp
+ 1
− 0
View file @ b604d39b
......@@ -14,6 +14,7 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
case 0: return polynomialBases::find(MSH_LIN_1);
case 1: return polynomialBases::find(MSH_LIN_2);
case 2: return polynomialBases::find(MSH_LIN_3);
case 3: return polynomialBases::find(MSH_LIN_4);
......
Geo/MPrism.cpp
+ 2
− 0
View file @ b604d39b
......@@ -38,6 +38,7 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
case 0: return polynomialBases::find(MSH_PRI_1);
case 1: return polynomialBases::find(MSH_PRI_6);
case 2: return polynomialBases::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
......@@ -45,6 +46,7 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
}
else {
switch (order) {
case 0: return polynomialBases::find(MSH_PRI_1);
case 1: return polynomialBases::find(MSH_PRI_6);
case 2: return polynomialBases::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
......
Geo/MQuadrangle.cpp
+ 2
− 0
View file @ b604d39b
......@@ -24,6 +24,7 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
case 0: return polynomialBases::find(MSH_QUA_1);
case 1: return polynomialBases::find(MSH_QUA_4);
case 2: return polynomialBases::find(MSH_QUA_8);
case 3: return polynomialBases::find(MSH_QUA_12);
......@@ -37,6 +38,7 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
}
}
switch (order) {
case 0: return polynomialBases::find(MSH_QUA_1);
case 1: return polynomialBases::find(MSH_QUA_4);
case 2: return polynomialBases::find(MSH_QUA_9);
case 3: return polynomialBases::find(MSH_QUA_16);
......
Geo/MTetrahedron.cpp
+ 2
− 0
View file @ b604d39b
......@@ -125,6 +125,7 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
case 0: return polynomialBases::find(MSH_TET_1);
case 1: return polynomialBases::find(MSH_TET_4);
case 2: return polynomialBases::find(MSH_TET_10);
case 3: return polynomialBases::find(MSH_TET_20);
......@@ -140,6 +141,7 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
}
else {
switch (order) {
case 0: return polynomialBases::find(MSH_TET_1);
case 1: return polynomialBases::find(MSH_TET_4);
case 2: return polynomialBases::find(MSH_TET_10);
case 3: return polynomialBases::find(MSH_TET_20);
......
Geo/MTriangle.cpp
+ 4
− 0
View file @ b604d39b
......@@ -84,6 +84,7 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
case 0: return polynomialBases::find(MSH_TRI_1);
case 1: return polynomialBases::find(MSH_TRI_3);
case 2: return polynomialBases::find(MSH_TRI_6);
case 3: return polynomialBases::find(MSH_TRI_9);
......@@ -99,6 +100,7 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
}
else {
switch (order) {
case 0: return polynomialBases::find(MSH_TRI_1);
case 1: return polynomialBases::find(MSH_TRI_3);
case 2: return polynomialBases::find(MSH_TRI_6);
case 3: return polynomialBases::find(MSH_TRI_10);
......@@ -117,12 +119,14 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
const JacobianBasis* MTriangle::getJacobianFuncSpace(int o) const
{
int order = (o == -1) ? getPolynomialOrder() : o;
int nf = getNumFaceVertices();
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);
......
Numeric/polynomialBasis.cpp
+ 48
− 0
View file @ b604d39b
......@@ -14,6 +14,7 @@
static int getTriangleType (int order) {
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;
......@@ -29,6 +30,7 @@ static int getTriangleType (int order) {
}
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;
......@@ -44,6 +46,7 @@ static int getQuadType (int order) {
}
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;
......@@ -1018,6 +1021,16 @@ static void generate1dVertexClosure(polynomialBasis::clCont &closure)
closure[1].type = MSH_PNT;
}
static void generateClosureOrder0(polynomialBasis::clCont &closure, int nb)
{
closure.clear();
closure.resize(nb);
for (int i=0; i < nb; i++) {
closure[i].push_back(0);
closure[i].type = MSH_PNT;
}
}
std::map<int, polynomialBasis> polynomialBases::fs;
const polynomialBasis *polynomialBases::find(int tag)
......@@ -1034,6 +1047,13 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = generate1DPoints(0);
F.parentType = TYPE_PNT;
break;
case MSH_LIN_1 :
F.numFaces = 2;
F.monomials = generate1DMonomials(0);
F.points = generate1DPoints(0);
generateClosureOrder0(F.closures,2);
F.parentType = TYPE_LIN;
break;
case MSH_LIN_2 :
F.numFaces = 2;
F.monomials = generate1DMonomials(1);
......@@ -1104,6 +1124,13 @@ const polynomialBasis *polynomialBases::find(int tag)
F.parentType = TYPE_LIN;
generate1dVertexClosure(F.closures);
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);
break;
case MSH_TRI_3 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(1);
......@@ -1230,6 +1257,13 @@ const polynomialBasis *polynomialBases::find(int tag)
F.parentType = TYPE_TRI;
generate2dEdgeClosure(F.closures, 10);
break;
case MSH_TET_1 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(0);
F.points = gmshGeneratePointsTetrahedron(0, false);
F.parentType = TYPE_TET;
generateClosureOrder0(F.closures,24);
break;
case MSH_TET_4 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(1);
......@@ -1349,6 +1383,13 @@ const polynomialBasis *polynomialBases::find(int tag)
F.parentType = TYPE_TET;
generateFaceClosureTet(F.closures, 10);
break;
case MSH_QUA_1 :
F.numFaces = 4;
F.monomials = generatePascalQuad(0);
F.points = gmshGeneratePointsQuad(0,false);
F.parentType = TYPE_QUA;
generateClosureOrder0(F.closures,8);
break;
case MSH_QUA_4 :
F.numFaces = 4;
F.monomials = generatePascalQuad(1);
......@@ -1482,6 +1523,13 @@ const polynomialBasis *polynomialBases::find(int tag)
F.parentType = TYPE_QUA;
generate2dEdgeClosure(F.closures, 10, 4);
break;
case MSH_PRI_1 :
F.numFaces = 5;
F.monomials = generatePascalPrism(0);
F.points = gmshGeneratePointsPrism(1, false);
F.parentType = TYPE_PRI;
generateClosureOrder0(F.closures,48);
break;
case MSH_PRI_6 :
F.numFaces = 5;
F.monomials = generatePascalPrism(1);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment