Skip to content
Snippets Groups Projects
Commit f248e026 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

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
parent c3fa31d2
No related branches found
No related tags found
No related merge requests found
...@@ -36,10 +36,12 @@ static void change_completeness_cb(Fl_Widget *w, void *data) ...@@ -36,10 +36,12 @@ static void change_completeness_cb(Fl_Widget *w, void *data)
highOrderToolsWindow *o = FlGui::instance()->highordertools; highOrderToolsWindow *o = FlGui::instance()->highordertools;
bool onlyVisible = (bool)o->butt[1]->value(); bool onlyVisible = (bool)o->butt[1]->value();
if (!o->complete){ if (!o->complete){
// BOF BOF BOF -- CG
SetHighOrderComplete(GModel::current(), onlyVisible); SetHighOrderComplete(GModel::current(), onlyVisible);
o->complete = 1; o->complete = 1;
} }
else if (o->complete){ else if (o->complete){
// BOF BOF BOF -- CG
SetHighOrderInComplete(GModel::current(), onlyVisible); SetHighOrderInComplete(GModel::current(), onlyVisible);
o->complete = 0; o->complete = 0;
} }
...@@ -61,6 +63,7 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data) ...@@ -61,6 +63,7 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
else else
SetOrderN(GModel::current(), order, linear, incomplete, onlyVisible); SetOrderN(GModel::current(), order, linear, incomplete, onlyVisible);
/*
distanceFromMeshToGeometry_t dist; distanceFromMeshToGeometry_t dist;
computeDistanceFromMeshToGeometry (GModel::current(), dist); computeDistanceFromMeshToGeometry (GModel::current(), dist);
for (std::map<GEntity*, double> ::iterator it = dist.d2.begin(); for (std::map<GEntity*, double> ::iterator it = dist.d2.begin();
...@@ -68,6 +71,7 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data) ...@@ -68,6 +71,7 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
printf("GEntity %d of dim %d : dist %12.5E\n", printf("GEntity %d of dim %d : dist %12.5E\n",
it->first->tag(), it->first->dim(), it->second); it->first->tag(), it->first->dim(), it->second);
} }
*/
CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME); CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
drawContext::global()->draw(); drawContext::global()->draw();
......
...@@ -644,11 +644,9 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele, ...@@ -644,11 +644,9 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
const double t1 = points(k, 0); const double t1 = points(k, 0);
const double t2 = points(k, 1); const double t2 = points(k, 1);
const double t3 = points(k, 2); const double t3 = points(k, 2);
//printf("inside getRegionVertices, point is %g %g %g\n", t1, t2, t3);
SPoint3 pos; SPoint3 pos;
incomplete->pnt(t1, t2, t3, pos); incomplete->pnt(t1, t2, t3, pos);
v = new MVertex(pos.x(), pos.y(), pos.z(), gr); 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); gr->mesh_vertices.push_back(v);
vr.push_back(v); vr.push_back(v);
} }
...@@ -684,20 +682,16 @@ static MTriangle *setHighOrder(MTriangle *t, GFace *gf, ...@@ -684,20 +682,16 @@ static MTriangle *setHighOrder(MTriangle *t, GFace *gf,
ve[0], ve[1], ve[2], 0, t->getPartition()); ve[0], ve[1], ve[2], 0, t->getPartition());
} }
else{ else{
if(incomplete){ if(!incomplete){
return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
ve, nPts + 1, 0, t->getPartition());
}
else{
MTriangleN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2), MTriangleN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
ve, nPts + 1, 0, t->getPartition()); ve, nPts + 1, 0, t->getPartition());
getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts); getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end()); ve.insert(ve.end(), vf.begin(), vf.end());
}
return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2), return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
ve, nPts + 1, 0, t->getPartition()); ve, nPts + 1, 0, t->getPartition());
} }
} }
}
static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf, static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
edgeContainer &edgeVertices, edgeContainer &edgeVertices,
...@@ -743,8 +737,8 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices, ...@@ -743,8 +737,8 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
std::vector<MTriangle*> triangles2; std::vector<MTriangle*> triangles2;
for(unsigned int i = 0; i < gf->triangles.size(); i++){ for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i]; MTriangle *t = gf->triangles[i];
MTriangle *tNew = setHighOrder(t, gf,edgeVertices, faceVertices, linear, incomplete, MTriangle *tNew = setHighOrder(t, gf, edgeVertices, faceVertices, linear,
nPts); incomplete, nPts);
triangles2.push_back(tNew); triangles2.push_back(tNew);
delete t; delete t;
} }
...@@ -761,133 +755,133 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices, ...@@ -761,133 +755,133 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
gf->deleteVertexArrays(); gf->deleteVertexArrays();
} }
static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices, static MTetrahedron *setHighOrder(MTetrahedron *t, GRegion *gr,
faceContainer &faceVertices, bool linear, bool incomplete, edgeContainer &edgeVertices,
int nPts = 1) 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; std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts); getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts);
if(nPts == 1){ if(nPts == 1){
tetrahedra2.push_back return new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
(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], t->getVertex(3), ve[0], ve[1], ve[2], ve[3], ve[4], ve[5],
0, t->getPartition())); 0, t->getPartition());
} }
else{ else{
getFaceVertices(gr, t, vf, faceVertices, edgeVertices, linear, nPts); if(!incomplete){
ve.insert(ve.end(), vf.begin(), vf.end()); // 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), MTetrahedronN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
t->getVertex(3), ve, nPts + 1, 0, t->getPartition()); 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); getRegionVertices(gr, &incpl, t, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end()); ve.insert(ve.end(), vr.begin(), vr.end());
MTetrahedron *n = new MTetrahedronN(t->getVertex(0), t->getVertex(1), }
return new MTetrahedronN(t->getVertex(0), t->getVertex(1),
t->getVertex(2), t->getVertex(3), ve, nPts + 1, t->getVertex(2), t->getVertex(3), ve, nPts + 1,
0, t->getPartition()); 0, t->getPartition());
tetrahedra2.push_back(n);
} }
delete t;
} }
gr->tetrahedra = tetrahedra2;
std::vector<MHexahedron*> hexahedra2; static MHexahedron *setHighOrder(MHexahedron *h, GRegion *gr,
for(unsigned int i = 0; i < gr->hexahedra.size(); i++){ edgeContainer &edgeVertices,
MHexahedron *h = gr->hexahedra[i]; faceContainer &faceVertices,
bool linear, bool incomplete, int nPts)
{
std::vector<MVertex*> ve, vf, vr; std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts); getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts);
if(nPts == 1){
if(incomplete){ if(incomplete){
hexahedra2.push_back if(nPts == 1){
(new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2), return new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2),
h->getVertex(3), h->getVertex(4), h->getVertex(5), h->getVertex(3), h->getVertex(4), h->getVertex(5),
h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2], 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[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10],
ve[11], 0, h->getPartition())); ve[11], 0, h->getPartition());
}
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{ else{
if(nPts == 1){
getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts); getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
SPoint3 pc = h->barycenter(); SPoint3 pc = h->barycenter();
MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr); MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
gr->mesh_vertices.push_back(v); gr->mesh_vertices.push_back(v);
hexahedra2.push_back return new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2),
(new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2),
h->getVertex(3), h->getVertex(4), h->getVertex(5), h->getVertex(3), h->getVertex(4), h->getVertex(5),
h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2], 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[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, ve[11], vf[0], vf[1], vf[2], vf[3], vf[4], vf[5], v,
0, h->getPartition())); 0, h->getPartition());
}
} }
else { else {
getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts); // create serendipity element to place internal vertices (we used to
ve.insert(ve.end(), vf.begin(), vf.end()); // 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), MHexahedronN incpl(h->getVertex(0), h->getVertex(1), h->getVertex(2),
h->getVertex(3), h->getVertex(4), h->getVertex(5), 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); getRegionVertices(gr, &incpl, h, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end()); ve.insert(ve.end(), vr.begin(), vr.end());
MHexahedron *n = new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2), return new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
h->getVertex(3), h->getVertex(4), h->getVertex(5), h->getVertex(3), h->getVertex(4), h->getVertex(5),
h->getVertex(6), h->getVertex(7), ve, nPts + 1, h->getVertex(6), h->getVertex(7), ve, nPts + 1,
0, h->getPartition()); 0, h->getPartition());
hexahedra2.push_back(n);
} }
delete h;
} }
gr->hexahedra = hexahedra2; }
std::vector<MPrism*> prisms2; static MPrism *setHighOrder(MPrism *p, GRegion *gr,
for(unsigned int i = 0; i < gr->prisms.size(); i++){ edgeContainer &edgeVertices,
MPrism *p = gr->prisms[i]; faceContainer &faceVertices,
bool linear, bool incomplete, int nPts)
{
std::vector<MVertex*> ve, vf, vr; std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts); getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
if(incomplete){ if(incomplete){
prisms2.push_back if(nPts == 1){
(new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2), return new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5), 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], ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
0, p->getPartition())); 0, p->getPartition());
}
else{
Msg::Error("PrismN generation not implemented");
return 0;
}
} }
else{ else{
if (nPts == 1) { if (nPts == 1) {
getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts); getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
prisms2.push_back return new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2),
(new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2),
p->getVertex(3), p->getVertex(4), p->getVertex(5), 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], ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
vf[0], vf[1], vf[2], vf[0], vf[1], vf[2],
0, p->getPartition())); 0, p->getPartition());
} }
else { else {
Msg::Error("PrismN generation not implemented"); Msg::Error("PrismN generation not implemented");
/* return 0;
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);
*/
} }
} }
delete p;
}
gr->prisms = prisms2;
std::vector<MPyramid*> pyramids2; static MPyramid *setHighOrder(MPyramid *p, GRegion *gr,
for(unsigned int i = 0; i < gr->pyramids.size(); i++) { edgeContainer &edgeVertices,
MPyramid *p = gr->pyramids[i]; faceContainer &faceVertices,
bool linear, bool incomplete, int nPts)
{
std::vector<MVertex*> ve, vf, vr; std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts); getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts); getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
...@@ -1003,16 +997,56 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices, ...@@ -1003,16 +997,56 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
} }
} }
ve.insert(ve.end(), vr.begin(), vr.end()); ve.insert(ve.end(), vr.begin(), vr.end());
MPyramid *n = new MPyramidN(p->getVertex(0), p->getVertex(1), return new MPyramidN(p->getVertex(0), p->getVertex(1),
p->getVertex(2), p->getVertex(3), p->getVertex(2), p->getVertex(3),
p->getVertex(4), ve, nPts + 1, p->getVertex(4), ve, nPts + 1,
0, p->getPartition()); 0, p->getPartition());
pyramids2.push_back(n); }
SPoint3 test_pnt;
n->pnt(-1,-1,0, test_pnt); 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; delete p;
} }
gr->pyramids = pyramids2; gr->pyramids = pyramids2;
gr->deleteVertexArrays(); gr->deleteVertexArrays();
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment