Commit 1edd4845 authored by Christophe Geuzaine's avatar Christophe Geuzaine

Merge branch 'pyramids' into 'master'

new automatic hybrid mesh generation algorithm (tet mesh with quads on boundary,…

See merge request !216
parents 57d766c0 0c6433f0
Pipeline #2883 passed with stage
in 57 minutes and 41 seconds
This diff is collapsed.
......@@ -64,9 +64,22 @@ GFace *findInFaceSearchStructure(MVertex *p1, MVertex *p2, MVertex *p3,
GFace *findInFaceSearchStructure(const MFace &f, const fs_cont &search);
GEdge *findInEdgeSearchStructure(MVertex *p1, MVertex *p2,
const es_cont &search);
bool buildFaceSearchStructure(GModel *model, fs_cont &search);
bool buildFaceSearchStructure(GModel *model, fs_cont &search, bool onlyTriangles = false);
bool buildEdgeSearchStructure(GModel *model, es_cont &search);
// hybrid mesh recovery structure
class splitQuadRecovery {
private:
std::map<MFace, MVertex *, Less_Face> _quad;
std::map<MFace, GFace *, Less_Face> _tri;
public:
splitQuadRecovery() {}
void add(const MFace &f, MVertex *v, GFace *gf);
std::map<MFace, GFace *, Less_Face> &getTri() { return _tri; }
std::map<MFace, MVertex *, Less_Face> &getQuad() { return _quad; }
int buildPyramids(GModel *gm);
};
// adapt the mesh of a region
class adaptMeshGRegion {
public:
......
......@@ -12,6 +12,7 @@
#if defined(HAVE_TETGENBR)
#include "meshGRegion.h"
#include "meshGRegionDelaunayInsertion.h"
#include "robustPredicates.h"
#include "GModel.h"
......@@ -21,6 +22,7 @@
#include "MLine.h"
#include "MPoint.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "Context.h"
#include "OS.h"
......@@ -43,6 +45,11 @@ namespace tetgenBR {
#define REAL double
struct brdata{
GRegion *gr;
splitQuadRecovery *sqr;
};
// dummy tetgenio class
class tetgenio {
public:
......@@ -115,7 +122,8 @@ namespace tetgenBR {
bool tetgenmesh::reconstructmesh(void *p)
{
GRegion *_gr = (GRegion *)p;
GRegion *_gr = ((brdata *)p)->gr;
splitQuadRecovery *_sqr = ((brdata *)p)->sqr;
in = new tetgenio();
b = new tetgenbehavior();
......@@ -135,9 +143,33 @@ namespace tetgenBR {
++it) {
GFace *gf = *it;
for(unsigned int i = 0; i < gf->triangles.size(); i++) {
all.insert(gf->triangles[i]->getVertex(0));
all.insert(gf->triangles[i]->getVertex(1));
all.insert(gf->triangles[i]->getVertex(2));
MVertex *v0 = gf->triangles[i]->getVertex(0);
MVertex *v1 = gf->triangles[i]->getVertex(1);
MVertex *v2 = gf->triangles[i]->getVertex(2);
all.insert(v0);
all.insert(v1);
all.insert(v2);
}
if(_sqr){
for(unsigned int i = 0; i < gf->quadrangles.size(); i++) {
MVertex *v0 = gf->quadrangles[i]->getVertex(0);
MVertex *v1 = gf->quadrangles[i]->getVertex(1);
MVertex *v2 = gf->quadrangles[i]->getVertex(2);
MVertex *v3 = gf->quadrangles[i]->getVertex(3);
MVertex *newv = new MVertex
((v0->x() + v1->x() + v2->x() + v3->x()) * 0.25,
(v0->y() + v1->y() + v2->y() + v3->y()) * 0.25,
(v0->z() + v1->z() + v2->z() + v3->z()) * 0.25, gf);
// the extra vertex will be added in a GRegion (and reclassified
// correctly on that GRegion) when the pyramid is generated
MFace mf = gf->quadrangles[i]->getFace(0);
_sqr->add(mf, newv, gf);
all.insert(v0);
all.insert(v1);
all.insert(v2);
all.insert(v3);
all.insert(newv);
}
}
}
std::vector<GEdge *> const &e = _gr->embeddedEdges();
......@@ -425,7 +457,6 @@ namespace tetgenBR {
setpointtype(p[j], FACETVERTEX);
}
}
// Create an initial triangulation.
makeshellface(subfaces, &newsh);
setshvertices(newsh, p[0], p[1], p[2]);
setshellmark(newsh, gf->tag()); // the GFace's tag.
......@@ -438,9 +469,35 @@ namespace tetgenBR {
ssbond(newsh, newseg);
senextself(newsh);
}
} // i
}
} // it
if(_sqr){
std::map<MFace, GFace *, Less_Face> f = _sqr->getTri();
for(std::map<MFace, GFace *, Less_Face>::iterator it = f.begin();
it != f.end(); it++){
const MFace &mf = it->first;
for(int j = 0; j < 3; j++) {
p[j] = idx2verlist[mf.getVertex(j)->getIndex()];
if(pointtype(p[j]) == VOLVERTEX) {
setpointtype(p[j], FACETVERTEX);
}
}
makeshellface(subfaces, &newsh);
setshvertices(newsh, p[0], p[1], p[2]);
setshellmark(newsh, it->second->tag());
recentsh = newsh;
for(int j = 0; j < 3; j++) {
makeshellface(subsegs, &newseg);
setshvertices(newseg, sorg(newsh), sdest(newsh), NULL);
// Set the default segment marker '-1'.
setshellmark(newseg, -1);
ssbond(newsh, newseg);
senextself(newsh);
}
}
}
// Connecting triangles, removing redundant segments.
unifysegments();
......@@ -843,6 +900,11 @@ namespace tetgenBR {
delete gf->triangles[i];
gf->triangles.clear();
gf->deleteVertexArrays();
if(gf->quadrangles.size()){
Msg::Warning("Steiner points not handled for quad surface mesh");
}
// Create the new triangles.
subloop.shver = 0;
subfaces->traversalinit();
......@@ -1117,12 +1179,13 @@ namespace tetgenBR {
} // end namespace
bool meshGRegionBoundaryRecovery(GRegion *gr)
bool meshGRegionBoundaryRecovery(GRegion *gr, splitQuadRecovery *sqr)
{
bool ret = false;
try {
tetgenBR::tetgenmesh *m = new tetgenBR::tetgenmesh();
ret = m->reconstructmesh((void *)gr);
tetgenBR::brdata data = {gr, sqr};
ret = m->reconstructmesh((void *)&data);
delete m;
} catch(int err) {
if(err == 1) {
......@@ -1253,6 +1316,9 @@ bool meshGRegionBoundaryRecovery(GRegion *gr)
#else
bool meshGRegionBoundaryRecovery(GRegion *gr) { return false; }
bool meshGRegionBoundaryRecovery(GRegion *gr, splitQuadRecovery *sqr)
{
return false;
}
#endif
......@@ -7,7 +7,8 @@
#define _MESH_GREGION_BOUNDARY_RECOVERY_H_
class GRegion;
class splitQuadRecovery;
bool meshGRegionBoundaryRecovery(GRegion *gr);
bool meshGRegionBoundaryRecovery(GRegion *gr, splitQuadRecovery *sqr=0);
#endif
......@@ -663,7 +663,7 @@ void non_recursive_classify(MTet4 *t, std::list<MTet4 *> &theRegion,
t = _stackounette.top();
_stackounette.pop();
if(!t) {
Msg::Error("a tet is not connected by a boundary face");
Msg::Error("A tetrahedron is not connected to a boundary face");
touchesOutsideBox = true;
}
else if(!t->onWhat()) {
......@@ -1246,7 +1246,8 @@ int isCavityCompatibleWithEmbeddedFace(
return 1;
}
void insertVerticesInRegion(GRegion *gr, int maxVert, bool _classify)
void insertVerticesInRegion(GRegion *gr, int maxVert, bool _classify,
splitQuadRecovery *sqr)
{
// TEST_IF_BOUNDARY_IS_RECOVERED (gr);
......@@ -1335,7 +1336,10 @@ void insertVerticesInRegion(GRegion *gr, int maxVert, bool _classify)
if(_classify) {
fs_cont search;
buildFaceSearchStructure(gr->model(), search);
buildFaceSearchStructure(gr->model(), search, true); // only triangles
if(sqr)
search.insert(sqr->getTri().begin(), sqr->getTri().end());
for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end();
++it) {
if(!(*it)->onWhat()) {
......@@ -1362,8 +1366,11 @@ void insertVerticesInRegion(GRegion *gr, int maxVert, bool _classify)
// Make sure that Steiner points will end up in the right region
std::vector<MVertex *> vertices;
(*it2)->tet()->getVertices(vertices);
for(std::vector<MVertex *>::iterator itv = vertices.begin(); itv != vertices.end(); ++itv){
if ((*itv)->onWhat() != NULL && (*itv)->onWhat()->dim() == 3 && (*itv)->onWhat() != myGRegion){
for(std::vector<MVertex *>::iterator itv = vertices.begin();
itv != vertices.end(); ++itv){
if ((*itv)->onWhat() != NULL &&
(*itv)->onWhat()->dim() == 3 &&
(*itv)->onWhat() != myGRegion){
myGRegion->addMeshVertex((*itv));
(*itv)->setEntity(myGRegion);
}
......
......@@ -20,6 +20,7 @@
class GRegion;
class GFace;
class GModel;
class splitQuadRecovery;
double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
double circumcenter[3], double *xi, double *eta,
......@@ -219,7 +220,7 @@ void connectTets(std::vector<MTet4 *> &, const std::set<MFace, Less_Face> * = 0)
void delaunayMeshIn3D(std::vector<MVertex *> &, std::vector<MTetrahedron *> &,
bool removeBox = true);
void insertVerticesInRegion(GRegion *gr, int maxVert = 2000000000,
bool _classify = true);
bool _classify = true, splitQuadRecovery *sqr = 0);
void bowyerWatsonFrontalLayers(GRegion *gr, bool hex);
struct compareTet4Ptr {
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment