From f248e026e7a4da795e1ca4454c17be6d93068946 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 6 Sep 2013 13:48:26 +0000
Subject: [PATCH] use serendipity tetrahedron and haxahedron to compute high
order interior vertices (instead of "unusual"... and now inexistant
face-saturated incomplete element)
clean up high order routines by refactoring high-order tet/hex/prism/pyram generation code
---
Fltk/highOrderToolsWindow.cpp | 4 +
Mesh/HighOrder.cpp | 488 ++++++++++++++++++----------------
2 files changed, 265 insertions(+), 227 deletions(-)
diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index 141a39ab58..43c6093393 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -36,10 +36,12 @@ static void change_completeness_cb(Fl_Widget *w, void *data)
highOrderToolsWindow *o = FlGui::instance()->highordertools;
bool onlyVisible = (bool)o->butt[1]->value();
if (!o->complete){
+ // BOF BOF BOF -- CG
SetHighOrderComplete(GModel::current(), onlyVisible);
o->complete = 1;
}
else if (o->complete){
+ // BOF BOF BOF -- CG
SetHighOrderInComplete(GModel::current(), onlyVisible);
o->complete = 0;
}
@@ -61,6 +63,7 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
else
SetOrderN(GModel::current(), order, linear, incomplete, onlyVisible);
+ /*
distanceFromMeshToGeometry_t dist;
computeDistanceFromMeshToGeometry (GModel::current(), dist);
for (std::map<GEntity*, double> ::iterator it = dist.d2.begin();
@@ -68,6 +71,7 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
printf("GEntity %d of dim %d : dist %12.5E\n",
it->first->tag(), it->first->dim(), it->second);
}
+ */
CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
drawContext::global()->draw();
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 8de07ca289..69d6b8c09b 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -644,11 +644,9 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
const double t1 = points(k, 0);
const double t2 = points(k, 1);
const double t3 = points(k, 2);
- //printf("inside getRegionVertices, point is %g %g %g\n", t1, t2, t3);
SPoint3 pos;
incomplete->pnt(t1, t2, t3, pos);
v = new MVertex(pos.x(), pos.y(), pos.z(), gr);
- //printf("inside getRegionVertices, vertex is %g %g %g\n", v->x(), v->y(), v->z());
gr->mesh_vertices.push_back(v);
vr.push_back(v);
}
@@ -684,18 +682,14 @@ static MTriangle *setHighOrder(MTriangle *t, GFace *gf,
ve[0], ve[1], ve[2], 0, t->getPartition());
}
else{
- if(incomplete){
- return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
- ve, nPts + 1, 0, t->getPartition());
- }
- else{
+ if(!incomplete){
MTriangleN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
ve, nPts + 1, 0, t->getPartition());
getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
- return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
- ve, nPts + 1, 0, t->getPartition());
}
+ return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve, nPts + 1, 0, t->getPartition());
}
}
@@ -743,8 +737,8 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
std::vector<MTriangle*> triangles2;
for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
- MTriangle *tNew = setHighOrder(t, gf,edgeVertices, faceVertices, linear, incomplete,
- nPts);
+ MTriangle *tNew = setHighOrder(t, gf, edgeVertices, faceVertices, linear,
+ incomplete, nPts);
triangles2.push_back(tNew);
delete t;
}
@@ -761,258 +755,298 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
gf->deleteVertexArrays();
}
-static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
- faceContainer &faceVertices, bool linear, bool incomplete,
- int nPts = 1)
+static MTetrahedron *setHighOrder(MTetrahedron *t, GRegion *gr,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear, bool incomplete, int nPts)
{
- std::vector<MTetrahedron*> tetrahedra2;
- for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
- MTetrahedron *t = gr->tetrahedra[i];
- std::vector<MVertex*> ve, vf, vr;
- getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts);
- if(nPts == 1){
- tetrahedra2.push_back
- (new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
- t->getVertex(3), ve[0], ve[1], ve[2], ve[3], ve[4], ve[5],
- 0, t->getPartition()));
- }
- else{
- getFaceVertices(gr, t, vf, faceVertices, edgeVertices, linear, nPts);
- ve.insert(ve.end(), vf.begin(), vf.end());
+ std::vector<MVertex*> ve, vf, vr;
+ getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts);
+ if(nPts == 1){
+ return new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ t->getVertex(3), ve[0], ve[1], ve[2], ve[3], ve[4], ve[5],
+ 0, t->getPartition());
+ }
+ else{
+ if(!incomplete){
+ // create serendipity element to place internal vertices (we used to
+ // saturate face vertices also, but the corresponding function spaces do
+ // not exist anymore, and there is no theoretical justification for doing
+ // it either way)
MTetrahedronN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
t->getVertex(3), ve, nPts + 1, 0, t->getPartition());
+ getFaceVertices(gr, t, vf, faceVertices, edgeVertices, linear, nPts);
+ ve.insert(ve.end(), vf.begin(), vf.end());
getRegionVertices(gr, &incpl, t, vr, linear, nPts);
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,
- 0, t->getPartition());
- tetrahedra2.push_back(n);
}
- delete t;
+ return new MTetrahedronN(t->getVertex(0), t->getVertex(1),
+ t->getVertex(2), t->getVertex(3), ve, nPts + 1,
+ 0, t->getPartition());
}
- gr->tetrahedra = tetrahedra2;
+}
- std::vector<MHexahedron*> hexahedra2;
- for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
- MHexahedron *h = gr->hexahedra[i];
- std::vector<MVertex*> ve, vf, vr;
- getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts);
+static MHexahedron *setHighOrder(MHexahedron *h, GRegion *gr,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear, bool incomplete, int nPts)
+{
+ std::vector<MVertex*> ve, vf, vr;
+ getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts);
+ if(incomplete){
if(nPts == 1){
- if(incomplete){
- hexahedra2.push_back
- (new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2),
- h->getVertex(3), h->getVertex(4), h->getVertex(5),
- h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
- ve[11], 0, h->getPartition()));
- }
- else{
- getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
- SPoint3 pc = h->barycenter();
- MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
- gr->mesh_vertices.push_back(v);
- hexahedra2.push_back
- (new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2),
- h->getVertex(3), h->getVertex(4), h->getVertex(5),
- h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
- ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
- ve[11], vf[0], vf[1], vf[2], vf[3], vf[4], vf[5], v,
- 0, h->getPartition()));
- }
+ return new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
+ ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
+ ve[11], 0, h->getPartition());
}
- else {
+ else{
+ return new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve, nPts + 1, 0,
+ h->getPartition());
+ }
+ }
+ else{
+ if(nPts == 1){
getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
- ve.insert(ve.end(), vf.begin(), vf.end());
+ SPoint3 pc = h->barycenter();
+ MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
+ gr->mesh_vertices.push_back(v);
+ return new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2],
+ ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
+ ve[11], vf[0], vf[1], vf[2], vf[3], vf[4], vf[5], v,
+ 0, h->getPartition());
+ }
+ else {
+ // create serendipity element to place internal vertices (we used to
+ // saturate face vertices also, but the corresponding function spaces do
+ // not exist anymore, and there is no theoretical justification for doing
+ // it either way)
MHexahedronN incpl(h->getVertex(0), h->getVertex(1), h->getVertex(2),
h->getVertex(3), h->getVertex(4), h->getVertex(5),
- h->getVertex(6), h->getVertex(7), ve, nPts + 1, 0, h->getPartition());
+ h->getVertex(6), h->getVertex(7), ve, nPts + 1, 0,
+ h->getPartition());
+ getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
+ ve.insert(ve.end(), vf.begin(), vf.end());
getRegionVertices(gr, &incpl, h, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end());
- MHexahedron *n = new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
- h->getVertex(3), h->getVertex(4), h->getVertex(5),
- h->getVertex(6), h->getVertex(7), ve, nPts + 1,
- 0, h->getPartition());
- hexahedra2.push_back(n);
+ return new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve, nPts + 1,
+ 0, h->getPartition());
}
- delete h;
}
- gr->hexahedra = hexahedra2;
+}
- std::vector<MPrism*> prisms2;
- for(unsigned int i = 0; i < gr->prisms.size(); i++){
- MPrism *p = gr->prisms[i];
- std::vector<MVertex*> ve, vf, vr;
- getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
- if(incomplete){
- prisms2.push_back
- (new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), p->getVertex(5),
- ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
- 0, p->getPartition()));
+static MPrism *setHighOrder(MPrism *p, GRegion *gr,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear, bool incomplete, int nPts)
+{
+ std::vector<MVertex*> ve, vf, vr;
+ getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
+ if(incomplete){
+ if(nPts == 1){
+ return new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+ p->getVertex(3), p->getVertex(4), p->getVertex(5),
+ ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
+ 0, p->getPartition());
}
else{
- if (nPts == 1) {
- getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
- prisms2.push_back
- (new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), p->getVertex(5),
- ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
- vf[0], vf[1], vf[2],
- 0, p->getPartition()));
- }
- else {
- Msg::Error("PrismN generation not implemented");
- /*
- getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
- ve.insert(ve.end(), vf.begin(), vf.end());
- MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2), p->getVertex(3),
- p->getVertex(4), p->getVertex(5), ve, nPts + 1);
- getRegionVertices(gr, &incpl, p, vr, linear, nPts);
- if (nPts == 0) {
- printf("Screwed\n");
- }
- ve.insert(ve.end(), vr.begin(), vr.end());
- MPrism* n = new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- p->getVertex(3), p->getVertex(4), p->getVertex(5),
- ve, nPts+1);
- if (!mappingIsInvertible(n))
- Msg::Warning("Found invalid curved volume element (# %d in list)", i);
- prisms2.push_back(n);
- */
- }
+ Msg::Error("PrismN generation not implemented");
+ return 0;
}
- delete p;
}
- gr->prisms = prisms2;
-
- 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);
- 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{
+ if (nPts == 1) {
+ getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+ return new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+ p->getVertex(3), p->getVertex(4), p->getVertex(5),
+ ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
+ vf[0], vf[1], vf[2],
+ 0, p->getPartition());
}
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;
+ Msg::Error("PrismN generation not implemented");
+ return 0;
}
- int verts_lvl1[4];
- switch(nPts) {
- case(4):
- verts_lvl1[0] = 39;
- verts_lvl1[1] = 45;
- verts_lvl1[2] = 51;
- verts_lvl1[3] = 57;
- break;
- case(3):
- verts_lvl1[0] = 31;
- verts_lvl1[1] = 37;
- verts_lvl1[2] = 40;
- verts_lvl1[3] = 34;
- break;
- case(2):
- verts_lvl1[0] = 21;
- verts_lvl1[1] = 23;
- verts_lvl1[2] = 24;
- verts_lvl1[3] = 22;
- break;
+ }
+}
+
+static MPyramid *setHighOrder(MPyramid *p, GRegion *gr,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear, bool incomplete, int nPts)
+{
+ std::vector<MVertex*> ve, vf, vr;
+ getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
+ 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) {
+ case(4):
+ verts_lvl1[0] = 39;
+ verts_lvl1[1] = 45;
+ verts_lvl1[2] = 51;
+ verts_lvl1[3] = 57;
+ break;
+ case(3):
+ verts_lvl1[0] = 31;
+ verts_lvl1[1] = 37;
+ verts_lvl1[2] = 40;
+ verts_lvl1[3] = 34;
+ break;
+ case(2):
+ verts_lvl1[0] = 21;
+ verts_lvl1[1] = 23;
+ verts_lvl1[2] = 24;
+ verts_lvl1[3] = 22;
+ 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]);
+
+ 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);
}
- 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);
+ 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());
+ return new MPyramidN(p->getVertex(0), p->getVertex(1),
+ p->getVertex(2), p->getVertex(3),
+ p->getVertex(4), ve, nPts + 1,
+ 0, p->getPartition());
+}
+
+static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
+ faceContainer &faceVertices, bool linear, bool incomplete,
+ int nPts = 1)
+{
+ std::vector<MTetrahedron*> tetrahedra2;
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
+ MTetrahedron *t = gr->tetrahedra[i];
+ MTetrahedron *tNew = setHighOrder(t, gr, edgeVertices, faceVertices, linear,
+ incomplete, nPts);
+ tetrahedra2.push_back(tNew);
+ delete t;
+ }
+ gr->tetrahedra = tetrahedra2;
+
+ std::vector<MHexahedron*> hexahedra2;
+ for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
+ MHexahedron *h = gr->hexahedra[i];
+ MHexahedron *hNew = setHighOrder(h, gr, edgeVertices, faceVertices, linear,
+ incomplete, nPts);
+ hexahedra2.push_back(hNew);
+ delete h;
+ }
+ gr->hexahedra = hexahedra2;
+
+ std::vector<MPrism*> prisms2;
+ for(unsigned int i = 0; i < gr->prisms.size(); i++){
+ MPrism *p = gr->prisms[i];
+ MPrism *pNew = setHighOrder(p, gr, edgeVertices, faceVertices, linear,
+ incomplete, nPts);
+ prisms2.push_back(pNew);
+ delete p;
+ }
+ gr->prisms = prisms2;
+
+ std::vector<MPyramid*> pyramids2;
+ for(unsigned int i = 0; i < gr->pyramids.size(); i++) {
+ MPyramid *p = gr->pyramids[i];
+ MPyramid *pNew = setHighOrder(p, gr, edgeVertices, faceVertices, linear,
+ incomplete, nPts);
+ pyramids2.push_back(pNew);
delete p;
}
gr->pyramids = pyramids2;
+
gr->deleteVertexArrays();
}
--
GitLab