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

rewrote gmsh2DMeshGenerator to make it thread-safe

*** PLEASE TEST THAT THE 2D MESH GENERATORS BEHAVE EXACTLY AS BEFORE ***
parent be6c7619
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Mesh/meshGFace.cpp
+ 318
− 476
View file @ 82ad0b39
......@@ -31,59 +31,6 @@
#include "OS.h"
#include "HighOrder.h"
void computeEdgeLoops(const GFace *gf,
std::vector<MVertex*> &all_mvertices,
std::vector<int> &indices)
{
std::list<GEdge*> edges = gf->edges();
std::list<int> ori = gf->orientations();
std::list<GEdge*>::iterator it = edges.begin();
std::list<int>::iterator ito = ori.begin();
indices.push_back(0);
GVertex *start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
GVertex *v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
all_mvertices.push_back(start->mesh_vertices[0]);
if (*ito == 1)
for (unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
all_mvertices.push_back((*it)->mesh_vertices[i]);
else
for (int i = (*it)->mesh_vertices.size() - 1; i >= 0; i--)
all_mvertices.push_back((*it)->mesh_vertices[i]);
GVertex *v_start = start;
while(1){
++it;
++ito;
if(v_end == start){
indices.push_back(all_mvertices.size());
if(it == edges.end ()) break;
start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
v_start = start;
}
else{
if(it == edges.end ()){
Msg::Error("Something wrong in edge loop computation");
return;
}
v_start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
if(v_start != v_end){
Msg::Error("Something wrong in edge loop computation");
return;
}
v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
}
all_mvertices.push_back(v_start->mesh_vertices[0]);
if(*ito == 1)
for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
all_mvertices.push_back((*it)->mesh_vertices[i]);
else
for (int i = (*it)->mesh_vertices.size()-1; i >= 0; i--)
all_mvertices.push_back((*it)->mesh_vertices[i]);
}
}
void fourthPoint(double *p1, double *p2, double *p3, double *p4)
{
double c[3];
......@@ -113,7 +60,8 @@ static bool noseam(GFace *gf)
return true;
}
static void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv,
std::set<EdgeToRecover> &edgesNotRecovered,
std::list<GFace*> &facesToRemesh)
{
facesToRemesh.clear();
......@@ -122,7 +70,7 @@ static void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
for(; itr != edgesNotRecovered.end(); ++itr){
std::list<GFace*> l_faces = itr->ge->faces();
// Un-mesh model faces adjacent to the model edge
// un-mesh model faces adjacent to the model edge
for(std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end(); ++it){
if((*it)->triangles.size() || (*it)->quadrangles.size()){
facesToRemesh.push_back(*it);
......@@ -143,8 +91,9 @@ static void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
for(int i = 0; i < N; i++){
MVertex *v1 = itr->ge->lines[i]->getVertex(0);
MVertex *v2 = itr->ge->lines[i]->getVertex(1);
if ((v1->getIndex() == p1 && v2->getIndex() == p2) ||
(v1->getIndex() == p2 && v2->getIndex() == p1)){
BDS_Point *pp1 = recoverMapInv[v1];
BDS_Point *pp2 = recoverMapInv[v2];
if((pp1->iD == p1 && pp2->iD == p2) || (pp1->iD == p2 && pp2->iD == p1)){
double t1;
double lc1 = -1;
if(v1->onWhat() == g1) t1 = bb.low();
......@@ -223,7 +172,9 @@ void computeElementShapes(GFace *gf, double &worst, double &avg, double &best,
avg /= nT;
}
static bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
static bool recover_medge(BDS_Mesh *m, GEdge *ge,
std::map<MVertex*, BDS_Point*> &recoverMapInv,
std::set<EdgeToRecover> *e2r,
std::set<EdgeToRecover> *not_recovered, int pass_)
{
BDS_GeomEntity *g = 0;
......@@ -235,32 +186,35 @@ static bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
for(unsigned int i = 0; i < ge->lines.size(); i++){
MVertex *vstart = ge->lines[i]->getVertex(0);
MVertex *vend = ge->lines[i]->getVertex(1);
if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
BDS_Point *pstart = recoverMapInv[vstart];
BDS_Point *pend = recoverMapInv[vend];
if(pass_ == 1)
e2r->insert(EdgeToRecover(pstart->iD, pend->iD, ge));
else{
BDS_Edge *e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, e2r, not_recovered);
if(e) e->g = g;
else {
// else {
// Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
// vstart->x(), vstart->y(), vend->x(), vend->y(), i,
// ge->mesh_vertices.size());
// return false;
}
// }
}
}
if(pass_ == 2 && ge->getBeginVertex()){
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
BDS_Point *pstart = m->find_point(vstart->getIndex());
BDS_Point *pend = m->find_point(vend->getIndex());
BDS_Point *pstart = recoverMapInv[vstart];
BDS_Point *pend = recoverMapInv[vend];
if(!pstart->g){
m->add_geom (vstart->getIndex(), 0);
BDS_GeomEntity *g0 = m->get_geom(vstart->getIndex(), 0);
m->add_geom(pstart->iD, 0);
BDS_GeomEntity *g0 = m->get_geom(pstart->iD, 0);
pstart->g = g0;
}
if(!pend->g){
m->add_geom (vend->getIndex(), 0);
BDS_GeomEntity *g0 = m->get_geom(vend->getIndex(), 0);
m->add_geom(pend->iD, 0);
BDS_GeomEntity *g0 = m->get_geom(pend->iD, 0);
pend->g = g0;
}
}
......@@ -268,147 +222,37 @@ static bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
return true;
}
static bool recover_medge_old(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
std::set<EdgeToRecover> *not_recovered, int pass_)
{
BDS_GeomEntity *g = 0;
if (pass_ == 2){
m->add_geom (ge->tag(), 1);
g = m->get_geom(ge->tag(),1);
}
if(ge->mesh_vertices.size() == 0){
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
if(pass_ == 1){
e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
return true;
}
else{
BDS_Point *pstart = m->find_point(vstart->getIndex());
BDS_Point *pend = m->find_point(vend->getIndex());
if(!pstart->g){
m->add_geom (vstart->getIndex(), 0);
BDS_GeomEntity *g0 = m->get_geom(vstart->getIndex(), 0);
pstart->g = g0;
}
if(!pend->g){
m->add_geom(vend->getIndex(), 0);
BDS_GeomEntity *g0 = m->get_geom(vend->getIndex(), 0);
pend->g = g0;
}
BDS_Edge * e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e) e->g = g;
else {
// Msg::Error("The unrecoverable edge is on model edge %d", ge->tag());
return false;
}
return true;
}
}
BDS_Edge *e;
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->mesh_vertices.begin());
if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
else{
BDS_Point *pstart = m->find_point(vstart->getIndex());
if(!pstart->g){
m->add_geom (vstart->getIndex(), 0);
BDS_GeomEntity *g0 = m->get_geom(vstart->getIndex(), 0);
pstart->g = g0;
}
e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e) e->g = g;
else {
// Msg::Error("The unrecoverable edge is on model edge %d", ge->tag());
// return false;
}
}
for (unsigned int i = 1; i < ge->mesh_vertices.size(); i++){
vstart = ge->mesh_vertices[i - 1];
vend = ge->mesh_vertices[i];
if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
else{
e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e) e->g = g;
else {
// Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
// vstart->x(), vstart->y(), vend->x(), vend->y(), i,
// ge->mesh_vertices.size());
// return false;
}
}
}
vstart = vend;
vend = *(ge->getEndVertex()->mesh_vertices.begin());
if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
else{
e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
if (e)e->g = g;
else {
// Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
// vstart->x(), vstart->y(), vend->x(), vend->y(),
// ge->mesh_vertices.size(), ge->mesh_vertices.size());
// return false;
}
BDS_Point *pend = m->find_point(vend->getIndex());
if(!pend->g){
m->add_geom (vend->getIndex(), 0);
BDS_GeomEntity *g0 = m->get_geom(vend->getIndex(), 0);
pend->g = g0;
}
}
return true;
}
// Builds An initial triangular mesh that respects the boundaries of
// the domain, including embedded points and surfaces
// FIXME: to make the algorithm thread-safe, we need to change the way
// we renumber the boundary vertices (and the associated u/v reparam):
// using setIndex leads to unpredicatbel results if meshing
// concurrently two surfaces that share a common edge!
static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
bool repairSelfIntersecting1dMesh, bool debug = true)
bool repairSelfIntersecting1dMesh,
bool debug = true)
{
BDS_GeomEntity CLASS_F(1, 2);
typedef std::set<MVertex*> v_container;
v_container all_vertices;
std::map<int, MVertex*>numbered_vertices;
std::map<BDS_Point*, MVertex*> recoverMap;
std::map<MVertex*, BDS_Point*> recoverMapInv;
std::list<GEdge*> edges = gf->edges();
// here, we will replace edges by their compounds
// replace edges by their compounds
if(gf->geomType() == GEntity::CompoundSurface){
//printf("replace edges by compound lines \n");
std::set<GEdge*> mySet;
std::list<GEdge*>::iterator it = edges.begin();
while(it != edges.end()){
if ((*it)->getCompound()){
if((*it)->getCompound())
mySet.insert((*it)->getCompound());
//printf("compound edge %d found in edge %d\n",(*it)->getCompound()->tag(), (*it)->tag());
}
else
mySet.insert(*it);
++it;
}
//printf("replacing %d edges by %d in the GFaceCompound %d\n",edges.size(),mySet.size(),gf->tag());
edges.clear();
edges.insert(edges.begin(), mySet.begin(), mySet.end());
}
// build a set with all points of the boundaries
std::set<MVertex*> all_vertices;
std::list<GEdge*> emb_edges = gf->embeddedEdges();
std::list<GVertex*> emb_vertx = gf->embeddedVertices();
std::list<GEdge*>::iterator it = edges.begin();
std::list<GVertex*>::iterator itvx = emb_vertx.begin();
// build a set with all points of the boundaries
it = edges.begin();
while(it != edges.end()){
if((*it)->isSeam(gf)) return false;
if(!(*it)->isMeshDegenerated()){
......@@ -420,6 +264,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
++it;
}
std::list<GEdge*> emb_edges = gf->embeddedEdges();
it = emb_edges.begin();
while(it != emb_edges.end()){
if(!(*it)->isMeshDegenerated()){
......@@ -434,6 +279,8 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
}
// add embedded vertices
std::list<GVertex*> emb_vertx = gf->embeddedVertices();
std::list<GVertex*>::iterator itvx = emb_vertx.begin();
while(itvx != emb_vertx.end()){
all_vertices.insert((*itvx)->mesh_vertices.begin(),
(*itvx)->mesh_vertices.end());
......@@ -448,111 +295,90 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
return true;
}
double *U_ = new double[all_vertices.size()];
double *V_ = new double[all_vertices.size()];
v_container::iterator itv = all_vertices.begin();
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
m->scalingU = 1;
m->scalingV = 1;
int count = 0;
std::vector<BDS_Point*> points(all_vertices.size());
std::vector<MVertex*> indexed_vertices(all_vertices.size());
SBoundingBox3d bbox;
while(itv != all_vertices.end()){
MVertex *here = *itv;
int count = 0;
for(std::set<MVertex*>::iterator it = all_vertices.begin();
it != all_vertices.end(); it++){
MVertex *here = *it;
GEntity *ge = here->onWhat();
SPoint2 param;
reparamMeshVertexOnFace(here, gf, param);
U_[count] = param[0];
V_[count] = param[1];
//printf("-->> meshGFace : %g %g -> u,v = %g %g\n",here->x(),here->y(), param.x(),param.y() );
(*itv)->setIndex(count);
numbered_vertices[(*itv)->getIndex()] = *itv;
bbox += SPoint3(param.x(), param.y(), 0);
BDS_Point *pp = m->add_point(count, param[0], param[1], gf);
m->add_geom(ge->tag(), ge->dim());
BDS_GeomEntity *g = m->get_geom(ge->tag(), ge->dim());
pp->g = g;
recoverMap[pp] = here;
recoverMapInv[here] = pp;
points[count] = pp;
bbox += SPoint3(param[0], param[1], 0);
count++;
++itv;
}
all_vertices.clear();
// compute the bounding box in parametric space. I do not have
// SBoundingBox, so I use a 3D one... At the same time, number the
// vertices locally.
// compute the bounding box in parametric space
SVector3 dd(bbox.max(), bbox.min());
double LC2D = norm(dd);
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
m->scalingU = 1;
m->scalingV = 1;
// Use a divide & conquer type algorithm to create a triangulation.
// use a divide & conquer type algorithm to create a triangulation.
// We add to the triangulation a box with 4 points that encloses the
// domain.
{
DocRecord doc(all_vertices.size() + 4);
itv = all_vertices.begin();
int j = 0;
while(itv != all_vertices.end()){
MVertex *here = *itv;
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
doc.points[j].where.h = U_[j] + XX;
doc.points[j].where.v = V_[j] + YY;
doc.points[j].adjacent = NULL;
doc.points[j].data = here;
j++;
++itv;
}
// Increase the size of the bounding box
DocRecord doc(points.size() + 4);
for(unsigned int i = 0; i < points.size(); i++){
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
(double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
(double)RAND_MAX;
doc.points[i].where.h = points[i]->u + XX;
doc.points[i].where.v = points[i]->v + YY;
doc.points[i].data = points[i];
doc.points[i].adjacent = NULL;
}
// increase the size of the bounding box
bbox *= 2.5;
// add 4 points than encloses the domain (use negative number to
// distinguish thos fake vertices)
MVertex *bb[4];
bb[0] = new MVertex(bbox.min().x(), bbox.min().y(), 0, gf, -1);
bb[1] = new MVertex(bbox.min().x(), bbox.max().y(), 0, gf, -2);
bb[2] = new MVertex(bbox.max().x(), bbox.min().y(), 0, gf, -3);
bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, gf, -4);
// add 4 points than encloses the domain (use negative number to
// distinguish those fake vertices)
double bb[4][2] = {{bbox.min().x(), bbox.min().y()},
{bbox.min().x(), bbox.max().y()},
{bbox.max().x(), bbox.min().y()},
{bbox.max().x(), bbox.max().y()}};
for(int ip = 0; ip < 4; ip++){
doc.points[all_vertices.size() + ip].where.h = bb[ip]->x();
doc.points[all_vertices.size() + ip].where.v = bb[ip]->y();
doc.points[all_vertices.size() + ip].adjacent = 0;
doc.points[all_vertices.size() + ip].data = bb[ip];
BDS_Point *pp = m->add_point(-ip - 1, bb[ip][0], bb[ip][1], gf);
m->add_geom(gf->tag(), 2);
BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
pp->g = g;
doc.points[points.size() + ip].where.h = bb[ip][0];
doc.points[points.size() + ip].where.v = bb[ip][1];
doc.points[points.size() + ip].adjacent = 0;
doc.points[points.size() + ip].data = pp;
}
// Use "fast" inhouse recursive algo to generate the triangulation
// Use "fast" inhouse recursive algo to generate the triangulation.
// At this stage the triangulation is not what we need
// -) It does not necessary recover the boundaries
// -) It contains triangles outside the domain (the first edge
// loop is the outer one)
Msg::Debug("Meshing of the convex hull (%d points)", all_vertices.size());
Msg::Debug("Meshing of the convex hull (%d points)", points.size());
doc.MakeMeshWithPoints();
Msg::Debug("Meshing of the convex hull (%d points) done", all_vertices.size());
for(int i = 0; i < doc.numPoints; i++){
MVertex *here = (MVertex *)doc.points[i].data;
int num = here->getIndex();
double U, V;
if(num < 0){ // fake bbox points
U = bb[-1 - num]->x();
V = bb[-1 - num]->y();
}
else{
U = U_[num];
V = V_[num];
}
m->add_point(num, U, V, gf);
}
Msg::Debug("Meshing of the convex hull (%d points) done", points.size());
for(int i = 0; i < doc.numTriangles; i++) {
MVertex *V1 = (MVertex*)doc.points[doc.triangles[i].a].data;
MVertex *V2 = (MVertex*)doc.points[doc.triangles[i].b].data;
MVertex *V3 = (MVertex*)doc.points[doc.triangles[i].c].data;
m->add_triangle(V1->getIndex(), V2->getIndex(), V3->getIndex());
BDS_Point *p1 = (BDS_Point*)doc.points[doc.triangles[i].a].data;
BDS_Point *p2 = (BDS_Point*)doc.points[doc.triangles[i].b].data;
BDS_Point *p3 = (BDS_Point*)doc.points[doc.triangles[i].c].data;
m->add_triangle(p1->iD, p2->iD, p3->iD);
}
// Recover the boundary edges and compute characteristic lenghts
// using mesh edge spacing
if(debug && RECUR_ITER == 0){
char name[245];
sprintf(name, "surface%d-initial-real.pos", gf->tag());
......@@ -561,23 +387,24 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
outputScalarField(m->triangles, name, 1);
}
// Recover the boundary edges and compute characteristic lenghts
// using mesh edge spacing. If two of these edges intersect, then
// the 1D mesh have to be densified
Msg::Debug("Recovering %d model Edges", edges.size());
// build a structure with all mesh edges that have to be
// recovered. If two of these edges intersect, then the 1D mesh have
// to be densified
std::set<EdgeToRecover> edgesToRecover;
std::set<EdgeToRecover> edgesNotRecovered;
it = edges.begin();
while(it != edges.end()){
if(!(*it)->isMeshDegenerated())
recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 1);
recover_medge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
++it;
}
it = emb_edges.begin();
while(it != emb_edges.end()){
if(!(*it)->isMeshDegenerated())
recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 1);
recover_medge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
++it;
}
......@@ -587,7 +414,8 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
it = edges.begin();
while(it != edges.end()){
if(!(*it)->isMeshDegenerated()){
recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 2);
recover_medge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
}
++it;
}
......@@ -603,13 +431,14 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
if(debug){
char name[245];
sprintf(name, "surface%d-not_yet_recovered-real-%d.msh", gf->tag(),RECUR_ITER);
sprintf(name, "surface%d-not_yet_recovered-real-%d.msh", gf->tag(),
RECUR_ITER);
gf->model()->writeMSH(name);
}
std::list<GFace *> facesToRemesh;
if(repairSelfIntersecting1dMesh)
remeshUnrecoveredEdges(edgesNotRecovered, facesToRemesh);
remeshUnrecoveredEdges(recoverMapInv, edgesNotRecovered, facesToRemesh);
else{
std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
int *_error = new int[3 * edgesNotRecovered.size()];
......@@ -628,20 +457,20 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// delete the mesh
delete m;
delete [] U_;
delete [] V_;
if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
return gmsh2DMeshGenerator(gf, RECUR_ITER+1, repairSelfIntersecting1dMesh, debug);
return gmsh2DMeshGenerator
(gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, debug);
return false;
}
if(RECUR_ITER > 0)
Msg::Warning(":-) Gmsh was able to recover all edges after %d iterations", RECUR_ITER);
Msg::Warning(":-) Gmsh was able to recover all edges after %d iterations",
RECUR_ITER);
Msg::Debug("Boundary Edges recovered for surface %d", gf->tag());
// Look for an edge that is on the boundary for which one of the two
// neighbors has a negative number node. The other triangle is
// Look for an edge that is on the boundary for which one of the
// two neighbors has a negative number node. The other triangle is
// inside the domain and, because all edges were recovered,
// triangles inside the domain can be recovered using a simple
// recursive algorithm
......@@ -668,16 +497,19 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
it = emb_edges.begin();
while(it != emb_edges.end()){
if(!(*it)->isMeshDegenerated())
recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 2);
recover_medge
(m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
++it;
}
// compute characteristic lengths at vertices
Msg::Debug("Computing mesh size field at mesh vertices %d", edgesToRecover.size());
Msg::Debug("Computing mesh size field at mesh vertices %d",
edgesToRecover.size());
for(int i = 0; i < doc.numPoints; i++){
MVertex *here = (MVertex *)doc.points[i].data;
int num = here->getIndex();
BDS_Point *pp = (BDS_Point*)doc.points[i].data;
if(recoverMap.count(pp)){
MVertex *here = recoverMap[pp];
GEntity *ge = here->onWhat();
BDS_Point *pp = m->find_point(num);
if(ge->dim() == 0){
pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
}
......@@ -690,15 +522,14 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
pp->lcBGM() = MAX_LC;
pp->lc() = pp->lcBGM();
}
for(int ip = 0; ip < 4; ip++) delete bb[ip];
}
}
// delete useless stuff
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end()){
BDS_Face *t = *itt;
if (!t->g)
m->del_face (t);
if(!t->g) m->del_face(t);
++itt;
}
m->cleanup();
......@@ -712,8 +543,10 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
else{
if(!e->g)
e->g = &CLASS_F;
if (!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree)e->p1->g = e->g;
if (!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree)e->p2->g = e->g;
if(!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree)
e->p1->g = e->g;
if(!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree)
e->p2->g = e->g;
}
++ite;
}
......@@ -759,9 +592,10 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
while (itp != m->points.end()){
BDS_Point *p = *itp;
if (numbered_vertices.find(p->iD) == numbered_vertices.end()){
MVertex *v = new MFaceVertex (p->X, p->Y, p->Z, gf, p->u, p->v);
numbered_vertices[p->iD]=v;
if(recoverMap.find(p) == recoverMap.end()){
MVertex *v = new MFaceVertex
(p->X, p->Y, p->Z, gf, m->scalingU * p->u, m->scalingV * p->v);
recoverMap[p] = v;
gf->mesh_vertices.push_back(v);
}
++itp;
......@@ -774,13 +608,18 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
if(!t->deleted){
BDS_Point *n[4];
t->getNodes(n);
MVertex *v1 = numbered_vertices[n[0]->iD];
MVertex *v2 = numbered_vertices[n[1]->iD];
MVertex *v3 = numbered_vertices[n[2]->iD];
if (!n[3])
MVertex *v1 = recoverMap[n[0]];
MVertex *v2 = recoverMap[n[1]];
MVertex *v3 = recoverMap[n[2]];
if(!n[3]){
// when a singular point is present, degenerated triangles
// may be created, for example on a sphere that contains one
// pole
if(v1 != v2 && v1 != v3 && v2 != v3)
gf->triangles.push_back(new MTriangle(v1, v2, v3));
}
else{
MVertex *v4 = numbered_vertices[n[3]->iD];
MVertex *v4 = recoverMap[n[3]];
gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
}
}
......@@ -809,8 +648,6 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// delete the mesh
delete m;
delete [] U_;
delete [] V_;
if(gf->meshAttributes.recombine)
gmshRecombineIntoQuads(gf);
......@@ -846,7 +683,7 @@ static void printMesh1d(int iEdge, int seam, std::vector<SPoint2> &m)
static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
std::vector<BDS_Point*> &result,
SBoundingBox3d &bbox, BDS_Mesh *m,
std::map<BDS_Point*, MVertex*> &recover_map_global,
std::map<BDS_Point*, MVertex*> &recoverMap,
int &count, int countTot, double tol,
bool seam_the_first = false)
{
......@@ -859,7 +696,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
const int MYDEBUG = false;
std::map<BDS_Point*,MVertex*> recover_map;
std::map<BDS_Point*, MVertex*> recoverMapLocal;
result.clear();
count = 0;
......@@ -1000,8 +837,8 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
// It has not worked : either tolerance is wrong or the first seam edge
// has to be taken with the other parametric coordinates (because it is
// only present once in the closure of the domain).
for (std::map<BDS_Point*, MVertex*>::iterator it = recover_map.begin();
it != recover_map.end(); ++it){
for(std::map<BDS_Point*, MVertex*>::iterator it = recoverMapLocal.begin();
it != recoverMapLocal.end(); ++it){
m->del_point(it->first);
}
return false;
......@@ -1053,7 +890,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
pp->g->classif_degree);
bbox += SPoint3(U, V, 0);
edgeLoop_BDS.push_back(pp);
recover_map[pp] = here;
recoverMapLocal[pp] = here;
count++;
}
last_coord = coords[coords.size() - 1];
......@@ -1062,14 +899,14 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
}
// It has worked, so we add all the points to the recover map
recover_map_global.insert(recover_map.begin(),recover_map.end());
recoverMap.insert(recoverMapLocal.begin(), recoverMapLocal.end());
return true;
}
static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
{
std::map<BDS_Point*,MVertex*> recover_map;
std::map<BDS_Point*, MVertex*> recoverMap;
Range<double> rangeU = gf->parBounds(0);
Range<double> rangeV = gf->parBounds(1);
......@@ -1082,8 +919,9 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
m->scalingU = 1;//fabs(du);
m->scalingV = 1;//fabs(dv);
m->scalingU = 1;
m->scalingV = 1;
std::vector<std::vector<BDS_Point*> > edgeLoops_BDS;
SBoundingBox3d bbox;
int nbPointsTotal = 0;
......@@ -1096,13 +934,13 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
bool ok = false;
for(int i = 0; i < 4; i++){
if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
recover_map, nbPointsLocal,
recoverMap, nbPointsLocal,
nbPointsTotal, fact[i] * LC2D)){
ok = true;
break;
}
if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
recover_map, nbPointsLocal,
recoverMap, nbPointsLocal,
nbPointsTotal, fact[i] * LC2D, true)){
ok = true;
break;
......@@ -1129,14 +967,12 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
BDS_Point *pp = edgeLoop_BDS[j];
const double U = pp->u;
const double V = pp->v;
double XX = CTX::instance()->mesh.randFactor * LC2D *
(double)rand() / (double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D *
(double)rand() / (double)RAND_MAX;
doc.points[count].where.h = U + XX;
doc.points[count].where.v = V + YY;
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
(double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
(double)RAND_MAX;
doc.points[count].where.h = pp->u + XX;
doc.points[count].where.v = pp->v + YY;
doc.points[count].adjacent = NULL;
doc.points[count].data = pp;
count++;
......@@ -1196,11 +1032,11 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
BDS_Edge * e = m->recover_edge(edgeLoop_BDS[j]->iD,
edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
BDS_Edge * e = m->recover_edge
(edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
if(!e){
Msg::Error("Impossible to recover the edge %d %d",
edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
Msg::Error("Impossible to recover the edge %d %d", edgeLoop_BDS[j]->iD,
edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
gf->meshStatistics.status = GFace::FAILED;
return false;
}
......@@ -1257,8 +1093,10 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
else{
if(!e->g)
e->g = &CLASS_F;
if (!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree)e->p1->g = e->g;
if (!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree)e->p2->g = e->g;
if(!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree)
e->p1->g = e->g;
if(!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree)
e->p2->g = e->g;
}
++ite;
}
......@@ -1282,7 +1120,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
gmshRefineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true);
gmshOptimizeMeshBDS(gf, *m, 2);
gmshRefineMeshBDS (gf, *m, -CTX::instance()->mesh.refineSteps, false);
gmshOptimizeMeshBDS(gf, *m, 2, &recover_map);
gmshOptimizeMeshBDS(gf, *m, 2, &recoverMap);
// compute mesh statistics
computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index,
gf->meshStatistics.longest_edge_length,
......@@ -1297,9 +1135,10 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
while (itp != m->points.end()){
BDS_Point *p = *itp;
if (recover_map.find(p) == recover_map.end()){
MVertex *v = new MFaceVertex (p->X,p->Y,p->Z,gf,m->scalingU*p->u,m->scalingV*p->v);
recover_map[p] = v;
if(recoverMap.find(p) == recoverMap.end()){
MVertex *v = new MFaceVertex
(p->X, p->Y, p->Z, gf, m->scalingU * p->u, m->scalingV * p->v);
recoverMap[p] = v;
gf->mesh_vertices.push_back(v);
}
++itp;
......@@ -1313,9 +1152,9 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
if(!t->deleted){
BDS_Point *n[4];
t->getNodes(n);
MVertex *v1 = recover_map[n[0]];
MVertex *v2 = recover_map[n[1]];
MVertex *v3 = recover_map[n[2]];
MVertex *v1 = recoverMap[n[0]];
MVertex *v2 = recoverMap[n[1]];
MVertex *v3 = recoverMap[n[2]];
if(!n[3]){
// when a singular point is present, degenerated triangles
// may be created, for example on a sphere that contains one
......@@ -1324,7 +1163,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
gf->triangles.push_back(new MTriangle(v1, v2, v3));
}
else{
MVertex *v4 = recover_map[n[3]];
MVertex *v4 = recoverMap[n[3]];
gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
}
}
......@@ -1375,12 +1214,10 @@ void deMeshGFace::operator() (GFace *gf)
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
}
const int debugSurface = -1;//-100;
//const int debugSurface = -100;
const int debugSurface = -1;
void meshGFace::operator() (GFace *gf)
{
gf->model()->setCurrentMeshEntity(gf);
if(debugSurface >= 0 && gf->tag() != debugSurface){
......@@ -1403,7 +1240,8 @@ void meshGFace::operator() (GFace *gf)
const char *algo = "Unknown";
if(AlgoDelaunay2D(gf))
algo = (CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL) ? "Frontal" : "Delaunay";
algo = (CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL) ?
"Frontal" : "Delaunay";
else if(CTX::instance()->mesh.algo2d == ALGO_2D_MESHADAPT_OLD)
algo = "MeshAdapt (old)";
else
......@@ -1417,12 +1255,14 @@ void meshGFace::operator() (GFace *gf)
Msg::Debug("Computing edge loops");
Msg::Debug("Generating the mesh");
if(noseam(gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
//gmsh2DMeshGeneratorPeriodic(gf, debugSurface >= 0 || debugSurface == -100);
gmsh2DMeshGenerator(gf, 0, repairSelfIntersecting1dMesh, debugSurface >= 0 || debugSurface == -100);
if(noseam(gf) || gf->getNativeType() == GEntity::GmshModel ||
gf->edgeLoops.empty()){
gmsh2DMeshGenerator(gf, 0, repairSelfIntersecting1dMesh,
debugSurface >= 0 || debugSurface == -100);
}
else{
if(!gmsh2DMeshGeneratorPeriodic(gf, debugSurface >= 0 || debugSurface == -100))
if(!gmsh2DMeshGeneratorPeriodic
(gf, debugSurface >= 0 || debugSurface == -100))
Msg::Error("Impossible to mesh face %d", gf->tag());
}
......@@ -1493,9 +1333,11 @@ void orientMeshGFace::operator()(GFace *gf)
GEdge *ge = *edges.begin();
GVertex *beg = ge->getBeginVertex();
GVertex *end = ge->getEndVertex();
if(!beg || beg->mesh_vertices.empty() || !end || end->mesh_vertices.empty()) return;
if(!beg || beg->mesh_vertices.empty() || !end || end->mesh_vertices.empty())
return;
MVertex *v1 = beg->mesh_vertices[0];
MVertex *v2 = ge->mesh_vertices.empty() ? end->mesh_vertices[0] : ge->mesh_vertices[0];
MVertex *v2 = ge->mesh_vertices.empty() ? end->mesh_vertices[0] :
ge->mesh_vertices[0];
int sign = *ori.begin();
MEdge ref(sign > 0 ? v1 : v2, sign > 0 ? v2 : v1);
if(shouldRevert(ref, gf->triangles) || shouldRevert(ref, gf->quadrangles)){
......
Mesh/meshGFace.h
+ 0
− 6
View file @ 82ad0b39
......@@ -38,12 +38,6 @@ class orientMeshGFace {
};
void fourthPoint(double *p1, double *p2, double *p3, double *p4);
// Compute edge loops of the face, all_mvertices are the vertices of
// the
void computeEdgeLoops(const GFace *gf,
std::vector<MVertex*> &all_mvertices,
std::vector<int> &indices);
void findTransfiniteCorners(GFace *gf, std::vector<MVertex*> &corners);
int MeshTransfiniteSurface(GFace *gf);
int MeshExtrudedSurface(GFace *gf, std::set<std::pair<MVertex*, MVertex*> >
......
Mesh/meshGFaceTransfinite.cpp
+ 52
− 0
View file @ 82ad0b39
......@@ -80,6 +80,58 @@ void findTransfiniteCorners(GFace *gf, std::vector<MVertex*> &corners)
}
}
static void computeEdgeLoops(const GFace *gf, std::vector<MVertex*> &all_mvertices,
std::vector<int> &indices)
{
std::list<GEdge*> edges = gf->edges();
std::list<int> ori = gf->orientations();
std::list<GEdge*>::iterator it = edges.begin();
std::list<int>::iterator ito = ori.begin();
indices.push_back(0);
GVertex *start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
GVertex *v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
all_mvertices.push_back(start->mesh_vertices[0]);
if(*ito == 1)
for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
all_mvertices.push_back((*it)->mesh_vertices[i]);
else
for(int i = (*it)->mesh_vertices.size() - 1; i >= 0; i--)
all_mvertices.push_back((*it)->mesh_vertices[i]);
GVertex *v_start = start;
while(1){
++it;
++ito;
if(v_end == start){
indices.push_back(all_mvertices.size());
if(it == edges.end ()) break;
start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
v_start = start;
}
else{
if(it == edges.end ()){
Msg::Error("Something wrong in edge loop computation");
return;
}
v_start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
if(v_start != v_end){
Msg::Error("Something wrong in edge loop computation");
return;
}
v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
}
all_mvertices.push_back(v_start->mesh_vertices[0]);
if(*ito == 1)
for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
all_mvertices.push_back((*it)->mesh_vertices[i]);
else
for(int i = (*it)->mesh_vertices.size()-1; i >= 0; i--)
all_mvertices.push_back((*it)->mesh_vertices[i]);
}
}
int MeshTransfiniteSurface(GFace *gf)
{
if(gf->meshAttributes.Method != MESH_TRANSFINITE) return 0;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment