From 82ad0b39b964741a743f03a4bd1c488ae46ac4d1 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 27 Aug 2009 21:42:29 +0000
Subject: [PATCH] rewrote gmsh2DMeshGenerator to make it thread-safe
*** PLEASE TEST THAT THE 2D MESH GENERATORS BEHAVE EXACTLY AS BEFORE ***
---
Mesh/meshGFace.cpp | 794 ++++++++++++++--------------------
Mesh/meshGFace.h | 6 -
Mesh/meshGFaceTransfinite.cpp | 52 +++
3 files changed, 370 insertions(+), 482 deletions(-)
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index aec0afa95e..b959b082c9 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -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];
@@ -107,24 +54,25 @@ static bool noseam(GFace *gf)
while (it != edges.end()){
GEdge *ge = *it ;
bool seam = ge->isSeam(gf);
- if (seam) return false;
+ if(seam) return false;
++it;
}
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();
deMeshGFace dem;
std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
- for (; itr != edgesNotRecovered.end(); ++itr){
+ 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()){
+ if((*it)->triangles.size() || (*it)->quadrangles.size()){
facesToRemesh.push_back(*it);
dem(*it);
}
@@ -134,21 +82,22 @@ static void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
int p1 = itr->p1;
int p2 = itr->p2;
int N = itr->ge->lines.size();
- GVertex * g1 = itr->ge->getBeginVertex();
- GVertex * g2 = itr->ge->getEndVertex();
+ GVertex *g1 = itr->ge->getBeginVertex();
+ GVertex *g2 = itr->ge->getEndVertex();
Range<double> bb = itr->ge->parBounds(0);
std::vector<MLine*> newLines;
- for (int i = 0; i < N; i++){
+ 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();
- else if (v1->onWhat() == g2) t1 = bb.high();
+ if(v1->onWhat() == g1) t1 = bb.low();
+ else if(v1->onWhat() == g2) t1 = bb.high();
else {
MEdgeVertex *ev1 = (MEdgeVertex*)v1;
lc1 = ev1->getLc();
@@ -156,8 +105,8 @@ static void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
}
double t2;
double lc2 = -1;
- if (v2->onWhat() == g1) t2 = bb.low();
- else if (v2->onWhat() == g2) t2 = bb.high();
+ if(v2->onWhat() == g1) t2 = bb.low();
+ else if(v2->onWhat() == g2) t2 = bb.high();
else {
MEdgeVertex *ev2 = (MEdgeVertex*)v2;
lc2 = ev2->getLc();
@@ -165,14 +114,14 @@ static void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
}
// periodic
- if (v1->onWhat() == g1 && v1->onWhat() == g2)
+ if(v1->onWhat() == g1 && v1->onWhat() == g2)
t1 = fabs(t2-bb.low()) < fabs(t2-bb.high()) ? bb.low() : bb.high();
- if (v2->onWhat() == g1 && v2->onWhat() == g2)
+ if(v2->onWhat() == g1 && v2->onWhat() == g2)
t2 = fabs(t1-bb.low()) < fabs(t1-bb.high()) ? bb.low() : bb.high();
- if (lc1 == -1)
+ if(lc1 == -1)
lc1 = BGM_MeshSize(v1->onWhat(), 0, 0, v1->x(), v1->y(), v1->z());
- if (lc2 == -1)
+ if(lc2 == -1)
lc2 = BGM_MeshSize(v2->onWhat(), 0, 0, v2->x(), v2->y(), v2->z());
// should be better, i.e. equidistant
double t = 0.5 * (t2 + t1);
@@ -190,7 +139,7 @@ static void remeshUnrecoveredEdges(std::set<EdgeToRecover> &edgesNotRecovered,
itr->ge->lines = newLines;
itr->ge->mesh_vertices.clear();
N = itr->ge->lines.size();
- for (int i = 1; i < N; i++){
+ for(int i = 1; i < N; i++){
itr->ge->mesh_vertices.push_back(itr->ge->lines[i]->getVertex(0));
}
}
@@ -223,44 +172,49 @@ 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;
- if (pass_ == 2){
- m->add_geom (ge->tag(), 1);
- g = m->get_geom(ge->tag(),1);
+ if(pass_ == 2){
+ m->add_geom(ge->tag(), 1);
+ g = m->get_geom(ge->tag(), 1);
}
- for (unsigned int i = 0; i < ge->lines.size(); i++){
+ 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));
+ MVertex *vend = ge->lines[i]->getVertex(1);
+ 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);
- 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;
- }
+ BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, 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;
+ // }
}
}
- if (pass_ == 2 && ge->getBeginVertex()){
+ 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,151 +222,41 @@ 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;
+ BDS_GeomEntity CLASS_F(1, 2);
+ 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
- if (gf->geomType() == GEntity::CompoundSurface){
- //printf("replace edges by compound lines \n");
+ // replace edges by their compounds
+ if(gf->geomType() == GEntity::CompoundSurface){
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());
}
-
-
- 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();
+ std::set<MVertex*> all_vertices;
+
+ std::list<GEdge*>::iterator it = edges.begin();
while(it != edges.end()){
- if ((*it)->isSeam(gf)) return false;
+ if((*it)->isSeam(gf)) return false;
if(!(*it)->isMeshDegenerated()){
- for (unsigned int i = 0; i< (*it)->lines.size(); i++){
+ for(unsigned int i = 0; i< (*it)->lines.size(); i++){
all_vertices.insert((*it)->lines[i]->getVertex(0));
all_vertices.insert((*it)->lines[i]->getVertex(1));
}
@@ -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,13 +279,15 @@ 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() );
+ (*itvx)->mesh_vertices.end());
++itvx;
}
- if (all_vertices.size() < 3){
+ if(all_vertices.size() < 3){
Msg::Warning("Mesh Generation of Model Face %d Skipped: "
"Only %d Mesh Vertices on The Contours",
gf->tag(), all_vertices.size());
@@ -448,136 +295,116 @@ 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.
- SVector3 dd (bbox.max(), bbox.min());
+ // 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;
+ 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
+ // 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);
-
+ // 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){
+ if(debug && RECUR_ITER == 0){
char name[245];
sprintf(name, "surface%d-initial-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0);
sprintf(name, "surface%d-initial-param.pos", gf->tag());
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,40 +414,42 @@ 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;
}
- if (edgesNotRecovered.size()){
+ if(edgesNotRecovered.size()){
std::ostringstream sstream;
- for (std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
+ for(std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
itr != edgesNotRecovered.end(); ++itr)
sstream << " " << itr->ge->tag();
Msg::Warning(":-( There are %d intersections in the 1D mesh (curves%s)",
edgesNotRecovered.size(), sstream.str().c_str());
Msg::Warning("8-| Gmsh splits those edges and tries again");
- if (debug){
+ 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);
+ if(repairSelfIntersecting1dMesh)
+ remeshUnrecoveredEdges(recoverMapInv, edgesNotRecovered, facesToRemesh);
else{
std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
- int *_error = new int[3*edgesNotRecovered.size()];
+ int *_error = new int[3 * edgesNotRecovered.size()];
int I = 0;
- for (; itr != edgesNotRecovered.end(); ++itr){
+ for(; itr != edgesNotRecovered.end(); ++itr){
int p1 = itr->p1;
int p2 = itr->p2;
int tag = itr->ge->tag();
- _error[3*I+0] = p1;
- _error[3*I+1] = p2;
- _error[3*I+2] = tag;
+ _error[3 * I + 0] = p1;
+ _error[3 * I + 1] = p2;
+ _error[3 * I + 2] = tag;
I++;
}
throw _error;
@@ -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);
+ if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
+ 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
@@ -652,11 +481,11 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
if(e->g && e->numfaces() == 2){
BDS_Point *oface[2];
e->oppositeof(oface);
- if (oface[0]->iD < 0){
+ if(oface[0]->iD < 0){
recur_tag(e->faces(1), &CLASS_F);
break;
}
- else if (oface[1]->iD < 0){
+ else if(oface[1]->iD < 0){
recur_tag(e->faces(0), &CLASS_F);
break;
}
@@ -668,37 +497,39 @@ 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();
- 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());
- }
- else if(ge->dim() == 1){
- double u;
- here->getParameter(0,u);
- pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
+ BDS_Point *pp = (BDS_Point*)doc.points[i].data;
+ if(recoverMap.count(pp)){
+ MVertex *here = recoverMap[pp];
+ GEntity *ge = here->onWhat();
+ if(ge->dim() == 0){
+ pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
+ }
+ else if(ge->dim() == 1){
+ double u;
+ here->getParameter(0, u);
+ pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
+ }
+ else
+ pp->lcBGM() = MAX_LC;
+ pp->lc() = pp->lcBGM();
}
- else
- 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();
@@ -707,13 +538,15 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
std::list<BDS_Edge*>::iterator ite = m->edges.begin();
while (ite != m->edges.end()){
BDS_Edge *e = *ite;
- if (e->numfaces() == 0)
+ if(e->numfaces() == 0)
m->del_edge(e);
else{
- if (!e->g)
+ 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;
}
@@ -724,7 +557,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
- if (debug){
+ if(debug){
char name[245];
sprintf(name, "surface%d-recovered-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0);
@@ -756,12 +589,13 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// fill the small gmsh structures
{
- std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
+ 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;
@@ -771,16 +605,21 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end()){
BDS_Face *t = *itt;
- if (!t->deleted){
+ 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])
- gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ 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));
}
}
@@ -792,14 +631,14 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// BDS mesh is passed in order not to recompute local coordinates of
// vertices
if(AlgoDelaunay2D(gf)){
- if (CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
+ if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
gmshBowyerWatsonFrontal(gf);
else
gmshBowyerWatson(gf);
- for (int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
+ for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
laplaceSmoothing(gf);
}
- else if (debug){
+ else if(debug){
char name[256];
sprintf(name, "real%d.pos", gf->tag());
outputScalarField(m->triangles, name, 0, gf);
@@ -809,10 +648,8 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// delete the mesh
delete m;
- delete [] U_;
- delete [] V_;
- if (gf->meshAttributes.recombine)
+ if(gf->meshAttributes.recombine)
gmshRecombineIntoQuads(gf);
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
@@ -832,21 +669,21 @@ static inline double dist2(const SPoint2 &p1, const SPoint2 &p2)
{
const double dx = p1.x() - p2.x();
const double dy = p1.y() - p2.y();
- return dx*dx+dy*dy;
+ return dx * dx + dy * dy;
}
static void printMesh1d(int iEdge, int seam, std::vector<SPoint2> &m)
{
printf("Mesh1D for edge %d seam %d\n", iEdge, seam);
- for (unsigned int i = 0; i < m.size(); i++){
+ for(unsigned int i = 0; i < m.size(); i++){
printf("%12.5E %12.5E\n", m[i].x(), m[i].y());
}
}
-static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
+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)
{
@@ -854,19 +691,19 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
// the edge points on the face for seams, we build the list for
// every side for closed loops, we build it on both senses
- std::map<GEntity*,std::vector<SPoint2> > meshes;
- std::map<GEntity*,std::vector<SPoint2> > meshes_seam;
+ std::map<GEntity*, std::vector<SPoint2> > meshes;
+ std::map<GEntity*, std::vector<SPoint2> > meshes_seam;
const int MYDEBUG = false;
- std::map<BDS_Point*,MVertex*> recover_map;
+ std::map<BDS_Point*, MVertex*> recoverMapLocal;
result.clear();
count = 0;
GEdgeLoop::iter it = gel.begin();
- if (MYDEBUG)
+ if(MYDEBUG)
printf("face %d with %d edges case %d\n", gf->tag(),
(int)gf->edges().size(), seam_the_first);
@@ -880,23 +717,23 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
Range<double> range = ges.ge->parBounds(0);
MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0];
- mesh1d.push_back(ges.ge->reparamOnFace(gf,range.low(),1));
- if (seam) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, range.low(), -1));
- for (unsigned int i = 0; i < ges.ge->mesh_vertices.size(); i++){
+ mesh1d.push_back(ges.ge->reparamOnFace(gf, range.low(), 1));
+ if(seam) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, range.low(), -1));
+ for(unsigned int i = 0; i < ges.ge->mesh_vertices.size(); i++){
double u;
here = ges.ge->mesh_vertices[i];
here->getParameter(0, u);
mesh1d.push_back(ges.ge->reparamOnFace(gf, u, 1));
- if ( seam ) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, u, -1));
+ if(seam) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, u, -1));
}
here = ges.ge->getEndVertex()->mesh_vertices[0];
mesh1d.push_back(ges.ge->reparamOnFace(gf, range.high(), 1));
- if ( seam ) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, range.high(), -1));
+ if(seam) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, range.high(), -1));
meshes.insert(std::pair<GEntity*,std::vector<SPoint2> >(ges.ge, mesh1d));
if(seam) meshes_seam.insert(std::pair<GEntity*,std::vector<SPoint2> >
(ges.ge, mesh1d_seam));
- // printMesh1d (ges.ge->tag(), seam, mesh1d);
- // if (seam)printMesh1d (ges.ge->tag(), seam, mesh1d_seam);
+ // printMesh1d(ges.ge->tag(), seam, mesh1d);
+ // if(seam) printMesh1d (ges.ge->tag(), seam, mesh1d_seam);
it++;
}
@@ -908,7 +745,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
int counter = 0;
while (unordered.size()){
- if (MYDEBUG)
+ if(MYDEBUG)
printf("unordered.size() = %d\n", (int)unordered.size());
std::list<GEdgeSigned>::iterator it = unordered.begin();
std::vector<SPoint2> coords;
@@ -921,13 +758,13 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
GEdge *ge = (*it).ge;
bool seam = ge->isSeam(gf);
mesh1d = meshes[ge];
- if (seam){ mesh1d_seam = meshes_seam[ge]; }
+ if(seam){ mesh1d_seam = meshes_seam[ge]; }
mesh1d_reversed.insert(mesh1d_reversed.begin(), mesh1d.rbegin(), mesh1d.rend());
- if (seam) mesh1d_seam_reversed.insert(mesh1d_seam_reversed.begin(),
- mesh1d_seam.rbegin(),mesh1d_seam.rend());
- if (!counter){
+ if(seam) mesh1d_seam_reversed.insert(mesh1d_seam_reversed.begin(),
+ mesh1d_seam.rbegin(),mesh1d_seam.rend());
+ if(!counter){
counter++;
- if (seam && seam_the_first){
+ if(seam && seam_the_first){
coords = ((*it)._sign == 1) ? mesh1d_seam : mesh1d_seam_reversed;
found = (*it);
Msg::Info("This test case would have failed in previous Gmsh versions ;-)");
@@ -937,19 +774,19 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
found = (*it);
}
unordered.erase(it);
- if (MYDEBUG)
+ if(MYDEBUG)
printf("Starting with edge = %d seam %d\n", (*it).ge->tag(), seam);
break;
}
else{
- if (MYDEBUG)
+ if(MYDEBUG)
printf("Followed by edge = %d\n", (*it).ge->tag());
SPoint2 first_coord = mesh1d[0];
double d = -1, d_reversed = -1, d_seam = -1, d_seam_reversed = -1;
d = dist2(last_coord, first_coord);
- if (MYDEBUG)
+ if(MYDEBUG)
printf("%g %g dist = %12.5E\n", first_coord.x(), first_coord.y(), d);
- if (d < tol){
+ if(d < tol){
coords.clear();
coords = mesh1d;
found = GEdgeSigned(1,ge);
@@ -958,22 +795,22 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
}
SPoint2 first_coord_reversed = mesh1d_reversed[0];
d_reversed = dist2(last_coord, first_coord_reversed);
- if (MYDEBUG)
+ if(MYDEBUG)
printf("%g %g dist_reversed = %12.5E\n",
- first_coord_reversed.x(), first_coord_reversed.y(),d_reversed);
- if (d_reversed < tol){
+ first_coord_reversed.x(), first_coord_reversed.y(), d_reversed);
+ if(d_reversed < tol){
coords.clear();
coords = mesh1d_reversed;
found = (GEdgeSigned(-1,ge));
unordered.erase(it);
goto Finalize;
}
- if (seam){
+ if(seam){
SPoint2 first_coord_seam = mesh1d_seam[0];
SPoint2 first_coord_seam_reversed = mesh1d_seam_reversed[0];
d_seam = dist2(last_coord,first_coord_seam);
- if (MYDEBUG) printf("dist_seam = %12.5E\n", d_seam);
- if (d_seam < tol){
+ if(MYDEBUG) printf("dist_seam = %12.5E\n", d_seam);
+ if(d_seam < tol){
coords.clear();
coords = mesh1d_seam;
found = (GEdgeSigned(1,ge));
@@ -981,8 +818,8 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
goto Finalize;
}
d_seam_reversed = dist2(last_coord, first_coord_seam_reversed);
- if (MYDEBUG) printf("dist_seam_reversed = %12.5E\n", d_seam_reversed);
- if (d_seam_reversed < tol){
+ if(MYDEBUG) printf("dist_seam_reversed = %12.5E\n", d_seam_reversed);
+ if(d_seam_reversed < tol){
coords.clear();
coords = mesh1d_seam_reversed;
found = GEdgeSigned(-1, ge);
@@ -995,27 +832,27 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
++it;
}
Finalize:
- if (MYDEBUG) printf("Finalize, found %d points\n", (int)coords.size());
- if (coords.size() == 0){
+ if(MYDEBUG) printf("Finalize, found %d points\n", (int)coords.size());
+ if(coords.size() == 0){
// 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;
}
std::vector<MVertex*> edgeLoop;
- if (found._sign == 1){
+ if(found._sign == 1){
edgeLoop.push_back(found.ge->getBeginVertex()->mesh_vertices[0]);
- for (unsigned int i = 0; i <found.ge->mesh_vertices.size(); i++)
+ for(unsigned int i = 0; i <found.ge->mesh_vertices.size(); i++)
edgeLoop.push_back(found.ge->mesh_vertices[i]);
}
else{
edgeLoop.push_back(found.ge->getEndVertex()->mesh_vertices[0]);
- for (int i = found.ge->mesh_vertices.size() - 1; i >= 0; i--)
+ for(int i = found.ge->mesh_vertices.size() - 1; i >= 0; i--)
edgeLoop.push_back(found.ge->mesh_vertices[i]);
}
@@ -1024,7 +861,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
found.ge->tag(), (int)edgeLoop.size(), (int)coords.size());
std::vector<BDS_Point*> edgeLoop_BDS;
- for (unsigned int i = 0; i < edgeLoop.size(); i++){
+ for(unsigned int i = 0; i < edgeLoop.size(); i++){
MVertex *here = edgeLoop[i];
GEntity *ge = here->onWhat();
double U, V;
@@ -1035,7 +872,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
if(ge->dim() == 0){
pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
}
- else if (ge->dim() == 1){
+ else if(ge->dim() == 1){
double u;
here->getParameter(0, u);
pp->lcBGM() = BGM_MeshSize(ge, u, 0,here->x(), here->y(), here->z());
@@ -1047,29 +884,29 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
m->add_geom (ge->tag(), ge->dim());
BDS_GeomEntity *g = m->get_geom(ge->tag(), ge->dim());
pp->g = g;
- if (MYDEBUG)
+ if(MYDEBUG)
printf("point %3d (%8.5f %8.5f : %8.5f %8.5f) (%2d,%2d)\n",
count, pp->u, pp->v, param.x(), param.y(), pp->g->classif_tag,
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];
- if (MYDEBUG) printf("last coord %g %g\n", last_coord.x(), last_coord.y());
+ if(MYDEBUG) printf("last coord %g %g\n", last_coord.x(), last_coord.y());
result.insert(result.end(), edgeLoop_BDS.begin(), edgeLoop_BDS.end());
}
// 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,27 +919,28 @@ 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;
{
- for (std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin();
- it != gf->edgeLoops.end(); it++){
+ for(std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin();
+ it != gf->edgeLoops.end(); it++){
std::vector<BDS_Point* > edgeLoop_BDS;
int nbPointsLocal;
const double fact[4] = {1.e-12, 1.e-7, 1.e-5, 1.e-3};
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;
@@ -1125,18 +963,16 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
{
DocRecord doc(nbPointsTotal + 4);
int count = 0;
- for (unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
+ 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++){
+ 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++;
@@ -1152,7 +988,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
bb[1] = new MVertex(bbox.min().x(), bbox.max().y(), 0, 0, -2);
bb[2] = new MVertex(bbox.max().x(), bbox.min().y(), 0, 0, -3);
bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, 0, -4);
- for (int ip = 0; ip < 4; ip++){
+ for(int ip = 0; ip < 4; ip++){
BDS_Point *pp = m->add_point(-ip - 1, bb[ip]->x(), bb[ip]->y(), gf);
m->add_geom(gf->tag(), 2);
BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
@@ -1162,7 +998,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
doc.points[nbPointsTotal+ip].adjacent = 0;
doc.points[nbPointsTotal+ip].data = pp;
}
- for (int ip = 0; ip < 4; ip++) delete bb[ip];
+ for(int ip = 0; ip < 4; ip++) delete bb[ip];
// Use "fast" inhouse recursive algo to generate the triangulation
// At this stage the triangulation is not what we need
@@ -1185,7 +1021,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
BDS_GeomEntity CLASS_F(1, 2);
BDS_GeomEntity CLASS_E(1, 1);
- if (debug){
+ if(debug){
char name[245];
sprintf(name, "surface%d-initial-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0);
@@ -1193,14 +1029,14 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
outputScalarField(m->triangles, name, 1);
}
- for (unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
+ 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);
- if (!e){
- Msg::Error("Impossible to recover the edge %d %d",
- edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
+ 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);
+ if(!e){
+ 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;
}
@@ -1221,11 +1057,11 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
if(e->g && e->numfaces () == 2){
BDS_Point *oface[2];
e->oppositeof(oface);
- if (oface[0]->iD < 0){
+ if(oface[0]->iD < 0){
recur_tag(e->faces(1), &CLASS_F);
break;
}
- else if (oface[1]->iD < 0){
+ else if(oface[1]->iD < 0){
recur_tag(e->faces(0), &CLASS_F);
break;
}
@@ -1239,7 +1075,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end()){
BDS_Face *t = *itt;
- if (!t->g){
+ if(!t->g){
m->del_face (t);
}
++itt;
@@ -1252,13 +1088,15 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
std::list<BDS_Edge*>::iterator ite = m->edges.begin();
while (ite != m->edges.end()){
BDS_Edge *e = *ite;
- if (e->numfaces() == 0)
+ if(e->numfaces() == 0)
m->del_edge(e);
else{
- if (!e->g)
+ 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;
}
@@ -1269,7 +1107,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
- if (debug){
+ if(debug){
char name[245];
sprintf(name, "surface%d-recovered-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0);
@@ -1278,11 +1116,11 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
}
// start mesh generation
- if (!AlgoDelaunay2D(gf)){
+ if(!AlgoDelaunay2D(gf)){
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;
@@ -1310,21 +1149,21 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end()){
BDS_Face *t = *itt;
- if (!t->deleted){
+ 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]];
- 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)
+ if(v1 != v2 && v1 != v3 && v2 != v3)
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));
}
}
@@ -1332,7 +1171,7 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
- if (debug){
+ if(debug){
char name[245];
sprintf(name, "surface%d-final-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0, gf);
@@ -1340,19 +1179,19 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
outputScalarField(m->triangles, name, 1);
}
- if (AlgoDelaunay2D(gf)){
- if (CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
+ if(AlgoDelaunay2D(gf)){
+ if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)
gmshBowyerWatsonFrontal(gf);
else
gmshBowyerWatson(gf);
- for (int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
+ for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
laplaceSmoothing(gf);
}
// delete the mesh
delete m;
- if (gf->meshAttributes.recombine)
+ if(gf->meshAttributes.recombine)
gmshRecombineIntoQuads(gf);
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
@@ -1375,15 +1214,13 @@ 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){
+ if(debugSurface >= 0 && gf->tag() != debugSurface){
gf->meshStatistics.status = GFace::DONE;
return;
}
@@ -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,16 +1255,18 @@ 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());
}
- // gmshQMorph(gf);
+ // gmshQMorph(gf);
Msg::Debug("Type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
@@ -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)){
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index bff73b2626..9ad9da65a0 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -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*> >
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
index fb1125f127..e2e9579625 100644
--- a/Mesh/meshGFaceTransfinite.cpp
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -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;
--
GitLab