diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index cc96dc3bec6c7e9d4455b62b29a9d14406c49801..4b7c68809b496d0bd4b528f1715ee2a82598afd5 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -15,6 +15,14 @@
GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
: GEdge(m, tag, 0 , 0), _compound(compound)
{
+
+ for (std::vector<GEdge*>::iterator it = compound.begin(); it != compound.end(); ++it){
+ if (!(*it)) {
+ Msg::Error("Incorrect edge in compound line %d\n", tag);
+ Msg::Exit(1);
+ }
+ }
+
orderEdges ();
int N = _compound.size();
v0 = _orientation[0] ? _compound[0]->getBeginVertex() : _compound[0]->getEndVertex();
@@ -75,7 +83,7 @@ void GEdgeCompound::orderEdges()
}
else{
Msg::Error("EdgeCompound %d is wrong (it has %d end points)",tag(),tempv.size());
- exit(1);
+ Msg::Exit(1);
}
//loop over all segments to order segments and store it in the list _c
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 3eb3687f4e0f936878a3d24550b0ae03a5f165a0..c162e107833c7cff4a8e17e923f5d2119c1fb376 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -174,14 +174,20 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U1, std::list<GEdge*> &V1)
: GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0)
{
+
+ for (std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
+ if (!(*it)){
+ Msg::Error("Incorrect face in compound surface %d\n", tag);
+ Msg::Exit(1);
+ }
+ }
+
getBoundingEdges();
if (!_U0.size()) _type = UNITCIRCLE;
else if (!_V1.size()) _type = UNITCIRCLE;
else _type = SQUARE;
- for (std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
- if (!(*it)) Msg::Error("Incorrect face in compound surface %d\n", tag);
- }
+
}
@@ -307,16 +313,16 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
vR = ge->mesh_vertices[j];
vR->getParameter(0,tR);
if (!vR->getParameter(0,tR)) {
- printf("vertex vr %p not medgevertex \n", vR);
- exit(1);
+ Msg::Error("vertex vr %p not medgevertex \n", vR);
+ Msg::Exit(1);
}
//printf("***tLoc=%g tL=%g, tR=%g L=%p (%g,%g,%g) et R= %p (%g,%g,%g)\n", tLoc, tL, tR, vL, vL->x(),vL->y(),vL->z(),vR, vR->x(), vR->y(), vR->z());
if (tLoc >= tL && tLoc <= tR){
found = true;
itR = coordinates.find(vR);
if (itR == coordinates.end()){
- printf("vertex %p (%g %g %g) not found\n", vR, vR->x(), vR->y(), vR->z());
- exit(1);
+ Msg::Error("vertex %p (%g %g %g) not found\n", vR, vR->x(), vR->y(), vR->z());
+ Msg::Exit(1);
}
break;
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index a61e647fd1482d0bdee43c08d4e1d4a3be0ad4e2..c1826521367104f6e5f1a921d0a41459fe74afdc 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -378,7 +378,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
std::list<GEdge*> edges = gf->edges();
// here, we will replace edges by their compounds
- //printf("***** In meshGFace: \n");
+ printf("***** In meshGFace: \n");
if (gf->geomType() == GEntity::CompoundSurface){
//printf("replace edges by compound lines \n");
std::set<GEdge*> mySet;
@@ -448,7 +448,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
double *V_ = new double[all_vertices.size()];
v_container::iterator itv = all_vertices.begin();
- //printf("boundary vertices size = %d \n", all_vertices.size());
+ printf("boundary vertices size = %d \n", all_vertices.size());
int count = 0;
SBoundingBox3d bbox;
@@ -458,7 +458,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
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());
+ //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);
@@ -539,7 +539,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
}
m->add_point(num, U, V, gf);
}
-
+
for(int i = 0; i < doc.numTriangles; i++) {
MVertex *V1 = (MVertex*)doc.points[doc.triangles[i].a].data;
@@ -548,10 +548,9 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
m->add_triangle(V1->getIndex(), V2->getIndex(), V3->getIndex());
}
-
// 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);
@@ -632,11 +631,12 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
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::Info("Boundary Edges recovered for surface %d", gf->tag());
+ Msg::Info("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
// inside the domain and, because all edges were recovered,
@@ -669,7 +669,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
++it;
}
// compute characteristic lengths at vertices
- Msg::Debug("Computing mesh size field at mesh vertices", 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();
@@ -730,15 +730,15 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
}
int nb_swap;
- // outputScalarField(m->triangles, "beforeswop.pos",1);
+ //outputScalarField(m->triangles, "beforeswop.pos",1);
Msg::Debug("Delaunizing the initial mesh");
gmshDelaunayizeBDS(gf, *m, nb_swap);
- // outputScalarField(m->triangles, "afterswop.pos",0)
+ //outputScalarField(m->triangles, "afterswop.pos",0);
Msg::Debug("Starting to add internal points");
// start mesh generation
if(!AlgoDelaunay2D(gf)){
- gmshRefineMeshBDS (gf,*m, CTX::instance()->mesh.refineSteps, 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);
@@ -789,7 +789,7 @@ 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);
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index ccac813653c78a1da65a1a6c044c3e21c5b80b9e..491ef1912f50bcc9dd20c91cc3c3e80950a3d058 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -449,9 +449,9 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
if (!e->deleted){
const double coord = 0.5;
BDS_Point *mid ;
- mid = m.add_point(++m.MAXPOINTNUMBER,
- coord * e->p1->u + (1 - coord) * e->p2->u,
- coord * e->p1->v + (1 - coord) * e->p2->v,gf);
+ double U = coord * e->p1->u + (1 - coord) * e->p2->u;
+ double V = coord * e->p1->v + (1 - coord) * e->p2->v;
+ mid = m.add_point(++m.MAXPOINTNUMBER, U , V, gf);
mid->lcBGM() = BGM_MeshSize
(gf,
(coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
@@ -534,11 +534,11 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const bool computeNodalSizeField)
{
+
int IT = 0;
int MAXNP = m.MAXPOINTNUMBER;
-
// classify correctly the embedded vertices
// use a negative model face number to avoid
// mesh motion
@@ -553,6 +553,7 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
p->lcBGM() = (*itvx)->prescribedMeshSizeAtVertex();
++itvx;
}
+
// IF ASKED , compute nodal size field using 1D Mesh
if (computeNodalSizeField){
@@ -597,6 +598,7 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
+
while (1){
if (NN2++ >= NN1)break;
if (!(*it)->deleted){
@@ -608,6 +610,8 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
}
++it;
}
+
+
if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT))) break;
double maxE = MAXE_;
double minE = MINE_;
@@ -623,7 +627,7 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
double t5 = Cpu();
smoothVertexPass(gf, m, nb_smooth, false);
double t6 = Cpu();
- // swapEdgePass ( gf, m, nb_swap);
+ swapEdgePass ( gf, m, nb_swap);
double t7 = Cpu();
// clean up the mesh
t_spl += t2 - t1;
@@ -632,13 +636,16 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
t_sw += t7 - t6;
t_col += t4 - t3;
t_sm += t6 - t5;
- m.cleanup();
+ m.cleanup();
+
IT++;
Msg::Debug(" iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : "
"%6d splits, %6d swaps, %6d collapses, %6d moves",
IT, minL, minE, maxL, maxE, nb_split, nb_swap, nb_collaps, nb_smooth);
if (nb_split == 0 && nb_collaps == 0) break;
+
}
+
double t_total = t_spl + t_sw + t_col + t_sm;
if(!t_total) t_total = 1.e-6;
diff --git a/benchmarks/2d/circle2.geo b/benchmarks/2d/circle2.geo
new file mode 100644
index 0000000000000000000000000000000000000000..8cd01468281a5e00901e0d008ac97b303a18f6f6
--- /dev/null
+++ b/benchmarks/2d/circle2.geo
@@ -0,0 +1,27 @@
+lc = 0.5;
+
+Point(1) = {0, 0, 0, lc};
+Point(2) = {0, 2, 0, lc};
+Point(3) = {2, 0, 0, lc};
+Point(4) = {0, -2, 0, lc};
+Point(5) = {-2, 0, 0, lc};
+Circle(1) = {2, 1, 5};
+Circle(2) = {5, 1, 4};
+Circle(3) = {4, 1, 3};
+Circle(4) = {3, 1, 2};
+
+Point(22) = {0, 1, 0, lc};
+Point(33) = {1, 0, 0, lc};
+Point(44) = {0, -1, 0, lc};
+Point(55) = {-1, 0, 0, lc};
+Circle(11) = {22, 1, 55};
+Circle(22) = {55, 1, 44};
+Circle(33) = {44, 1, 33};
+Circle(44) = {33, 1, 22};
+
+
+Line Loop(45) = {44, 11, 22, 33};
+Line Loop(46) = {4, 1, 2, 3};
+Plane Surface(47) = {45, 46};
+
+
diff --git a/benchmarks/2d/circle2_GEO.geo b/benchmarks/2d/circle2_GEO.geo
new file mode 100644
index 0000000000000000000000000000000000000000..2d493a82fd0275143b127b2bc084df7814ffd4d2
--- /dev/null
+++ b/benchmarks/2d/circle2_GEO.geo
@@ -0,0 +1,9 @@
+Mesh.CharacteristicLengthFactor=0.8;
+
+Merge "circle2.msh";
+CreateTopology;
+
+Compound Line (100) = {1,2,3,4};
+Compound Line (200) = {11,22,33,44};
+
+Compound Surface(1000)={47} Boundary {{1,2,3,4}, {11,22,33,44}};
diff --git a/benchmarks/2d/square_CLASS_GEO.geo b/benchmarks/2d/square_CLASS_GEO.geo
index 52d4b52fca2ec34d99a57b1b338bdf5ba01ab088..c584ddcd10174ea338150ec30000144211832c7d 100644
--- a/benchmarks/2d/square_CLASS_GEO.geo
+++ b/benchmarks/2d/square_CLASS_GEO.geo
@@ -3,11 +3,6 @@ Mesh.CharacteristicLengthFactor=6;
Merge "square_CLASS.msh";
CreateTopology;
-//Compound Line(100)={1};
-//Compound Line(200)={2};
-//Compound Line(300)={3};
-//Compound Line(400)={4};
-
Compound Line(150)={1,2};
Compound Line(160)={3,4};