diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 7a89cdff8c44fa5cf6574a4bb57640f3d1db82c2..b07ae2fdc2f4a4dae308232a098a3b8e7dd84fa8 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1605,9 +1605,6 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_PRI_6: return new MPrism(v, num, part);
case MSH_PRI_15: return new MPrism15(v, num, part);
case MSH_PRI_18: return new MPrism18(v, num, part);
- case MSH_PYR_5: return new MPyramid(v, num, part);
- case MSH_PYR_13: return new MPyramid13(v, num, part);
- case MSH_PYR_14: return new MPyramid14(v, num, part);
case MSH_TET_20: return new MTetrahedronN(v, 3, num, part);
case MSH_TET_34: return new MTetrahedronN(v, 3, num, part);
case MSH_TET_35: return new MTetrahedronN(v, 4, num, part);
@@ -1631,6 +1628,9 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_LIN_SUB: return new MSubLine(v, num, part, owner, parent);
case MSH_TRI_SUB: return new MSubTriangle(v, num, part, owner, parent);
case MSH_TET_SUB: return new MSubTetrahedron(v, num, part, owner, parent);
+ case MSH_PYR_5: return new MPyramid(v, num, part);
+ case MSH_PYR_13: return new MPyramidN(v, 2, num, part);
+ case MSH_PYR_14: return new MPyramidN(v, 2, num, part);
case MSH_PYR_30: return new MPyramidN(v, 3, num, part);
case MSH_PYR_55: return new MPyramidN(v, 4, num, part);
case MSH_PYR_91: return new MPyramidN(v, 5, num, part);
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index ffae41f7d1cfd15a0983c2339711bb155be5446d..e59ff396e4458d69ba8705ea00aa4c0ab3450b4c 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -193,101 +193,6 @@ class MPyramid : public MElement {
* 1--------8-------2
*
*/
-class MPyramid13 : public MPyramid {
- protected:
- MVertex *_vs[8];
- public :
- MPyramid13(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
- MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, MVertex *v9,
- MVertex *v10, MVertex *v11, MVertex *v12, int num=0, int part=0)
- : MPyramid(v0, v1, v2, v3, v4, num, part)
- {
- _vs[0] = v5; _vs[1] = v6; _vs[2] = v7; _vs[3] = v8; _vs[4] = v9;
- _vs[5] = v10; _vs[6] = v11; _vs[7] = v12;
- for(int i = 0; i < 8; i++) _vs[i]->setPolynomialOrder(2);
- }
- MPyramid13(const std::vector<MVertex*> &v, int num=0, int part=0)
- : MPyramid(v, num, part)
- {
- for(int i = 0; i < 8; i++) _vs[i] = v[5 + i];
- for(int i = 0; i < 8; i++) _vs[i]->setPolynomialOrder(2);
- }
- ~MPyramid13(){}
- virtual int getPolynomialOrder() const { return 2; }
- virtual int getNumVertices() const { return 13; }
- virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
- virtual int getNumEdgeVertices() const { return 8; }
- virtual int getNumEdgesRep(){ return 16; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[16][2] = {
- {0, 5}, {5, 1},
- {0, 6}, {6, 3},
- {0, 7}, {7, 4},
- {1, 8}, {8, 2},
- {1, 9}, {9, 4},
- {2, 10}, {10, 3},
- {2, 11}, {11, 4},
- {3, 12}, {12, 4}
- };
- static const int f[8] = {0, 1, 1, 2, 0, 3, 2, 3};
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
- }
- virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
- {
- v.resize(3);
- MPyramid::_getEdgeVertices(num, v);
- v[2] = _vs[num];
- }
- virtual int getNumFacesRep(){ return 22; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[22][3] = {
- {0, 5, 7}, {1, 9, 5}, {4, 7, 9}, {5, 9, 7},
- {3, 6, 12}, {0, 7, 6}, {4, 12, 7}, {6, 7, 12},
- {1, 8, 9}, {2, 11, 8}, {4, 9, 11}, {8, 11, 9},
- {2, 10, 11}, {3, 12, 10}, {4, 11, 12}, {10, 12, 11},
- {0, 6, 5}, {3, 10, 6}, {2, 8, 10}, {1, 5, 8}, {5, 6, 10}, {5, 10, 8}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
- virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
- {
- v.resize((num < 4) ? 6 : 8);
- MPyramid::_getFaceVertices(num, v);
- static const int f[4][3] = {
- {0, 4, 2},
- {1, 2, 7},
- {3, 6, 4},
- {5, 7, 6}
- };
- if(num < 4) {
- v[3] = _vs[f[num][0]];
- v[4] = _vs[f[num][1]];
- v[5] = _vs[f[num][2]];
- }
- else {
- v[4] = _vs[1];
- v[5] = _vs[5];
- v[6] = _vs[3];
- v[7] = _vs[0];
- }
- }
- virtual int getTypeForMSH() const { return MSH_PYR_13; }
- virtual void revert()
- {
- MVertex *tmp;
- tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
- tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp;
- tmp = _vs[1]; _vs[1] = _vs[5]; _vs[5] = tmp;
- tmp = _vs[2]; _vs[2] = _vs[6]; _vs[6] = tmp;
- }
- virtual void getNode(int num, double &u, double &v, double &w)
- {
- num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
- }
-};
/*
* MPyramid14
@@ -308,106 +213,6 @@ class MPyramid13 : public MPyramid {
* 1--------8-------2
*
*/
-class MPyramid14 : public MPyramid {
- protected:
- MVertex *_vs[9];
- public :
- MPyramid14(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
- MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, MVertex *v9,
- MVertex *v10, MVertex *v11, MVertex *v12, MVertex *v13,
- int num=0, int part=0)
- : MPyramid(v0, v1, v2, v3, v4, num, part)
- {
- _vs[0] = v5; _vs[1] = v6; _vs[2] = v7; _vs[3] = v8; _vs[4] = v9;
- _vs[5] = v10; _vs[6] = v11; _vs[7] = v12; _vs[8] = v13;
- for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
- }
- MPyramid14(const std::vector<MVertex*> &v, int num=0, int part=0)
- : MPyramid(v, num, part)
- {
- for(int i = 0; i < 9; i++) _vs[i] = v[5 + i];
- for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
- }
- ~MPyramid14(){}
- virtual int getPolynomialOrder() const { return 2; }
- virtual int getNumVertices() const { return 14; }
- virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
- virtual int getNumEdgeVertices() const { return 8; }
- virtual int getNumFaceVertices() const { return 1; }
- virtual int getNumEdgesRep(){ return 16; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[16][2] = {
- {0, 5}, {5, 1},
- {0, 6}, {6, 3},
- {0, 7}, {7, 4},
- {1, 8}, {8, 2},
- {1, 9}, {9, 4},
- {2, 10}, {10, 3},
- {2, 11}, {11, 4},
- {3, 12}, {12, 4}
- };
- static const int f[8] = {0, 1, 1, 2, 0, 3, 2, 3};
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
- }
- virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
- {
- v.resize(3);
- MPyramid::_getEdgeVertices(num, v);
- v[2] = _vs[num];
- }
- virtual int getNumFacesRep(){ return 24; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[24][3] = {
- {0, 5, 7}, {1, 9, 5}, {4, 7, 9}, {5, 9, 7},
- {3, 6, 12}, {0, 7, 6}, {4, 12, 7}, {6, 7, 12},
- {1, 8, 9}, {2, 11, 8}, {4, 9, 11}, {8, 11, 9},
- {2, 10, 11}, {3, 12, 10}, {4, 11, 12}, {10, 12, 11},
- {0, 6, 13}, {0, 13, 5}, {3, 10, 13}, {3, 13, 6},
- {2, 8, 13}, {2, 13, 10}, {1, 5, 13}, {1, 13, 8}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
- virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
- {
- v.resize((num < 4) ? 6 : 9);
- MPyramid::_getFaceVertices(num, v);
- static const int f[4][3] = {
- {0, 4, 2},
- {1, 2, 7},
- {3, 6, 4},
- {5, 7, 6}
- };
- if(num < 4) {
- v[3] = _vs[f[num][0]];
- v[4] = _vs[f[num][1]];
- v[5] = _vs[f[num][2]];
- }
- else {
- v[4] = _vs[1];
- v[5] = _vs[5];
- v[6] = _vs[3];
- v[7] = _vs[0];
- v[8] = _vs[8];
- }
- }
- virtual int getTypeForMSH() const { return MSH_PYR_14; }
- virtual const char *getStringForPOS() const { return "SY2"; }
- virtual void revert()
- {
- MVertex *tmp;
- tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
- tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp;
- tmp = _vs[1]; _vs[1] = _vs[5]; _vs[5] = tmp;
- tmp = _vs[2]; _vs[2] = _vs[6]; _vs[6] = tmp;
- }
- virtual void getNode(int num, double &u, double &v, double &w)
- {
- num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
- }
-};
//------------------------------------------------------------------------------
@@ -483,6 +288,8 @@ class MPyramidN : public MPyramid {
}
virtual int getTypeForMSH() const
{
+ if(_order == 2 && _vs.size() + 5 == 13) return MSH_PYR_13;
+ if(_order == 2 && _vs.size() + 5 == 14) return MSH_PYR_14;
if(_order == 3 && _vs.size() + 5 == 29) return MSH_PYR_29;
if(_order == 3 && _vs.size() + 5 == 30) return MSH_PYR_30;
if(_order == 4 && _vs.size() + 5 == 50) return MSH_PYR_50;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 07e7077804248a5581b791a37c15838d43a5d233..9d6fd10e25d856ab42257909151c1af1f289fea8 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1023,61 +1023,35 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
}
gr->prisms = prisms2;
+/* * * * * * * * * * * * * * * * * PYRAMIDS * * * * * * * * * * * * * * * * * */
+
std::vector<MPyramid*> pyramids2;
+
for(unsigned int i = 0; i < gr->pyramids.size(); i++) {
MPyramid *p = gr->pyramids[i];
std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
- if (nPts == 1) {
- if(incomplete) {
- pyramids2.push_back(new MPyramid13(p->getVertex(0), p->getVertex(1),
- p->getVertex(2), p->getVertex(3),
- p->getVertex(4), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7],
- 0, p->getPartition()));
- }
- else {
- getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
- pyramids2.push_back
- (new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7], vf[0],
- 0, p->getPartition()));
- }
- }
- else {
- getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
- ve.insert(ve.end(), vf.begin(), vf.end());
- // Creating quads to get the internal vertices
- /*
- * P3 : q1 - 21 22 23 24
- * P4 : q1 - 29 30 32 33 35 36 38 39
- * q2 - 31 34 37 40
- * P5 : q1 - 37 40 38 43 46 44 49 52 50 55 58 56
- * q2 - 42 41 48 47 54 53 60 59
- * q3 - 39 45 51 57
- */
-
- vr.reserve((nPts-1)*(nPts)*(2*(nPts-1)+1)/6);
- /*for (int tmp = 0; tmp < (nPts-1)*(nPts)*(2*(nPts-1)+1)/6; tmp++) {
- vr.push_back(0);
- }*/
-
- int verts_lvl3[12] = {37,40,38,43,46,44,49,52,50,55,58,56};
-
- int verts_lvl2[8];
- if (nPts == 4) {
- verts_lvl2[0] = 42; verts_lvl2[1] = 41;
- verts_lvl2[2] = 48; verts_lvl2[3] = 47;
- verts_lvl2[4] = 54; verts_lvl2[5] = 53;
- verts_lvl2[6] = 60; verts_lvl2[7] = 59;
- } else {
- verts_lvl2[0] = 29; verts_lvl2[1] = 30;
- verts_lvl2[2] = 35; verts_lvl2[3] = 36;
- verts_lvl2[4] = 38; verts_lvl2[5] = 39;
- verts_lvl2[6] = 32; verts_lvl2[7] = 33;
- }
+ getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+
+ ve.insert(ve.end(), vf.begin(), vf.end());
+
+ vr.reserve((nPts-1)*(nPts)*(2*(nPts-1)+1)/6);
+
+ int verts_lvl3[12] = {37,40,38,43,46,44,49,52,50,55,58,56};
+
+ int verts_lvl2[8];
+ if (nPts == 4) {
+ verts_lvl2[0] = 42; verts_lvl2[1] = 41;
+ verts_lvl2[2] = 48; verts_lvl2[3] = 47;
+ verts_lvl2[4] = 54; verts_lvl2[5] = 53;
+ verts_lvl2[6] = 60; verts_lvl2[7] = 59;
+ } else {
+ verts_lvl2[0] = 29; verts_lvl2[1] = 30;
+ verts_lvl2[2] = 35; verts_lvl2[3] = 36;
+ verts_lvl2[4] = 38; verts_lvl2[5] = 39;
+ verts_lvl2[6] = 32; verts_lvl2[7] = 33;
+ }
int verts_lvl1[4];
switch(nPts) {
@@ -1101,100 +1075,99 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
break;
}
- for (int q = 0; q < nPts - 1; q++) {
- std::vector<MVertex*> vq, veq;
- vq.push_back(ve[2*nPts + q]);
- vq.push_back(ve[4*nPts + q]);
- vq.push_back(ve[6*nPts + q]);
- vq.push_back(ve[7*nPts + q]);
-
- //int triverts = nPts*(nPts-1)/2;
-
- if (nPts-q == 4)
- for (int f = 0; f < 12; f++)
- veq.push_back(ve[verts_lvl3[f]-5]);
- else if (nPts-q == 3)
- for (int f = 0; f < 8; f++)
- veq.push_back(ve[verts_lvl2[f]-5]);
- else if (nPts-q == 2)
- for (int f = 0; f < 4; f++)
- veq.push_back(ve[verts_lvl1[f]-5]);
-
- if (nPts-q == 2) {
- MQuadrangle8 incpl2(vq[0], vq[1], vq[2], vq[3],
- veq[0], veq[1], veq[2], veq[3]);
- SPoint3 pointz;
- incpl2.pnt(0,0,0,pointz);
- MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
+ for (int q = 0; q < nPts - 1; q++) {
+ std::vector<MVertex*> vq, veq;
+ vq.push_back(ve[2*nPts + q]);
+ vq.push_back(ve[4*nPts + q]);
+ vq.push_back(ve[6*nPts + q]);
+ vq.push_back(ve[7*nPts + q]);
+
+ if (nPts-q == 4)
+ for (int f = 0; f < 12; f++)
+ veq.push_back(ve[verts_lvl3[f]-5]);
+ else if (nPts-q == 3)
+ for (int f = 0; f < 8; f++)
+ veq.push_back(ve[verts_lvl2[f]-5]);
+ else if (nPts-q == 2)
+ for (int f = 0; f < 4; f++)
+ veq.push_back(ve[verts_lvl1[f]-5]);
+
+ if (nPts-q == 2) {
+ MQuadrangle8 incpl2(vq[0], vq[1], vq[2], vq[3],
+ veq[0], veq[1], veq[2], veq[3]);
+ SPoint3 pointz;
+ incpl2.pnt(0,0,0,pointz);
+ MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
+ gr->mesh_vertices.push_back(v);
+ std::vector<MVertex*>::iterator cursor = vr.begin();
+ cursor += nPts == 2 ? 0 : 4;
+ vr.insert(cursor, v);
+ }
+ else if (nPts-q == 3) {
+
+ MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 3);
+ int offsets[4] = {nPts == 4 ? 7 : 0,
+ nPts == 4 ? 9 : 1,
+ nPts == 4 ? 11 : 2,
+ nPts == 4 ? 12 : 3};
+ double quad_v [4][2] = {{-1.0/3.0, -1.0/3.0},
+ { 1.0/3.0, -1.0/3.0},
+ { 1.0/3.0, 1.0/3.0},
+ {-1.0/3.0, 1.0/3.0}};
+ SPoint3 pointz;
+ for (int k = 0; k<4; k++) {
+ incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
+ MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
gr->mesh_vertices.push_back(v);
std::vector<MVertex*>::iterator cursor = vr.begin();
- cursor += nPts == 2 ? 0 : 4;
+ cursor += offsets[k];
vr.insert(cursor, v);
}
- else if (nPts-q == 3) {
-
- MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 3);
- int offsets[4] = {nPts == 4 ? 7 : 0,
- nPts == 4 ? 9 : 1,
- nPts == 4 ? 11 : 2,
- nPts == 4 ? 12 : 3};
- double quad_v [4][2] = {{-1.0/3.0, -1.0/3.0},
- { 1.0/3.0, -1.0/3.0},
- { 1.0/3.0, 1.0/3.0},
- {-1.0/3.0, 1.0/3.0}};
- SPoint3 pointz;
- for (int k = 0; k<4; k++) {
- incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
- MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
- gr->mesh_vertices.push_back(v);
- std::vector<MVertex*>::iterator cursor = vr.begin();
- cursor += offsets[k];
- vr.insert(cursor, v);
- }
- }
- else if (nPts-q == 4) {
- MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 4);
- int offsets[9] = {0, 1, 2, 3, 5, 8, 10, 6, 13};
- double quad_v [9][2] = {
- { -0.5, -0.5},
- { 0.5, -0.5},
- { 0.5, 0.5},
- { -0.5, 0.5},
- { 0.0, -0.5},
- { 0.5, 0.0},
- { 0.0, 0.5},
- { -0.5, 0.0},
- { 0.0, 0.0}
- };
- SPoint3 pointz;
- for (int k = 0; k<9; k++) {
- incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
- MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
- gr->mesh_vertices.push_back(v);
- std::vector<MVertex*>::iterator cursor = vr.begin();
- cursor += offsets[k];
- vr.insert(cursor, v);
- }
-
-
+ }
+ else if (nPts-q == 4) {
+ MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 4);
+ int offsets[9] = {0, 1, 2, 3, 5, 8, 10, 6, 13};
+ double quad_v [9][2] = {
+ { -0.5, -0.5},
+ { 0.5, -0.5},
+ { 0.5, 0.5},
+ { -0.5, 0.5},
+ { 0.0, -0.5},
+ { 0.5, 0.0},
+ { 0.0, 0.5},
+ { -0.5, 0.0},
+ { 0.0, 0.0}
+ };
+ SPoint3 pointz;
+ for (int k = 0; k<9; k++) {
+ incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
+ MVertex *v = new MVertex(pointz.x(), pointz.y(), pointz.z(), gr);
+ gr->mesh_vertices.push_back(v);
+ std::vector<MVertex*>::iterator cursor = vr.begin();
+ cursor += offsets[k];
+ vr.insert(cursor, v);
}
-
}
- ve.insert(ve.end(), vr.begin(), vr.end());
- MPyramid *n = new MPyramidN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), ve, nPts + 1,
- 0, p->getPartition());
- pyramids2.push_back(n);
- SPoint3 test_pnt;
- n->pnt(-1,-1,0, test_pnt);
+
}
+ ve.insert(ve.end(), vr.begin(), vr.end());
+ MPyramid *n = new MPyramidN(p->getVertex(0), p->getVertex(1),
+ p->getVertex(2), p->getVertex(3),
+ p->getVertex(4), ve, nPts + 1,
+ 0, p->getPartition());
+ pyramids2.push_back(n);
+ SPoint3 test_pnt;
+ n->pnt(-1,-1,0, test_pnt);
+
delete p;
}
gr->pyramids = pyramids2;
gr->deleteVertexArrays();
}
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+
template<class T>
static void setFirstOrder(GEntity *e, std::vector<T*> &elements, bool onlyVisible)
{