From b604d39b69adf8f60529b9eb5bf84fd352872d95 Mon Sep 17 00:00:00 2001
From: Karim Slaoui <karim.slaoui@uclouvain.be>
Date: Tue, 2 Nov 2010 11:03:38 +0000
Subject: [PATCH] Added P0 elements
---
Common/GmshDefines.h | 6 +++++
Geo/MLine.cpp | 1 +
Geo/MPrism.cpp | 2 ++
Geo/MQuadrangle.cpp | 2 ++
Geo/MTetrahedron.cpp | 2 ++
Geo/MTriangle.cpp | 4 ++++
Numeric/polynomialBasis.cpp | 48 +++++++++++++++++++++++++++++++++++++
7 files changed, 65 insertions(+)
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index adc4bae068..b954314d4d 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -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
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 4771a12613..66e84489c2 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -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);
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 22586f844d..ea385278fb 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -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);
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 13eb65e2cc..0563ae6d29 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -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);
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index aea4336b8d..83785d50e3 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -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);
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 721fa4de22..0c4f6322ba 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -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);
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 79d31aca83..da15f93cde 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -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);
--
GitLab