From 232b3f7a7e325507b68657baa32a9489b2bd1d58 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Fri, 24 Sep 2010 15:16:41 +0000
Subject: [PATCH] missing bits for tets up to order 9
---
Common/GmshDefines.h | 8 ++-
Geo/MElement.cpp | 10 ++++
Geo/MTetrahedron.cpp | 106 ++++++++++++++++++++++++++++++++++++
Geo/MTetrahedron.h | 45 ++++++---------
Mesh/HighOrder.cpp | 98 +++++++++++++++++++++++----------
Numeric/polynomialBasis.cpp | 30 ++++++++++
6 files changed, 238 insertions(+), 59 deletions(-)
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 64f8d3233c..adc4bae068 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -136,8 +136,12 @@
#define MSH_HEX_64 76
#define MSH_HEX_125 77
#define MSH_HEX_196 78
-
-#define MSH_NUM_TYPE 75
+#define MSH_TET_74 79
+#define MSH_TET_100 80
+#define MSH_TET_130 81
+#define MSH_TET_164 82
+#define MSH_TET_202 83
+#define MSH_NUM_TYPE 84
// Geometric entities
#define ENT_NONE 0
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 48ee94d70d..9e23906295 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -866,6 +866,11 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_TET_35 : if(name) *name = "Tetrahedron 35"; return 4 + 18 + 12 + 1;
case MSH_TET_52 : if(name) *name = "Tetrahedron 52"; return 4 + 24 + 24 + 0;
case MSH_TET_56 : if(name) *name = "Tetrahedron 56"; return 4 + 24 + 24 + 4;
+ case MSH_TET_84 : if(name) *name = "Tetrahedron 84"; return (7*8*9)/6;
+ case MSH_TET_120 : if(name) *name = "Tetrahedron 120"; return (8*9*10)/6;
+ case MSH_TET_165 : if(name) *name = "Tetrahedron 165"; return (9*10*11)/6;
+ case MSH_TET_220 : if(name) *name = "Tetrahedron 220"; return (10*11*12)/6;
+ case MSH_TET_286 : if(name) *name = "Tetrahedron 286"; return (11*12*13)/6;
case MSH_HEX_8 : if(name) *name = "Hexahedron 8"; return 8;
case MSH_HEX_20 : if(name) *name = "Hexahedron 20"; return 8 + 12;
case MSH_HEX_27 : if(name) *name = "Hexahedron 27"; return 8 + 12 + 6 + 1;
@@ -1024,6 +1029,11 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_TET_35: return new MTetrahedronN(v, 4, num, part);
case MSH_TET_52: return new MTetrahedronN(v, 5, num, part);
case MSH_TET_56: return new MTetrahedronN(v, 5, num, part);
+ case MSH_TET_84: return new MTetrahedronN(v, 6, num, part);
+ case MSH_TET_120: return new MTetrahedronN(v, 7, num, part);
+ case MSH_TET_165: return new MTetrahedronN(v, 8, num, part);
+ case MSH_TET_220: return new MTetrahedronN(v, 9, num, part);
+ case MSH_TET_286: return new MTetrahedronN(v, 10, num, part);
case MSH_POLYH_: return new MPolyhedron(v, num, part, owner, parent);
default: return 0;
}
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 145fbd6fa2..cc44a7c6fb 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -130,6 +130,11 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
case 3: return polynomialBases::find(MSH_TET_20);
case 4: return polynomialBases::find(MSH_TET_34);
case 5: return polynomialBases::find(MSH_TET_52);
+ case 6: return polynomialBases::find(MSH_TET_74);
+ case 7: return polynomialBases::find(MSH_TET_100);
+ case 8: return polynomialBases::find(MSH_TET_130);
+ case 9: return polynomialBases::find(MSH_TET_164);
+ case 10: return polynomialBases::find(MSH_TET_202);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
@@ -356,3 +361,104 @@ void MTetrahedron::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &
}
Msg::Error("Could not get face information for tetrahedron %d", getNum());
}
+
+static std::vector<std::vector<int> > tetReverseIndices(20);
+
+const std::vector<int> &MTetrahedronN::_getReverseIndices (int order) {
+ if(order>=tetReverseIndices.size())
+ tetReverseIndices.resize(order+1);
+ std::vector<int> &r = tetReverseIndices[order];
+ if (r.size() != 0) return r;
+ //
+ // not the funniest code ever ... (guaranteed correct only up to order 5)
+ //
+ int nb = (order+1)*(order+2)*(order+3)/6;
+ r.resize(nb);
+ int p=0;
+ for (int layerOrder = order; layerOrder>=0; layerOrder-=4) {
+ //principal vertices
+ r[p+0] = p+0;
+ if (layerOrder ==0) break;
+ r[p+1] = p+2;
+ r[p+2] = p+1;
+ r[p+3] = p+3;
+ p+=4;
+ for (int i = 0; i<layerOrder-1; i++) {
+ //E2 reversed switches with E0
+ r[p+i] = p+3*(layerOrder-1)-(i+1);
+ r[p+3*(layerOrder-1)-(i+1)] = p+i;
+ //E1 is reversed
+ r[p+(layerOrder-1)+i] = p+2*(layerOrder-1)-(i+1);
+ //E3 is preserved
+ r[p+3*(layerOrder-1)+i] = p+3*(layerOrder-1)+i;
+ //E4 switches with E5
+ r[p+4*(layerOrder-1)+i] = p+5*(layerOrder-1)+i;
+ r[p+5*(layerOrder-1)+i] = p+4*(layerOrder-1)+i;
+ }
+ p+=6*(layerOrder-1);
+ //F0(=012) switches its nodes 1 and 2
+ for (int of = layerOrder-3; of >= 0; of -= 3) {
+ r[p] = p;
+ if (of == 0) {
+ p+=1;
+ break;
+ }
+ r[p+1] = p+2;
+ r[p+2] = p+1;
+ for (int i = 0; i < of-1; i++) {
+ //switch edges 0 and 2
+ r[p+3+i] = p+3+3*(of-1)-(i+1);
+ r[p+3+3*(of-1)-(i+1)] = p+3+i;
+ //reverse edge 1
+ r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1);
+ }
+ p += 3*of;
+ }
+ //F1 (=013) reversed switches with F2 (=032)
+ int nf = (layerOrder-2)*(layerOrder-1)/2;
+ for (int of = layerOrder-3; of >= 0; of -= 3) {
+ r[p] = p+nf;
+ r[p+nf] = p;
+ if (of == 0) {
+ p += 1;
+ break;
+ }
+ r[p+1] = p+nf+2;
+ r[p+nf+2] = p+1;
+ r[p+2] = p+nf+1;
+ r[p+nf+1] = p+2;
+ for (int i = 0; i < of-1; i++) {
+ //switch edges 0 and 2
+ r[p+3+i] = p+3+3*(of-1)-(i+1)+nf;
+ r[p+3+3*(of-1)-(i+1)] = p+3+i+nf;
+ r[p+3+i+nf] = p+3+3*(of-1)-(i+1);
+ r[p+3+3*(of-1)-(i+1)+nf] = p+3+i;
+ //reverse edge 1
+ r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1)+nf;
+ r[p+3+(of-1)+i+nf] = p+3+2*(of-1)-(i+1);
+ }
+ p += 3*of;
+ }
+ p+=nf;
+
+ //F3(=312) switches its nodes 1 and 2
+ for (int of = layerOrder-3; of >= 0; of -= 3) {
+ r[p] = p;
+ if (of == 0) {
+ p += 1;
+ break;
+ }
+ r[p+1] = p+2;
+ r[p+2] = p+1;
+ for (int i = 0; i < of-1; i++) {
+ //switch edges 0 and 2
+ r[p+3+i] = p+3+3*(of-1)-(i+1);
+ r[p+3+3*(of-1)-(i+1)] = p+3+i;
+ //reverse edge 1
+ r[p+3+(of-1)+i] = p+3+2*(of-1)-(i+1);
+ }
+ p += 3*of;
+ }
+ }
+ return r;
+}
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index dde4c80ff1..6ae19d0f94 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -347,6 +347,7 @@ static int reverseTet34[34] = {0,2,1,3, // principal vertices
31,33,32};// F3 is uvw plane, reverts normal
class MTetrahedronN : public MTetrahedron {
+ const static std::vector<int> &_getReverseIndices(int order);
protected:
std::vector<MVertex *> _vs;
const char _order;
@@ -396,6 +397,11 @@ class MTetrahedronN : public MTetrahedron {
switch(getTypeForMSH()){
case MSH_TET_35 : return 1;
case MSH_TET_56 : return 4;
+ case MSH_TET_84 : return 10;
+ case MSH_TET_120 : return 20;
+ case MSH_TET_165 : return 35;
+ case MSH_TET_220 : return 56;
+ case MSH_TET_286 : return 84;
default : return 0;
}
}
@@ -407,41 +413,22 @@ class MTetrahedronN : public MTetrahedron {
if(_order == 4 && _vs.size() + 4 == 35) return MSH_TET_35;
if(_order == 5 && _vs.size() + 4 == 56) return MSH_TET_56;
if(_order == 5 && _vs.size() + 4 == 52) return MSH_TET_52;
+ if(_order == 6 && _vs.size() + 4 == 84) return MSH_TET_84;
+ if(_order == 7 && _vs.size() + 4 == 120) return MSH_TET_120;
+ if(_order == 8 && _vs.size() + 4 == 165) return MSH_TET_165;
+ if(_order == 9 && _vs.size() + 4 == 220) return MSH_TET_220;
+ if(_order == 10 && _vs.size() + 4 == 286) return MSH_TET_286;
return 0;
}
virtual void revert()
{
MVertex *tmp;
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
- switch (getTypeForMSH()) {
- case MSH_TET_20:
- {
- std::vector<MVertex*> inv(16);
- for (int i=0;i<16;i++) inv[i] = _vs[reverseTet20[i+4]-4];
- _vs = inv;
- break;
- }
- case MSH_TET_35:
- {
- std::vector<MVertex*> inv(31);
- for (int i=0;i<31;i++) inv[i] = _vs[reverseTet35[i+4]-4];
- _vs = inv;
- break;
- }
- case MSH_TET_34:
- {
- std::vector<MVertex*> inv(30);
- for (int i=0;i<30;i++) inv[i] = _vs[reverseTet34[i+4]-4];
- _vs = inv;
- break;
- }
- default:
- {
- Msg::Error("Reversion of %d order tetrahedron (type %d) not implemented\n",
- _order, getTypeForMSH());
- break;
- }
- }
+ std::vector<MVertex*> inv(_vs.size());
+ std::vector<int> reverseIndices = _getReverseIndices(_order);
+ for (int i = 0; i<_vs.size(); i++)
+ inv[i] = _vs[reverseIndices[i+4]-4];
+ _vs = inv;
}
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual int getNumEdgesRep();
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 1ad6aaba7c..d10ff18ba3 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -479,29 +479,29 @@ static void reorientTrianglePoints(std::vector<MVertex*> &vtcs, int orientation,
bool swap)
{
int nbPts = vtcs.size();
-
- if(nbPts <= 1) return;
-
- if(nbPts > 3){
- Msg::Error("Interior face nodes reorientation not supported for order > 4");
- return;
- }
-
+ if (nbPts <= 1) return;
std::vector<MVertex*> tmp(nbPts);
-
- // rotation
- // --- interior "principal vertices"
- for (int i = 0; i < 3; i++) tmp[(i + orientation) % 3] = vtcs[i];
-
- // normal swap
- if (swap) {
- // --- interior "principal vertices"
- vtcs[orientation] = tmp[orientation];
- vtcs[(orientation + 1) % 3] = tmp[(orientation + 2) % 3];
- vtcs[(orientation + 2) % 3] = tmp[(orientation + 1) % 3];
+ int interiorOrder = (int)((sqrt(1+8*nbPts)-3)/2);
+ int pos = 0;
+ for (int o = interiorOrder; o>0; o-=3) {
+ if (swap) {
+ tmp[pos] = vtcs[pos];
+ tmp[pos+1] = vtcs[pos+2];
+ tmp[pos+2] = vtcs[pos+1];
+ for (int i = 0; i < 3*(o-1); i++)
+ tmp[pos+3+i] = vtcs[pos+3*o-i-1];
+ } else {
+ for (int i=0; i< 3*o; i++)
+ tmp[pos+i] = vtcs[pos+i];
+ }
+ for (int i = 0; i < 3; i++) {
+ int ri = (i+orientation)%3;
+ vtcs[pos+ri] = tmp[pos+i];
+ for (int j = 0; j < o-1; j++)
+ vtcs[pos+3+(o-1)*ri+j] = tmp[pos+3+(o-1)*i+j];
+ }
+ pos += 3*o;
}
- // no swap
- else vtcs = tmp;
}
// KH: check face orientation wrt element ...
@@ -512,8 +512,11 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
{
fullMatrix<double> points;
int start = 0;
-
switch (nPts){
+ case 0 :
+ case 1 :
+ // do nothing (e.g. for 2nd order tri faces or for quad faces)
+ break;
case 2 :
points = polynomialBases::find(MSH_TRI_10)->points;
start = 9;
@@ -526,8 +529,28 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
points = polynomialBases::find(MSH_TRI_21)->points;
start = 15;
break;
+ case 5 :
+ points = polynomialBases::find(MSH_TRI_28)->points;
+ start = 18;
+ break;
+ case 6 :
+ points = polynomialBases::find(MSH_TRI_36)->points;
+ start = 21;
+ break;
+ case 7 :
+ points = polynomialBases::find(MSH_TRI_45)->points;
+ start = 24;
+ break;
+ case 8 :
+ points = polynomialBases::find(MSH_TRI_55)->points;
+ start = 27;
+ break;
+ case 9 :
+ points = polynomialBases::find(MSH_TRI_66)->points;
+ start = 30;
+ break;
default :
- // do nothing (e.g. for 2nd order tri faces or for quad faces)
+ Msg::Error("getFaceVertices not implemented for order %i\n",nPts +1);
break;
}
@@ -607,18 +630,36 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
int start = 0;
switch (nPts){
+ case 0:
+ case 1:
+ case 2:
+ // done: return!
+ return;
case 3 :
points = polynomialBases::find(MSH_TET_35)->points;
- start = 34;
break;
case 4 :
points = polynomialBases::find(MSH_TET_56)->points;
- start = 52;
+ break;
+ case 5 :
+ points = polynomialBases::find(MSH_TET_84)->points;
+ break;
+ case 6 :
+ points = polynomialBases::find(MSH_TET_120)->points;
+ break;
+ case 7 :
+ points = polynomialBases::find(MSH_TET_165)->points;
+ break;
+ case 8 :
+ points = polynomialBases::find(MSH_TET_220)->points;
+ break;
+ case 9 :
+ points = polynomialBases::find(MSH_TET_286)->points;
break;
default :
- // done: return!
- return;
+ Msg::Error("getRegionVertices not implemented for order %i\n", nPts+1);
}
+ start = ((nPts+2)*(nPts+3)*(nPts+4)-(nPts-2)*(nPts-1)*(nPts))/6;
for(int k = start; k < points.size1(); k++){
MVertex *v;
@@ -782,8 +823,9 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
ve.insert(ve.end(), vr.begin(), vr.end());
MTetrahedron* n = new MTetrahedronN(t->getVertex(0), t->getVertex(1),
t->getVertex(2), t->getVertex(3), ve, nPts + 1);
- if (!mappingIsInvertible(n))
+ if (!mappingIsInvertible(n)){
Msg::Warning("Found invalid curved volume element (# %d in list)", i);
+ }
tetrahedra2.push_back(n);
}
delete t;
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 85fbc1adb9..629060dedd 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -1159,24 +1159,54 @@ const polynomialBasis *polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(5, false);
generate3dFaceClosure(F.closures, 5);
break;
+ case MSH_TET_74 :
+ F.numFaces = 4;
+ F.monomials = generatePascalSerendipityTetrahedron(6);
+ F.points = gmshGeneratePointsTetrahedron(6, true);
+ generate3dFaceClosure(F.closures, 6);
+ break;
case MSH_TET_84 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(6);
F.points = gmshGeneratePointsTetrahedron(6, false);
generate3dFaceClosure(F.closures, 6);
break;
+ case MSH_TET_100 :
+ F.numFaces = 4;
+ F.monomials = generatePascalSerendipityTetrahedron(7);
+ F.points = gmshGeneratePointsTetrahedron(7, true);
+ generate3dFaceClosure(F.closures, 7);
+ break;
case MSH_TET_120 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(7);
F.points = gmshGeneratePointsTetrahedron(7, false);
generate3dFaceClosure(F.closures, 7);
break;
+ case MSH_TET_130 :
+ F.numFaces = 4;
+ F.monomials = generatePascalSerendipityTetrahedron(8);
+ F.points = gmshGeneratePointsTetrahedron(8, true);
+ generate3dFaceClosure(F.closures, 8);
+ break;
+ case MSH_TET_164 :
+ F.numFaces = 4;
+ F.monomials = generatePascalSerendipityTetrahedron(9);
+ F.points = gmshGeneratePointsTetrahedron(9, true);
+ generate3dFaceClosure(F.closures, 9);
+ break;
case MSH_TET_165 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(8);
F.points = gmshGeneratePointsTetrahedron(8, false);
generate3dFaceClosure(F.closures, 8);
break;
+ case MSH_TET_202 :
+ F.numFaces = 4;
+ F.monomials = generatePascalSerendipityTetrahedron(10);
+ F.points = gmshGeneratePointsTetrahedron(10, true);
+ generate3dFaceClosure(F.closures, 10);
+ break;
case MSH_TET_220 :
F.numFaces = 4;
F.monomials = generatePascalTetrahedron(9);
--
GitLab