From a3711e5e8f73221b8d04be8506cc8c91946cd0b1 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 27 May 2014 15:34:43 +0000
Subject: [PATCH] deep modification in the 2D mesh. Those modifs have been
prudently put between #if 0 #endif Those should allow to produce better
quality quad meshes and more robust triangular meshes (random factor has
disappeared) : Technically, the divide and conquer algo has been replaced by
a bowyer-watson with robust predicates.
---
Common/CommandLine.cpp | 12 +-
Common/DefaultOptions.h | 3 -
Common/Options.cpp | 7 -
Common/Options.h | 1 -
Mesh/BackgroundMesh.cpp | 11 +-
Mesh/BackgroundMesh.h | 2 +
Mesh/Generator.cpp | 72 +--
Mesh/meshGEdge.cpp | 14 +-
Mesh/meshGFace.cpp | 708 +++++++++++++++++---------
Mesh/meshGFace.h | 5 +-
Mesh/meshGFaceDelaunayInsertion.cpp | 448 ++++++++++++++--
Mesh/meshGFaceDelaunayInsertion.h | 15 +-
Mesh/meshGFaceOptimize.cpp | 51 +-
Mesh/meshGFaceOptimize.h | 3 +-
Mesh/meshGRegionDelaunayInsertion.cpp | 2 +-
Numeric/Numeric.cpp | 4 +-
Numeric/Numeric.h | 4 +-
benchmarks/2d/Square-01.geo | 2 +-
18 files changed, 940 insertions(+), 424 deletions(-)
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 65df25d005..60eefc288b 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -78,7 +78,7 @@ std::vector<std::pair<std::string, std::string> > GetUsage()
s.push_back(mp("-saveall", "Save all elements (discard physical group definitions)"));
s.push_back(mp("-parametric", "Save vertices with their parametric coordinates"));
s.push_back(mp("-algo string", "Select mesh algorithm (meshadapt, del2d, front2d, "
- "delquad, del3d, front3d, mmg3d)"));
+ "delquad, del3d, front3d, mmg3d, pack)"));
s.push_back(mp("-smooth int", "Set number of mesh smoothing steps"));
s.push_back(mp("-order int", "Set mesh order (1, ..., 5)"));
s.push_back(mp("-optimize[_netgen]", "Optimize quality of tetrahedral elements"));
@@ -709,16 +709,6 @@ void GetOptions(int argc, char *argv[])
else
Msg::Fatal("Missing number");
}
- else if(!strcmp(argv[i] + 1, "mpass")) {
- i++;
- if(argv[i]) {
- CTX::instance()->mesh.multiplePasses = atoi(argv[i++]);
- if(CTX::instance()->mesh.multiplePasses <= 0)
- Msg::Fatal("Number of Mesh Passes must be > 0");
- }
- else
- Msg::Fatal("Missing number");
- }
else if(!strcmp(argv[i] + 1, "ignorePartBound")) {
i++;
opt_mesh_ignore_part_bound(0, GMSH_SET, 1);
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 21a292c26a..55233156ca 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1037,9 +1037,6 @@ StringXNumber MeshOptions_Number[] = {
"Version of the MSH file format to use" },
{ F|O, "MshFilePartitioned" , opt_mesh_msh_file_partitioned , 0. ,
"Split MSH file by mesh partition (0: no, 1: yes, 2: create physicals by partition)" },
- { F|O, "MultiplePassesMeshes" , opt_mesh_multiple_passes , 0. ,
- "Do a first simple mesh and use it for complex background meshes (curvatures...)" },
-
{ F|O, "PartitionHexWeight" , opt_mesh_partition_hex_weight , 1 ,
"Weight of hexahedral element for METIS load balancing" },
{ F|O, "PartitionPrismWeight" , opt_mesh_partition_pri_weight , 1 ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 47165007fc..4be3ff536a 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5717,13 +5717,6 @@ double opt_mesh_second_order_experimental(OPT_ARGS_NUM)
return CTX::instance()->mesh.secondOrderExperimental;
}
-double opt_mesh_multiple_passes(OPT_ARGS_NUM)
-{
- if(action & GMSH_SET)
- CTX::instance()->mesh.multiplePasses = (int)val;
- return CTX::instance()->mesh.multiplePasses;
-}
-
double opt_mesh_second_order_linear(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index fd7c9ce9ab..9f79eed690 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -452,7 +452,6 @@ double opt_mesh_mesh_only_visible(OPT_ARGS_NUM);
double opt_mesh_min_circ_points(OPT_ARGS_NUM);
double opt_mesh_allow_swap_edge_angle(OPT_ARGS_NUM);
double opt_mesh_min_curv_points(OPT_ARGS_NUM);
-double opt_mesh_multiple_passes(OPT_ARGS_NUM);
double opt_mesh_order(OPT_ARGS_NUM);
double opt_mesh_ho_optimize(OPT_ARGS_NUM);
double opt_mesh_ho_nlayers(OPT_ARGS_NUM);
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index b215630f68..59c31b377c 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -472,6 +472,7 @@ void backgroundMesh::unset()
_current = 0;
}
+double backgroundMesh::sizeFactor = 1.0;
backgroundMesh::backgroundMesh(GFace *_gf, bool cfd)
#if defined(HAVE_ANN)
: _octree(0), uv_kdtree(0), nodes(0), angle_nodes(0), angle_kdtree(0)
@@ -836,21 +837,21 @@ void backgroundMesh::updateSizes(GFace *_gf)
MVertex *v = _2Dto3D[itv->first];
double lc;
if (v->onWhat()->dim() == 0){
- lc = BGM_MeshSize(v->onWhat(), 0,0,v->x(),v->y(),v->z());
+ lc = sizeFactor * BGM_MeshSize(v->onWhat(), 0,0,v->x(),v->y(),v->z());
}
else if (v->onWhat()->dim() == 1){
double u;
v->getParameter(0, u);
- lc = BGM_MeshSize(v->onWhat(), u, 0, v->x(), v->y(), v->z());
+ lc = sizeFactor * BGM_MeshSize(v->onWhat(), u, 0, v->x(), v->y(), v->z());
}
else{
reparamMeshVertexOnFace(v, _gf, p);
- lc = BGM_MeshSize(_gf, p.x(), p.y(), v->x(), v->y(), v->z());
+ lc = sizeFactor * BGM_MeshSize(_gf, p.x(), p.y(), v->x(), v->y(), v->z());
}
// printf("2D -- %g %g 3D -- %g %g\n",p.x(),p.y(),v->x(),v->y());
itv->second = std::min(lc,itv->second);
- itv->second = std::max(itv->second, CTX::instance()->mesh.lcMin);
- itv->second = std::min(itv->second, CTX::instance()->mesh.lcMax);
+ itv->second = std::max(itv->second, sizeFactor * CTX::instance()->mesh.lcMin);
+ itv->second = std::min(itv->second, sizeFactor * CTX::instance()->mesh.lcMax);
}
// do not allow large variations in the size field
// (Int. J. Numer. Meth. Engng. 43, 1143-1165 (1998) MESH GRADATION
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 6550bf69d9..a46f5b1ebc 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -41,6 +41,7 @@ struct crossField2d
class backgroundMesh : public simpleFunction<double>
{
+ static double sizeFactor;
MElementOctree *_octree;
std::vector<MVertex*> _vertices;
std::vector<MElement*> _triangles;
@@ -66,6 +67,7 @@ class backgroundMesh : public simpleFunction<double>
static void setCrossFieldsByDistance(GFace *);
static void unset();
static backgroundMesh *current () { return _current; }
+ static void setSizeFactor (double s) {sizeFactor = s;}
void propagate1dMesh(GFace *);
void propagateCrossField(GFace *);
void propagateCrossFieldByDistance(GFace *);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index a746e4136e..1c02a709a9 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -106,71 +106,6 @@ static void buildASetOfEquivalentMeshVertices(GFace *gf,
}
}
-struct geomThresholdVertexEquivalence
-{
- // Initial MVertex associated to one given MVertex
- std::map<GVertex*, MVertex*> backward_map;
- // initiate the forward and backward maps
- geomThresholdVertexEquivalence(GModel *g);
- // restores the initial state
- ~geomThresholdVertexEquivalence ();
-};
-
-geomThresholdVertexEquivalence::geomThresholdVertexEquivalence(GModel *g)
-{
- std::multimap<MVertex*, MVertex*> equivalenceMap;
- for (GModel::fiter it = g->firstFace(); it != g->lastFace(); ++it)
- buildASetOfEquivalentMeshVertices(*it, equivalenceMap, backward_map);
- // build the structure that identifiate geometrically equivalent
- // mesh vertices.
- for (std::map<GVertex*, MVertex*>::iterator it = backward_map.begin();
- it != backward_map.end(); ++it){
- GVertex *g = it->first;
- MVertex *v = it->second;
- MVertex *other = isEquivalentTo(equivalenceMap, v);
- if (v != other){
- printf("Finally : %d equivalent to %d\n", v->getNum(), other->getNum());
- g->mesh_vertices.clear();
- g->mesh_vertices.push_back(other);
- std::list<GEdge*> ed = g->edges();
- for (std::list<GEdge*>::iterator ite = ed.begin() ; ite != ed.end() ; ++ite){
- std::vector<MLine*> newl;
- for (unsigned int i = 0; i < (*ite)->lines.size(); ++i){
- MLine *l = (*ite)->lines[i];
- MVertex *v1 = l->getVertex(0);
- MVertex *v2 = l->getVertex(1);
- if (v1 == v && v2 != other){
- delete l;
- l = new MLine(other,v2);
- newl.push_back(l);
- }
- else if (v1 != other && v2 == v){
- delete l;
- l = new MLine(v1,other);
- newl.push_back(l);
- }
- else if (v1 != v && v2 != v)
- newl.push_back(l);
- else
- delete l;
- }
- (*ite)->lines = newl;
- }
- }
- }
-}
-
-geomThresholdVertexEquivalence::~geomThresholdVertexEquivalence()
-{
- // restore the initial data
- for (std::map<GVertex*, MVertex*>::iterator it = backward_map.begin();
- it != backward_map.end() ; ++it){
- GVertex *g = it->first;
- MVertex *v = it->second;
- g->mesh_vertices.clear();
- g->mesh_vertices.push_back(v);
- }
-}
template<class T>
static void GetQualityMeasure(std::vector<T*> &ele,
@@ -466,9 +401,6 @@ static void Mesh2D(GModel *m)
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
(*it)->meshStatistics.status = GFace::PENDING;
- // skip short mesh edges
- //geomThresholdVertexEquivalence inst(m);
-
// boundary layers are special: their generation (including vertices
// and curve meshes) is global as it depends on a smooth normal
// field generated from the surface mesh of the source surfaces
@@ -493,7 +425,7 @@ static void Mesh2D(GModel *m)
for(size_t K = 0 ; K < temp.size() ; K++){
if (temp[K]->meshStatistics.status == GFace::PENDING){
backgroundMesh::current()->unset();
- meshGFace mesher(true, CTX::instance()->mesh.multiplePasses);
+ meshGFace mesher(true);
mesher(temp[K]);
#if defined(HAVE_BFGS)
@@ -532,7 +464,7 @@ static void Mesh2D(GModel *m)
it != cf.end(); ++it){
if ((*it)->meshStatistics.status == GFace::PENDING){
backgroundMesh::current()->unset();
- meshGFace mesher(true, CTX::instance()->mesh.multiplePasses);
+ meshGFace mesher(true);
mesher(*it);
#if defined(HAVE_BFGS)
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 2ef1464367..617ce1ffc6 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -317,6 +317,12 @@ static void printFandPrimitive(int tag, std::vector<IntPoint> &Points)
}
*/
+// new algo for recombining + splitting
+static int increaseN (int N){
+ if (((N+1)/2 - 1) % 2 != 0) return N+2;
+ return N;
+}
+
void meshGEdge::operator() (GEdge *ge)
{
#if defined(HAVE_ANN)
@@ -417,15 +423,17 @@ void meshGEdge::operator() (GEdge *ge)
// force odd number of points if blossom is used for recombination
if((ge->meshAttributes.method != MESH_TRANSFINITE ||
CTX::instance()->mesh.flexibleTransfinite) &&
- CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0){
+ CTX::instance()->mesh.algoRecombine == 1){
if(CTX::instance()->mesh.recombineAll){
- N++;
+ if (N % 2 == 0) N++;
+ N = increaseN(N);
}
else{
std::list<GFace*> faces = ge->faces();
for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
if((*it)->meshAttributes.recombine){
- N++;
+ if (N % 2 == 0) N ++;
+ N = increaseN(N);
break;
}
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index beae02b54c..57e38900b6 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -46,6 +46,30 @@
#include "boundaryLayersData.h"
#include "filterElements.h"
+// new quad code
+// new initial delaunay
+#define OLD_CODE 1
+
+static void computeElementShapes(GFace *gf, double &worst, double &avg,
+ double &best, int &nT, int &greaterThan)
+{
+ worst = 1.e22;
+ avg = 0.0;
+ best = 0.0;
+ nT = 0;
+ greaterThan = 0;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ double q = qmTriangle(gf->triangles[i], QMTRI_RHO);
+ if(q > .9) greaterThan++;
+ avg += q;
+ worst = std::min(worst, q);
+ best = std::max(best, q);
+ nT++;
+ }
+ avg /= nT;
+}
+
+
inline double myAngle(const SVector3 &a, const SVector3 &b, const SVector3 &d)
{
double cosTheta = dot(a, b);
@@ -53,6 +77,147 @@ inline double myAngle(const SVector3 &a, const SVector3 &b, const SVector3 &d)
return atan2(sinTheta, cosTheta);
}
+class quadMeshRemoveHalfOfOneDMesh
+{
+ GFace *_gf;
+public:
+ std::map<GEdge*,std::vector<MLine*> > _backup;
+ std::map<MEdge, MVertex*,Less_Edge> _middle;
+ // remove one point every two and remember middle points
+ quadMeshRemoveHalfOfOneDMesh (GFace* gf) : _gf(gf){
+ // only do it if a recombination has to be done
+ if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine){
+ // printf("GFace %d removing half of the points in the 1D mesh\n",gf->tag());
+ std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*>::iterator ite = edges.begin();
+ while(ite != edges.end()){
+ if(!(*ite)->isMeshDegenerated()){
+ std::vector<MLine*> temp;
+ (*ite)->mesh_vertices.clear();
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i+=2){
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ MVertex *v3 = (*ite)->lines[i+1]->getVertex(1);
+ temp.push_back(new MLine(v1,v3));
+ if (v1->onWhat() == *ite) (*ite)->mesh_vertices.push_back(v1);
+ _middle[MEdge(v1,v3)] = v2;
+ }
+ _backup[*ite] = (*ite)->lines;
+ // printf("line %d goes from %d to %d\n",(*ite)->tag(), (*ite)->lines.size()-1,temp.size()-1);
+ (*ite)->lines = temp;
+ }
+ ++ite;
+ }
+ }
+ backgroundMesh::setSizeFactor(2.0);
+ }
+ void subdivide (){
+ std::vector<MQuadrangle*> qnew;
+
+ std::map<MEdge,MVertex*,Less_Edge> eds;
+
+ for(unsigned int i=0;i<_gf->triangles.size();i++){
+ MVertex *v[3];
+ SPoint2 m[3];
+ for (int j=0;j<3;j++){
+ MEdge E = _gf->triangles[i]->getEdge(j);
+ SPoint2 p1, p2;
+ reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
+ std::map<MEdge, MVertex *, Less_Edge>::iterator it = _middle.find(E);
+ std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
+ m[j] = p1;
+ if (it == _middle.end() && it2 == eds.end()){
+ GPoint gp = _gf->point((p1+p2)*0.5);
+ v[j] = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ _gf->mesh_vertices.push_back(v[j]);
+ eds[E] = v[j];
+ }
+ else if (it == _middle.end()){
+ v[j] = it2->second;
+ }
+ else {
+ v[j] = it->second;
+ v[j]->onWhat()->mesh_vertices.push_back(v[j]);
+ }
+ }
+ GPoint gp = _gf->point((m[0]+m[1]+m[2])*(1./3.));
+ MFaceVertex *vmid = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ _gf->mesh_vertices.push_back(vmid);
+ qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(0),v[0],vmid,v[2]));
+ qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(1),v[1],vmid,v[0]));
+ qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(2),v[2],vmid,v[1]));
+ delete _gf->triangles[i];
+ }
+ _gf->triangles.clear();
+ for(unsigned int i=0;i<_gf->quadrangles.size();i++){
+ MVertex *v[4];
+ SPoint2 m[4];
+ for (int j=0;j<4;j++){
+ MEdge E = _gf->quadrangles[i]->getEdge(j);
+ SPoint2 p1, p2;
+ reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
+ std::map<MEdge, MVertex *>::iterator it = _middle.find(E);
+ std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
+ m[j] = p1;
+ if (it == _middle.end() && it2 == eds.end()){
+ GPoint gp = _gf->point((p1+p2)*0.5);
+ v[j] = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ _gf->mesh_vertices.push_back(v[j]);
+ eds[E] = v[j];
+ }
+ else if (it == _middle.end()){
+ v[j] = it2->second;
+ }
+ else {
+ v[j] = it->second;
+ v[j]->onWhat()->mesh_vertices.push_back(v[j]);
+ }
+ }
+ GPoint gp = _gf->point((m[0]+m[1]+m[2]+m[3])*0.25);
+ MVertex *vmid = new MFaceVertex (gp.x(),gp.y(),gp.z(),_gf,gp.u(),gp.v());
+ _gf->mesh_vertices.push_back(vmid);
+ qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(0),v[0],vmid,v[3]));
+ qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(1),v[1],vmid,v[0]));
+ qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(2),v[2],vmid,v[1]));
+ qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(3),v[3],vmid,v[2]));
+ delete _gf->quadrangles[i];
+ }
+ _gf->quadrangles = qnew;
+ // printf("%d triangles %d quads\n",_gf->triangles.size(),_gf->quadrangles.size());
+ }
+ void finish (){
+ backgroundMesh::setSizeFactor(1.0);
+ if(CTX::instance()->mesh.recombineAll || _gf->meshAttributes.recombine){
+ // recombine the elements on the half mesh
+ recombineIntoQuads(_gf,true,true,.1);
+ Msg::Info("subdividing");
+ // _gf->model()->writeMSH("hop1.msh");
+ subdivide();
+ // _gf->model()->writeMSH("hop2.msh");
+ restore();
+ recombineIntoQuads(_gf,true,true,0.1);
+ computeElementShapes(_gf,
+ _gf->meshStatistics.worst_element_shape,
+ _gf->meshStatistics.average_element_shape,
+ _gf->meshStatistics.best_element_shape,
+ _gf->meshStatistics.nbTriangle,
+ _gf->meshStatistics.nbGoodQuality);
+ }
+ }
+ void restore (){
+ std::list<GEdge*> edges = _gf->edges();
+ std::list<GEdge*>::iterator ite = edges.begin();
+ while(ite != edges.end()){
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ delete (*ite)->lines[i];
+ }
+ (*ite)->lines = _backup[*ite];
+ ++ite;
+ }
+ }
+};
+
+
struct myPlane {
SPoint3 p;
SVector3 n;
@@ -464,25 +629,6 @@ static bool algoDelaunay2D(GFace *gf)
return false;
}
-static void computeElementShapes(GFace *gf, double &worst, double &avg,
- double &best, int &nT, int &greaterThan)
-{
- worst = 1.e22;
- avg = 0.0;
- best = 0.0;
- nT = 0;
- greaterThan = 0;
- for(unsigned int i = 0; i < gf->triangles.size(); i++){
- double q = qmTriangle(gf->triangles[i], QMTRI_RHO);
- if(q > .9) greaterThan++;
- avg += q;
- worst = std::min(worst, q);
- best = std::max(best, q);
- nT++;
- }
- avg /= nT;
-}
-
static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
std::map<MVertex*, BDS_Point*> &recoverMapInv,
std::set<EdgeToRecover> *e2r,
@@ -923,6 +1069,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
std::map<MVertex*, BDS_Point*> recoverMapInv;
std::list<GEdge*> edges = replacement_edges ? *replacement_edges : gf->edges();
std::list<int> dir = gf->edgeOrientations();
+ std::vector<MEdge> medgesToRecover;
// replace edges by their compounds
// if necessary split compound and remesh parts
@@ -938,6 +1085,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
while(ite != edges.end()){
if((*ite)->isSeam(gf)) return false;
if(!(*ite)->isMeshDegenerated()){
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i++)
+ medgesToRecover.push_back(MEdge((*ite)->lines[i]->getVertex(0),
+ (*ite)->lines[i]->getVertex(1)));
for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
MVertex *v1 = (*ite)->lines[i]->getVertex(0);
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
@@ -967,6 +1117,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
std::list<GEdge*> emb_edges = gf->embeddedEdges();
ite = emb_edges.begin();
while(ite != emb_edges.end()){
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i++)
+ medgesToRecover.push_back(MEdge((*ite)->lines[i]->getVertex(0),
+ (*ite)->lines[i]->getVertex(1)));
if(!(*ite)->isMeshDegenerated()){
all_vertices.insert((*ite)->mesh_vertices.begin(),
(*ite)->mesh_vertices.end() );
@@ -1036,7 +1189,6 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
bbox += SPoint3(param[0], param[1], 0);
count++;
}
- all_vertices.clear();
// here check if some boundary layer nodes should be added
@@ -1049,30 +1201,31 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// 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(points.size() + 4);
+#ifdef OLD_CODE
{
+ 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)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
// printf("%22.15E %22.15E \n",XX,YY);
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 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()}};
+ {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++){
BDS_Point *pp = m->add_point(-ip - 1, bb[ip][0], bb[ip][1], gf);
m->add_geom(gf->tag(), 2);
@@ -1083,12 +1236,12 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
doc.points[points.size() + ip].adjacent = 0;
doc.points[points.size() + ip].data = pp;
}
-
+
// 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)
+ // -) It contains triangles outside the domain (the first edge
+ // loop is the outer one)
Msg::Debug("Meshing of the convex hull (%d points)", points.size());
try{
doc.MakeMeshWithPoints();
@@ -1097,198 +1250,245 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
Msg::Error("%s", err);
}
Msg::Debug("Meshing of the convex hull (%d points) done", points.size());
-
+
for(int i = 0; i < doc.numTriangles; i++) {
int a = doc.triangles[i].a;
int b = doc.triangles[i].b;
int c = doc.triangles[i].c;
int n = doc.numPoints;
if(a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n){
- Msg::Warning("Skipping bad triangle %d", i);
- continue;
+ Msg::Warning("Skipping bad triangle %d", i);
+ continue;
}
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);
}
-
- 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());
- std::set<EdgeToRecover> edgesToRecover;
- std::set<EdgeToRecover> edgesNotRecovered;
- ite = edges.begin();
- while(ite != edges.end()){
- if(!(*ite)->isMeshDegenerated())
- recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
- ++ite;
+ }
+#else
+ {
+ std::vector<MVertex*> v;
+ std::vector<MTriangle*> result;
+ v.insert(v.end(),all_vertices.begin(),all_vertices.end());
+
+ std::map<MVertex*,SPoint3> pos;
+ for(unsigned int i = 0; i < v.size(); i++) {
+ MVertex *v0 = v[i];
+ BDS_Point *p0 = recoverMapInv[v0];
+ pos[v0] = SPoint3(v0->x(),v0->y(),v0->z());
+ v0->setXYZ(p0->u,p0->v,0.0);
+ }
+ delaunayMeshIn2D(v, result, 0);
+ // delaunayMeshIn2D(v, result, 0, & medgesToRecover);
+
+ for(unsigned int i = 0; i < v.size()-4; i++) {
+ MVertex *v0 = v[i];
+ SPoint3 pp = pos[v0];
+ v0->setXYZ(pp.x(),pp.y(),pp.z());
+ }
+
+ // add the four corners
+ for(int ip = 0; ip < 4; ip++){
+ MVertex *vv = v[v.size()-ip-1];
+ BDS_Point *pp = m->add_point(-ip - 1, vv->x(),vv->y(), gf);
+ m->add_geom(gf->tag(), 2);
+ recoverMapInv[vv] = pp;
+ BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
+ pp->g = g;
}
- ite = emb_edges.begin();
- while(ite != emb_edges.end()){
- if(!(*ite)->isMeshDegenerated())
- recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
- ++ite;
+ // add the triangles
+ for(unsigned int i = 0; i < result.size(); i++) {
+ MVertex *v0 = result[i]->getVertex(0);
+ MVertex *v1 = result[i]->getVertex(1);
+ MVertex *v2 = result[i]->getVertex(2);
+ BDS_Point *p0 = recoverMapInv[v0];
+ BDS_Point *p1 = recoverMapInv[v1];
+ BDS_Point *p2 = recoverMapInv[v2];
+ m->add_triangle(p0->iD, p1->iD, p2->iD);
}
+ }
+#endif
-
+ 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());
+ std::set<EdgeToRecover> edgesToRecover;
+ std::set<EdgeToRecover> edgesNotRecovered;
+ ite = edges.begin();
+ while(ite != edges.end()){
+ if(!(*ite)->isMeshDegenerated())
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ ++ite;
+ }
+ ite = emb_edges.begin();
+ while(ite != emb_edges.end()){
+ if(!(*ite)->isMeshDegenerated())
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ ++ite;
+ }
+
+
// effectively recover the medge
- ite = edges.begin();
- while(ite != edges.end()){
- if(!(*ite)->isMeshDegenerated()){
- if (!recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2)){
- delete m;
- gf->meshStatistics.status = GFace::FAILED;
- return false;
- }
+ ite = edges.begin();
+ while(ite != edges.end()){
+ if(!(*ite)->isMeshDegenerated()){
+ if (!recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2)){
+ delete m;
+ gf->meshStatistics.status = GFace::FAILED;
+ return false;
}
- ++ite;
}
-
- Msg::Debug("Recovering %d mesh Edges (%d not recovered)", edgesToRecover.size(),
- edgesNotRecovered.size());
-
- if(edgesNotRecovered.size()){
- std::ostringstream sstream;
- 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());
- if (repairSelfIntersecting1dMesh)
- Msg::Warning("8-| Gmsh splits those edges and tries again");
-
- if(debug){
- char name[245];
- 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(recoverMapInv, edgesNotRecovered, facesToRemesh);
- else{
- std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
- //int *_error = new int[3 * edgesNotRecovered.size()];
- int I = 0;
- for(; itr != edgesNotRecovered.end(); ++itr){
- int p1 = itr->p1;
- int p2 = itr->p2;
- int tag = itr->ge->tag();
- Msg::Error("Edge not recovered: %d %d %d", p1, p2, tag);
- //_error[3 * I + 0] = p1;
- //_error[3 * I + 1] = p2;
- //_error[3 * I + 2] = tag;
- I++;
- }
- //throw _error;
+ ++ite;
+ }
+
+ Msg::Debug("Recovering %d mesh Edges (%d not recovered)", edgesToRecover.size(),
+ edgesNotRecovered.size());
+
+ if(edgesNotRecovered.size()){
+ std::ostringstream sstream;
+ 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());
+ if (repairSelfIntersecting1dMesh)
+ Msg::Warning("8-| Gmsh splits those edges and tries again");
+
+ if(debug){
+ char name[245];
+ 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(recoverMapInv, edgesNotRecovered, facesToRemesh);
+ else{
+ std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
+ //int *_error = new int[3 * edgesNotRecovered.size()];
+ int I = 0;
+ for(; itr != edgesNotRecovered.end(); ++itr){
+ int p1 = itr->p1;
+ int p2 = itr->p2;
+ int tag = itr->ge->tag();
+ Msg::Error("Edge not recovered: %d %d %d", p1, p2, tag);
+ //_error[3 * I + 0] = p1;
+ //_error[3 * I + 1] = p2;
+ //_error[3 * I + 2] = tag;
+ I++;
}
-
- // delete the mesh
- delete m;
- if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
- return meshGenerator
- (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh,
- debug, replacement_edges);
- return false;
+ //throw _error;
}
-
- if(RECUR_ITER > 0)
- 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 a triangle that has a negative node and recursively
- // tag all exterior triangles
- {
- std::list<BDS_Face*>::iterator itt = m->triangles.begin();
- while (itt != m->triangles.end()){
- (*itt)->g = 0;
- ++itt;
+
+ // delete the mesh
+ delete m;
+ if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
+ return meshGenerator
+ (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh,
+ debug, replacement_edges);
+ return false;
+ }
+
+ if(RECUR_ITER > 0)
+ 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 a triangle that has a negative node and recursively
+ // tag all exterior triangles
+ {
+ std::list<BDS_Face*>::iterator itt = m->triangles.begin();
+ while (itt != m->triangles.end()){
+ (*itt)->g = 0;
+ ++itt;
+ }
+ itt = m->triangles.begin();
+ while (itt != m->triangles.end()){
+ BDS_Face *t = *itt;
+ BDS_Point *n[4];
+ t->getNodes(n);
+ if (n[0]->iD < 0 || n[1]->iD < 0 ||
+ n[2]->iD < 0 ) {
+ recur_tag(t, &CLASS_EXTERIOR);
+ break;
}
- itt = m->triangles.begin();
- while (itt != m->triangles.end()){
- BDS_Face *t = *itt;
- BDS_Point *n[4];
- t->getNodes(n);
- if (n[0]->iD < 0 || n[1]->iD < 0 ||
- n[2]->iD < 0 ) {
- recur_tag(t, &CLASS_EXTERIOR);
+ ++itt;
+ }
+ }
+
+ // now find an edge that has belongs to one of the exterior
+ // triangles
+ {
+ std::list<BDS_Edge*>::iterator ite = m->edges.begin();
+ while (ite != m->edges.end()){
+ BDS_Edge *e = *ite;
+ if(e->g && e->numfaces() == 2){
+ if(e->faces(0)->g == &CLASS_EXTERIOR){
+ recur_tag(e->faces(1), &CLASS_F);
+ break;
+ }
+ else if(e->faces(1)->g == &CLASS_EXTERIOR){
+ recur_tag(e->faces(0), &CLASS_F);
break;
}
- ++itt;
}
+ ++ite;
}
-
- // now find an edge that has belongs to one of the exterior
- // triangles
- {
- std::list<BDS_Edge*>::iterator ite = m->edges.begin();
- while (ite != m->edges.end()){
- BDS_Edge *e = *ite;
- if(e->g && e->numfaces() == 2){
- if(e->faces(0)->g == &CLASS_EXTERIOR){
- recur_tag(e->faces(1), &CLASS_F);
- break;
- }
- else if(e->faces(1)->g == &CLASS_EXTERIOR){
- recur_tag(e->faces(0), &CLASS_F);
- break;
- }
- }
- ++ite;
- }
- std::list<BDS_Face*>::iterator itt = m->triangles.begin();
- while (itt != m->triangles.end()){
- if ((*itt)->g == &CLASS_EXTERIOR) (*itt)->g = 0;
- ++itt;
- }
+ std::list<BDS_Face*>::iterator itt = m->triangles.begin();
+ while (itt != m->triangles.end()){
+ if ((*itt)->g == &CLASS_EXTERIOR) (*itt)->g = 0;
+ ++itt;
}
-
- {
- std::list<BDS_Edge*>::iterator ite = m->edges.begin();
- while (ite != m->edges.end()){
- BDS_Edge *e = *ite;
- if(e->g && e->numfaces() == 2){
- BDS_Point *oface[2];
- e->oppositeof(oface);
- if(oface[0]->iD < 0){
- recur_tag(e->faces(1), &CLASS_F);
- break;
- }
- else if(oface[1]->iD < 0){
- recur_tag(e->faces(0), &CLASS_F);
- break;
- }
- }
- ++ite;
+ }
+
+ {
+ std::list<BDS_Edge*>::iterator ite = m->edges.begin();
+ while (ite != m->edges.end()){
+ BDS_Edge *e = *ite;
+ if(e->g && e->numfaces() == 2){
+ BDS_Point *oface[2];
+ e->oppositeof(oface);
+ if(oface[0]->iD < 0){
+ recur_tag(e->faces(1), &CLASS_F);
+ break;
+ }
+ else if(oface[1]->iD < 0){
+ recur_tag(e->faces(0), &CLASS_F);
+ break;
+ }
}
- }
-
- ite = emb_edges.begin();
- while(ite != emb_edges.end()){
- if(!(*ite)->isMeshDegenerated())
- recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
++ite;
}
-
- // compute characteristic lengths at vertices
- if (!onlyInitialMesh){
+ }
+
+ ite = emb_edges.begin();
+ while(ite != emb_edges.end()){
+ if(!(*ite)->isMeshDegenerated())
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
+ ++ite;
+ }
+
+ // compute characteristic lengths at vertices
+ if (!onlyInitialMesh){
Msg::Debug("Computing mesh size field at mesh vertices %d",
edgesToRecover.size());
- for(int i = 0; i < doc.numPoints; i++){
- BDS_Point *pp = (BDS_Point*)doc.points[i].data;
+ std::set<BDS_Point*, PointLessThan>::iterator it = m->points.begin();
+ for(; it != m->points.end();++it){
+ // for(int i = 0; i < doc.numPoints; i++){
+ // BDS_Point *pp = (BDS_Point*)doc.points[i].data;
+ BDS_Point *pp = *it;
std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itv = recoverMap.find(pp);
if(itv != recoverMap.end()){
MVertex *here = itv->second;
@@ -1306,9 +1506,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
pp->lc() = pp->lcBGM();
}
}
- }
}
-
+
// delete useless stuff
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end()){
@@ -1317,7 +1516,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
++itt;
}
m->cleanup();
-
+
{
std::list<BDS_Edge*>::iterator ite = m->edges.begin();
while (ite != m->edges.end()){
@@ -1372,18 +1571,6 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
}
}
- if (Msg::GetVerbosity() == 10){
- GEdge *ge = new discreteEdge(gf->model(), 1000, 0, 0);
- MElementOctree octree(gf->model());
- Msg::Info("Writing voronoi and skeleton.pos");
- doc.Voronoi();
- doc.makePosView("voronoi.pos", gf);
- doc.printMedialAxis(octree.getInternalOctree(), "skeleton.pos", gf, ge);
- //todo add corners with lines with closest point
- ge->addPhysicalEntity(1000);
- gf->model()->add(ge);
- }
-
{
int nb_swap;
Msg::Debug("Delaunizing the initial mesh");
@@ -1487,9 +1674,11 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// delete the mesh
delete m;
- if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
- !CTX::instance()->mesh.optimizeLloyd && !onlyInitialMesh)
- recombineIntoQuads(gf);
+ if (1){
+ if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
+ !CTX::instance()->mesh.optimizeLloyd && !onlyInitialMesh)
+ recombineIntoQuads(gf);
+ }
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
@@ -1861,6 +2050,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// Use a divide & conquer type algorithm to create a triangulation.
// We add to the triangulation a box with 4 points that encloses the
// domain.
+#if 1 //OLD_CODE
{
DocRecord doc(nbPointsTotal + 4);
int count = 0;
@@ -1934,7 +2124,67 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
m->add_triangle(p1->iD, p2->iD, p3->iD);
}
}
-
+#else
+ {
+ /// FIXME FOR PERIODIC : SOME MVERTices SHOULD BE DUPLICATED ...
+ /// Still to be done...
+ printf("coucou1\n");
+ std::vector<MVertex*> v;
+ std::map<MVertex*, BDS_Point*> recoverMapInv;
+ 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_Point *pp = edgeLoop_BDS[j];
+ v.push_back(recoverMap[pp]);
+ recoverMapInv[recoverMap[pp]] = pp;
+ }
+ }
+
+ printf("coucou2 %d verices\n",v.size());
+ std::map<MVertex*,SPoint3> pos;
+ for(unsigned int i = 0; i < v.size(); i++) {
+ MVertex *v0 = v[i];
+ BDS_Point *p0 = recoverMapInv[v0];
+ pos[v0] = SPoint3(v0->x(),v0->y(),v0->z());
+ v0->setXYZ(p0->u,p0->v,0.0);
+ }
+ printf("coucou3\n");
+ std::vector<MTriangle*> result;
+ delaunayMeshIn2D(v, result, 0);
+ printf("coucou4\n");
+ // delaunayMeshIn2D(v, result, 0, & medgesToRecover);
+
+ for(unsigned int i = 0; i < v.size()-4; i++) {
+ MVertex *v0 = v[i];
+ SPoint3 pp = pos[v0];
+ v0->setXYZ(pp.x(),pp.y(),pp.z());
+ }
+ printf("coucou5\n");
+
+ // add the four corners
+ for(int ip = 0; ip < 4; ip++){
+ MVertex *vv = v[v.size()-ip-1];
+ BDS_Point *pp = m->add_point(-ip - 1, vv->x(),vv->y(), gf);
+ m->add_geom(gf->tag(), 2);
+ recoverMapInv[vv] = pp;
+ BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
+ pp->g = g;
+ }
+ printf("coucou6\n");
+ // add the triangles
+ for(unsigned int i = 0; i < result.size(); i++) {
+ MVertex *v0 = result[i]->getVertex(0);
+ MVertex *v1 = result[i]->getVertex(1);
+ MVertex *v2 = result[i]->getVertex(2);
+ BDS_Point *p0 = recoverMapInv[v0];
+ BDS_Point *p1 = recoverMapInv[v1];
+ BDS_Point *p2 = recoverMapInv[v2];
+ m->add_triangle(p0->iD, p1->iD, p2->iD);
+ }
+ printf("coucou7\n");
+ }
+#endif
+
// Recover the boundary edges and compute characteristic lenghts
// using mesh edge spacing
BDS_GeomEntity CLASS_F(1, 2);
@@ -2213,9 +2463,11 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// delete the mesh
delete m;
- if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
- !CTX::instance()->mesh.optimizeLloyd)
- recombineIntoQuads(gf);
+ if (1){
+ if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
+ !CTX::instance()->mesh.optimizeLloyd)
+ recombineIntoQuads(gf);
+ }
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
@@ -2309,6 +2561,10 @@ void meshGFace::operator() (GFace *gf, bool print)
return;
}
+#ifdef OLD_CODE
+ quadMeshRemoveHalfOfOneDMesh halfmesh (gf);
+#endif
+
if ((gf->getNativeType() != GEntity::AcisModel ||
(!gf->periodic(0) && !gf->periodic(1))) &&
(noSeam(gf) || gf->getNativeType() == GEntity::GmshModel ||
@@ -2325,25 +2581,9 @@ void meshGFace::operator() (GFace *gf, bool print)
Msg::Debug("Type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
- // do the 2D mesh in several passes. For second and other passes,
- // a background mesh is constructed with the previous mesh and
- // nodal values of the metric are computed that take into account
- // complex size fields that are tedious to evaluate on the fly
- if (!twoPassesMesh)return;
- twoPassesMesh--;
- if (backgroundMesh::current()){
- backgroundMesh::unset();
- //backgroundMesh::set(gf);
- }
- if (CTX::instance()->mesh.saveAll){
- backgroundMesh::set(gf);
- char name[256];
- sprintf(name,"bgm-%d.pos",gf->tag());
- backgroundMesh::current()->print(name,gf);
- sprintf(name,"cross-%d.pos",gf->tag());
- backgroundMesh::current()->print(name,gf,1);
- }
- (*this)(gf);
+#ifdef OLD_CODE
+ halfmesh.finish();
+#endif
}
bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index d1a2c2ae51..857e24e4a4 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -21,11 +21,10 @@ class GFaceCompound;
// Create the mesh of the face
class meshGFace {
const bool repairSelfIntersecting1dMesh;
- int twoPassesMesh;
bool onlyInitialMesh;
public :
- meshGFace(bool r = true, int t = 0)
- : repairSelfIntersecting1dMesh(r), twoPassesMesh(t), onlyInitialMesh(false)
+ meshGFace(bool r = true)
+ : repairSelfIntersecting1dMesh(r), onlyInitialMesh(false)
{
}
void operator()(GFace *, bool print=true);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 63f25b52c5..93e23b3bc6 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -24,6 +24,7 @@
#include "GFaceCompound.h"
#include "intersectCurveSurface.h"
#include "surfaceFiller.h"
+#include "HilbertCurve.h"
double LIMIT_ = 0.5 * sqrt(2.0) * 1;
int N_GLOBAL_SEARCH;
@@ -44,29 +45,29 @@ static bool isBoundary(MTri3 *t, double limit_, int &active)
*/
template <class ITERATOR>
-void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData & data, bool param=true)
+void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData * data)
{
FILE *ff = Fopen (name,"w");
fprintf(ff,"View\"test\"{\n");
while ( it != end ){
MTri3 *worst = *it;
if (!worst->isDeleted()){
- if (param)
+ if (data)
fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- data.Us[data.getIndex((worst)->tri()->getVertex(0))],
- data.Vs[data.getIndex((worst)->tri()->getVertex(0))],
+ data->Us[data->getIndex((worst)->tri()->getVertex(0))],
+ data->Vs[data->getIndex((worst)->tri()->getVertex(0))],
0.0,
- data.Us[data.getIndex((worst)->tri()->getVertex(1))],
- data.Vs[data.getIndex((worst)->tri()->getVertex(1))],
+ data->Us[data->getIndex((worst)->tri()->getVertex(1))],
+ data->Vs[data->getIndex((worst)->tri()->getVertex(1))],
0.0,
- data.Us[data.getIndex((worst)->tri()->getVertex(2))],
- data.Vs[data.getIndex((worst)->tri()->getVertex(2))],
+ data->Us[data->getIndex((worst)->tri()->getVertex(2))],
+ data->Vs[data->getIndex((worst)->tri()->getVertex(2))],
0.0,
(worst)->getRadius(),
(worst)->getRadius(),
(worst)->getRadius());
else
- fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
+ fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%d,%d,%d};\n",
(worst)->tri()->getVertex(0)->x(),
(worst)->tri()->getVertex(0)->y(),
(worst)->tri()->getVertex(0)->z(),
@@ -76,9 +77,9 @@ void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData & data, bo
(worst)->tri()->getVertex(2)->x(),
(worst)->tri()->getVertex(2)->y(),
(worst)->tri()->getVertex(2)->z(),
- (worst)->getRadius(),
- (worst)->getRadius(),
- (worst)->getRadius()
+ (worst)->tri()->getVertex(0)->getNum(),
+ (worst)->tri()->getVertex(1)->getNum(),
+ (worst)->tri()->getVertex(2)->getNum()
);
}
++it;
@@ -411,65 +412,47 @@ int inCircumCircle(MTriangle *base,
}
template <class ITER>
-void connectTris(ITER beg, ITER end)
+void connectTris(ITER beg, ITER end, std::vector<edgeXface> &conn)
{
- std::set<edgeXface> conn;
+ conn.clear();
while (beg != end){
if (!(*beg)->isDeleted()){
- for (int i = 0; i < 3; i++){
- edgeXface fxt(*beg, i);
- std::set<edgeXface>::iterator found = conn.find(fxt);
- if (found == conn.end())
- conn.insert(fxt);
- else if (found->t1 != *beg){
- found->t1->setNeigh(found->i1, *beg);
- (*beg)->setNeigh(i, found->t1);
- }
+ for (int j = 0; j < 3; j++){
+ conn.push_back(edgeXface(*beg, j));
}
}
++beg;
}
-}
-
-template <class ITER>
-void connectTris_vector(ITER beg, ITER end)
-{
- std::vector<edgeXface> conn;
- // std::set<edgeXface> conn;
- while (beg != end){
- if (!(*beg)->isDeleted()){
- for (int i = 0; i < 3; i++){
- edgeXface fxt(*beg, i);
- // std::set<edgeXface>::iterator found = conn.find(fxt);
- std::vector<edgeXface>::iterator found = std::find(conn.begin(), conn.end(), fxt);
- if (found == conn.end())
- conn.push_back(fxt);
- // conn.insert(fxt);
- else if (found->t1 != *beg){
- found->t1->setNeigh(found->i1, *beg);
- (*beg)->setNeigh(i, found->t1);
- }
- }
+ if (!conn.size())return;
+ std::sort(conn.begin(), conn.end());
+
+ for (unsigned int i=0;i<conn.size()-1;i++){
+ edgeXface &f1 = conn[i];
+ edgeXface &f2 = conn[i+1];
+ if (f1 == f2 && f1.t1 != f2.t1){
+ f1.t1->setNeigh(f1.i1, f2.t1);
+ f2.t1->setNeigh(f2.i1, f1.t1);
+ ++i;
}
- ++beg;
}
}
-
-
void connectTriangles(std::list<MTri3*> &l)
{
- connectTris(l.begin(), l.end());
+ std::vector<edgeXface> conn;
+ connectTris(l.begin(), l.end(),conn);
}
void connectTriangles(std::vector<MTri3*> &l)
{
- connectTris(l.begin(), l.end());
+ std::vector<edgeXface> conn;
+ connectTris(l.begin(), l.end(),conn);
}
void connectTriangles(std::set<MTri3*, compareTri3Ptr> &l)
{
- connectTris(l.begin(), l.end());
+ std::vector<edgeXface> conn;
+ connectTris(l.begin(), l.end(),conn);
}
int inCircumCircleTangentPlane(MTriangle *t,
@@ -490,6 +473,20 @@ int inCircumCircleTangentPlane(MTriangle *t,
return (result > 0) ? 1 : 0;
}
+int inCircumCircleXY(MTriangle *t, MVertex *v)
+{
+ MVertex *v1 = t->getVertex(0);
+ MVertex *v2 = t->getVertex(1);
+ MVertex *v3 = t->getVertex(2);
+ double p1[2] = {v1->x(),v1->y()};
+ double p2[2] = {v2->x(),v2->y()};
+ double p3[2] = {v3->x(),v3->y()};
+ double pp[2] = {v->x(),v->y()};
+ double result = robustPredicates::incircle(p1, p2, p3, pp) *
+ robustPredicates::orient2d(p1, p2, p3);
+ return (result > 0) ? 1 : 0;
+}
+
void recurFindCavityTangentPlane(std::list<edgeXface> &shell,
std::list<MTri3*> &cavity,
@@ -567,6 +564,27 @@ bool findCavityTangentPlane(GFace *gf, double *center,
}
+void recurFindCavity(std::vector<edgeXface> &shell, std::vector<MTri3*> &cavity,
+ MVertex *v, MTri3 *t)
+{
+ t->setDeleted(true);
+ // the cavity that has to be removed because it violates delaunay
+ // criterion
+ cavity.push_back(t);
+
+ for (int i = 0; i < 3; i++){
+ MTri3 *neigh = t->getNeigh(i) ;
+ if (!neigh)
+ shell.push_back(edgeXface(t, i));
+ else if (!neigh->isDeleted()){
+ int circ = inCircumCircleXY(neigh->tri(), v);
+ if (circ)
+ recurFindCavity(shell, cavity, v, neigh);
+ else
+ shell.push_back(edgeXface(t, i));
+ }
+ }
+}
void recurFindCavity(std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
double *v, double *param, MTri3 *t, bidimMeshData & data)
@@ -696,6 +714,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
if (shell.size() != cavity.size() + 2) return false;
std::list<MTri3*> new_cavity;
+ std::vector<edgeXface> conn;
// check that volume is conserved
double newVolume = 0;
@@ -756,7 +775,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && !onePointIsTooClose){
- connectTris_vector(new_cavity.begin(), new_cavity.end());
+ connectTris(new_cavity.begin(), new_cavity.end(),conn);
// printf("%d %d\n",shell.size(),cavity.size());
// clock_t t1 = clock();
// 30 % of the time is spent here !!!
@@ -828,6 +847,59 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
}
+bool invMapXY (MTriangle *t, MVertex *v){
+ MVertex *v0 = t->getVertex(0);
+ MVertex *v1 = t->getVertex(1);
+ MVertex *v2 = t->getVertex(2);
+ double mat[2][2],b[2],uv[2],tol=1.e-6;
+ mat[0][0] = v1->x() - v0->x();
+ mat[0][1] = v2->x() - v0->x();
+ mat[1][0] = v1->y() - v0->y();
+ mat[1][1] = v2->y() - v0->y();
+
+ b[0] = v->x() - v0->x();
+ b[1] = v->y() - v0->y();
+ sys2x2(mat, b, uv);
+
+ if(uv[0] >= -tol &&
+ uv[1] >= -tol &&
+ uv[0] <= 1. + tol &&
+ uv[1] <= 1. + tol &&
+ 1. - uv[0] - uv[1] > -tol) {
+ return true;
+ }
+ return false;
+
+}
+
+static MTri3* search4Triangle (MTri3 *t, MVertex *v, int maxx, int &ITER) {
+ bool inside = invMapXY(t->tri(), v);
+ SPoint3 q1 (v->x(),v->y(),0);
+ if (inside) return t;
+ while (1){
+ SPoint3 q2 = t->tri()->barycenter();
+ int i;
+ for (i = 0; i < 3; i++){
+ int i1 = i == 0 ? 2 : i-1;
+ int i2 = i;
+ MVertex *v1 = t->tri()->getVertex(i1);
+ MVertex *v2 = t->tri()->getVertex(i2);
+ SPoint3 p1 (v1->x(),v1->y(),0);
+ SPoint3 p2 (v2->x(),v2->y(),0);
+ double xcc[2];
+ if (intersection_segments(p1, p2, q1, q2, xcc)) break;
+ }
+ if (i >= 3) break;
+ t = t->getNeigh(i);
+ if (!t) break;
+ bool inside = invMapXY(t->tri(), v);
+ // printf("ITER %d %d\n",ITER,inside);
+ if (inside) return t;
+ if (ITER++ > .5*maxx) break;
+ }
+ return 0;
+}
+
static MTri3* search4Triangle (MTri3 *t, double pt[2], bidimMeshData & data,
std::set<MTri3*,compareTri3Ptr> &AllTris, double uv[2], bool force = false) {
@@ -1306,9 +1378,9 @@ void bowyerWatsonFrontal(GFace *gf,
if(ITERATION % 10== 0 && CTX::instance()->mesh.saveAll){
char name[245];
sprintf(name,"delFrontal_GFace_%d_Layer_%d.pos",gf->tag(),ITERATION);
- _printTris (name, AllTris.begin(), AllTris.end(), DATA,true);
+ _printTris (name, AllTris.begin(), AllTris.end(), &DATA);
sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
- _printTris (name, ActiveTris.begin(), ActiveTris.end(), DATA,true);
+ _printTris (name, ActiveTris.begin(), ActiveTris.end(), &DATA);
}
/* if(ITER % 100== 0){
char name[245];
@@ -1729,3 +1801,273 @@ void bowyerWatsonParallelograms(GFace *gf,
}
#endif
}
+
+
+static void initialSquare(std::vector<MVertex*> &v,
+ MVertex *box[4],
+ std::vector<MTri3*> &t){
+ SBoundingBox3d bbox ;
+ for (size_t i=0;i<v.size();i++){
+ MVertex *pv = v[i];
+ bbox += SPoint3(pv->x(),pv->y(),pv->z());
+ }
+ bbox *= 1.3;
+ box[0] = new MVertex (bbox.min().x(),bbox.min().y(),0);
+ box[1] = new MVertex (bbox.max().x(),bbox.min().y(),0);
+ box[2] = new MVertex (bbox.max().x(),bbox.max().y(),0);
+ box[3] = new MVertex (bbox.min().x(),bbox.max().y(),0);
+ std::vector<MTriangle*> t_box;
+ MTriangle *t0 = new MTriangle (box[0],box[1],box[2]);
+ MTriangle *t1 = new MTriangle (box[2],box[3],box[0]);
+ t.push_back(new MTri3(t0,0.0));
+ t.push_back(new MTri3(t1,0.0));
+ connectTriangles(t);
+}
+
+
+MTri3 * getTriToBreak (MVertex *v, std::vector<MTri3*> &t, int &NB_GLOBAL_SEARCH, int &ITER){
+ // last inserted is used as starting point
+ // we know it is not deleted
+ unsigned int k = t.size() - 1;
+ while(t[k]->isDeleted()){
+ k--;
+ }
+ MTri3 *start = t[k];
+ start = search4Triangle (start,v,(int)t.size(),ITER);
+ if (start)return start;
+ // printf("Global Search has to be done\n");
+ NB_GLOBAL_SEARCH++;
+ for (size_t i = 0;i<t.size();i++){
+ if (!t[i]->isDeleted() && inCircumCircleXY(t[i]->tri(),v))return t[i];
+ }
+ return 0;
+}
+
+bool triOnBox (MTriangle *t, MVertex *box[4]){
+ for (size_t i = 0;i<3;i++)
+ for (size_t j = 0;j<4;j++)
+ if (t->getVertex(i) == box[j])return true;
+ return false;
+}
+
+// vertices are supposed to be sitting in the XY plane !
+
+void recoverEdges (std::vector<MTri3*> &t, std::vector<MEdge> &edges);
+
+void delaunayMeshIn2D(std::vector<MVertex*> &v,
+ std::vector<MTriangle*> &result,
+ bool removeBox,
+ std::vector<MEdge> *edgesToRecover)
+{
+ std::vector<MTri3*> t;
+ t.reserve (v.size()*2);
+ std::vector<edgeXface> conn;
+ std::vector<edgeXface> shell;
+ std::vector<MTri3*> cavity;
+ MVertex *box[4];
+ initialSquare (v,box,t);
+
+ int NB_GLOBAL_SEARCH = 0;
+ int AVG_ITER = 0;
+ SortHilbert(v);
+ double t1 = Cpu();
+
+ double ta=0,tb=0,tc=0,td=0,T;
+ for (size_t i=0;i<v.size();i++){
+ MVertex *pv = v[i];
+
+ int NITER = 0;
+ T = Cpu();
+ MTri3 * found = getTriToBreak (pv,t,NB_GLOBAL_SEARCH,NITER);
+ ta += Cpu()-T;
+ AVG_ITER += NITER;
+ if(!found) {
+ Msg::Error("cannot insert a point in 3D Delaunay");
+ continue;
+ }
+ shell.clear();
+ cavity.clear();
+
+ T = Cpu();
+ recurFindCavity(shell, cavity, pv, found);
+ tb += Cpu()-T;
+ double V = 0.0;
+ for (unsigned int k=0;k<cavity.size();k++)V+=fabs(cavity[k]->tri()->getVolume());
+
+ std::vector<MTri3*> extended_cavity;
+ double Vb = 0.0;
+
+ T = Cpu();
+ for (unsigned int count = 0; count < shell.size(); count++){
+ const edgeXface &fxt = shell[count];
+ MTriangle *tr;
+ MTri3 *t3;
+ MVertex *v0 = fxt.v[0];
+ MVertex *v1 = fxt.v[1];
+ MTri3 *otherSide = fxt.t1->getNeigh(fxt.i1);
+ if (count < cavity.size()){
+ t3 = cavity[count];
+ tr = t3->tri() ;
+ tr->setVertex(0,v0);
+ tr->setVertex(1,v1);
+ tr->setVertex(2,pv);
+ }
+ else{
+ tr = new MTriangle(v0,v1,pv);
+ t3 = new MTri3(tr, 0.0);
+ t.push_back(t3);
+ }
+ Vb+= fabs(tr->getVolume());
+ extended_cavity.push_back(t3);
+ if (otherSide)
+ extended_cavity.push_back(otherSide);
+ }
+ tc += Cpu()-T;
+ if (fabs(Vb-V) > 1.e-8 * (Vb+V))printf("%12.5E %12.5E\n",Vb,V);
+
+ for (unsigned int k=0;k<std::min(cavity.size(),shell.size());k++){
+ cavity[k]->setDeleted(false);
+ for (unsigned int l=0;l<3;l++){
+ cavity[k]->setNeigh(l,0);
+ }
+ }
+ T = Cpu();
+ connectTris(extended_cavity.begin(),extended_cavity.end(),conn);
+ td += Cpu()-T;
+ }
+
+ double t2 = Cpu();
+ Msg::Debug("Delaunay 2D done for %d points : CPU = %g, %d global searches, AVG walk size %g",v.size(), t2-t1,NB_GLOBAL_SEARCH,1.+(double)AVG_ITER/v.size());
+ // printf("%g %g %g %g --> %g(%g)\n",ta,tb,tc,td,t2-t1,ta+tb+tc+td);
+
+ if (edgesToRecover)recoverEdges (t,*edgesToRecover);
+
+ // FILE *f = fopen ("tri.pos","w");
+ // fprintf(f,"View \"\"{\n");
+ for (size_t i = 0;i<t.size();i++){
+ if (t[i]->isDeleted() || (removeBox && triOnBox (t[i]->tri(),box))) delete t[i]->tri();
+ else {
+ result.push_back(t[i]->tri());
+ // t[i]->tri()->writePOS (f, false,false,true,false,false,false);
+ }
+ delete t[i];
+ }
+
+ if (removeBox)for (int i=0;i<4;i++)delete box[i];
+ else for (int i=0;i<4;i++)v.push_back(box[i]);
+
+ // fprintf(f,"};\n");
+ // fclose(f);
+}
+
+bool swapedge (MVertex *v1 ,MVertex *v2, MVertex *v3, MVertex *v4, MTri3* t1, int iLocalEdge){
+ MTri3 *t2 = t1->getNeigh(iLocalEdge);
+ if(!t2) return false;
+
+ MTriangle *t1b = new MTriangle(v2, v3, v4);
+ MTriangle *t2b = new MTriangle(v4, v3, v1);
+ double BEFORE = t1->tri()->getVolume()+t2->tri()->getVolume();
+ double AFTER = t1b->getVolume()+t2b->getVolume();
+ // printf("swapping %d %d %d %d %g %g\n",v1->getNum(),v2->getNum(),v3->getNum(),v4->getNum(),BEFORE,AFTER);
+ if (fabs(BEFORE-AFTER)/BEFORE > 1.e-8){
+ delete t1b;
+ delete t2b;
+ return false;
+ }
+ // printf("volumes ok\n");
+
+ delete t1->tri();
+ delete t2->tri();
+ t1->setTri(t1b);
+ t2->setTri(t2b);
+
+ std::set<MTri3*> cavity;
+ cavity.insert(t1);
+ cavity.insert(t2);
+ for(int i = 0; i < 3; i++){
+ if(t1->getNeigh(i))
+ cavity.insert(t1->getNeigh(i));
+ if(t2->getNeigh(i))
+ cavity.insert(t2->getNeigh(i));
+ }
+ std::vector<edgeXface> conn;
+ connectTris(cavity.begin(), cavity.end(), conn);
+ return true;
+}
+
+bool diffend (MVertex *v1, MVertex *v2, MVertex *p1, MVertex *p2){
+ if (v1 == p1 || v2 == p1 || v1 == p2 || v2 == p2)return false;
+ return true;
+}
+
+/*
+
+*/
+
+static bool recoverEdgeBySwaps (std::vector<MTri3*> &t, MVertex *mv1, MVertex *mv2, std::vector<MEdge> &edges){
+
+ SPoint3 pv1 (mv1->x(),mv1->y(),0);
+ SPoint3 pv2 (mv2->x(),mv2->y(),0);
+ double xcc[2];
+ for (unsigned int i=0;i<t.size();i++){
+ for (unsigned int j=0;j<3;j++){
+ MVertex *v1 = t[i]->tri()->getVertex((j+2)%3);
+ MVertex *v2 = t[i]->tri()->getVertex(j);
+ MVertex *v3 = t[i]->tri()->getVertex((j+1)%3);
+ MVertex *o = t[i]->otherSide(j);
+ if (o){
+ SPoint3 p1 (v1->x(),v1->y(),0);
+ SPoint3 p2 (v2->x(),v2->y(),0);
+ SPoint3 p3 (v3->x(),v3->y(),0);
+ SPoint3 po (o->x(),o->y(),0);
+ if (diffend(v1,v2,mv1,mv2)){
+ if (intersection_segments (p1, p2, pv1, pv2,xcc)){
+ // if (std::binary_search(edges.begin(),edges.end(),MEdge(v1,v2),Less_Edge)){
+ // Msg::Error("1D mesh self intersects");
+ // }
+ if (!intersection_segments(po, p3, pv1, pv2,xcc) || (v3 == mv1 || o == mv1 || v3 == mv2 || o == mv2)){
+ if(swapedge (v1,v2,v3,o,t[i],j))return true;
+ }
+ }
+ }
+ }
+ }
+ }
+ return false;
+}
+
+// recover the edges by edge swapping in the triangulation.
+// edges are not supposed to
+
+void recoverEdges (std::vector<MTri3*> &t, std::vector<MEdge> &edges)
+{
+ Less_Edge le;
+ std::sort(edges.begin(),edges.end(),le);
+ std::set<MEdge,Less_Edge> setOfMeshEdges;
+ for (size_t i = 0;i<t.size();i++){
+ for (int j=0;j<3;j++){
+ setOfMeshEdges.insert(t[i]->tri()->getEdge(j));
+ }
+ }
+
+ std::vector<MEdge> edgesToRecover;
+ for (unsigned int i=0;i<edges.size();i++){
+ if (setOfMeshEdges.find(edges[i])==setOfMeshEdges.end())
+ edgesToRecover.push_back(edges[i]);
+ }
+
+ Msg::Info("%d edges to recover among %d edges",edgesToRecover.size(),edges.size());
+ // int iter = 0;
+ // char name[256];
+ // sprintf(name,"iter%d.pos",iter++);
+ // _printTris (name, t.begin(),t.end(),0);
+ for (unsigned int i=0;i<edgesToRecover.size();i++){
+ MVertex *mstart = edgesToRecover[i].getVertex(0);
+ MVertex *mend = edgesToRecover[i].getVertex(1);
+ Msg::Info("recovering edge %d %d",mstart->getNum(),mend->getNum());
+ int iter;
+ while(recoverEdgeBySwaps (t, mstart, mend,edges)) {
+ iter ++;
+ }
+ }
+}
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 57f175b886..d2af45bf38 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -90,8 +90,17 @@ class MTri3
bool isDeleted() const { return deleted; }
void forceRadius(double r) { circum_radius = r; }
inline double getRadius() const { return circum_radius; }
-
+ inline MVertex *otherSide (int i){
+ MTri3 *n = neigh[i];
+ if (!n)return 0;
+ MVertex *v1 = base->getVertex((i+2)%3);
+ MVertex *v2 = base->getVertex(i);
+ for (int j=0;j<3;j++)
+ if (n->tri()->getVertex(j) != v1 && n->tri()->getVertex(j) != v2)return n->tri()->getVertex(j);
+ return 0;
+ }
MTri3(MTriangle *t, double lc, SMetric3 *m = 0, bidimMeshData * data = 0, GFace *gf = 0);
+ inline void setTri(MTriangle *t) { base = t; }
inline MTriangle *tri() const { return base; }
inline void setNeigh(int iN , MTri3 *n) { neigh[iN] = n; }
inline MTri3 *getNeigh(int iN ) const { return neigh[iN]; }
@@ -152,6 +161,10 @@ void buildBackGroundMesh (GFace *gf,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
+void delaunayMeshIn2D(std::vector<MVertex*> &,
+ std::vector<MTriangle*> &, bool removeBox = true,
+ std::vector<MEdge> *edgesToRecover = 0);
+
struct edgeXface
{
MVertex *v[2];
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e4b53917a7..a91f8b73d4 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1658,6 +1658,7 @@ static int _defectsRemovalBunin(GFace *gf, int maxCavitySize)
{
if (maxCavitySize == 0)return 0;
+ printf("ARGH\n");
v2t_cont adj;
std::vector<MElement*> c;
buildVertexToElement(gf->quadrangles, adj);
@@ -3294,6 +3295,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
{
// never recombine a face that is part of a compound!
if(gf->getCompound()) return 0;
+ if(gf->triangles.size() == 0) return 1;
int success = 1;
@@ -3586,15 +3588,16 @@ static double printStats(GFace *gf,const char *message)
Qav += Q;
Qmin = std::min(Q,Qmin);
}
- Msg::Info("%s : %5d quads %4d invalid quads %4d quads with Q < 0.1 "
- "Avg Q = %12.5E Min %12.5E", message, gf->quadrangles.size(),
+ Msg::Info("%s : %5d quads %5d triangles %4d invalid quads %4d quads with Q < 0.1 "
+ "Avg Q = %12.5E Min %12.5E", message, gf->quadrangles.size(), gf->triangles.size(),
nbInv, nbBad, Qav/gf->quadrangles.size(), Qmin);
return Qmin;
}
void recombineIntoQuads(GFace *gf,
bool topologicalOpti,
- bool nodeRepositioning)
+ bool nodeRepositioning,
+ double minqual)
{
double t1 = Cpu();
@@ -3608,7 +3611,7 @@ void recombineIntoQuads(GFace *gf,
// removeFourTrianglesNodes(gf, false);
if (saveAll) gf->model()->writeMSH("before.msh");
- int success = _recombineIntoQuads(gf, 0);
+ int success = _recombineIntoQuads(gf, minqual);
if (saveAll) gf->model()->writeMSH("raw.msh");
printStats (gf, "BEFORE OPTIMIZATION");
@@ -3621,30 +3624,27 @@ void recombineIntoQuads(GFace *gf,
if(haveParam){
if (saveAll) gf->model()->writeMSH("smoothed.msh");
int ITER=0;
- int ITERB = 0;
- int optistatus[6] = {0,0,0,0,0,0};
+ // int optistatus[6] = {0,0,0,0,0,0};
std::set<MEdge,Less_Edge> prioritory;
+ double exbad = -100;
while(1){
- int maxCavitySize = CTX::instance()->mesh.bunin;
- optistatus[0] = (ITERB == 1) ?splitFlatQuads(gf, .01, prioritory) : 0;
- optistatus[1] = removeTwoQuadsNodes(gf);
- optistatus[4] = _defectsRemovalBunin(gf,maxCavitySize);
- optistatus[2] = removeDiamonds(gf) ;
- laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
- optistatus[3] = edgeSwapQuadsForBetterQuality(gf,.01, prioritory);
- optistatus[5] = optiSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
- optistatus[5] = (ITERB == 1) ?
- untangleInvalidQuads(gf, CTX::instance()->mesh.nbSmoothing) : 0;
+ // int maxCavitySize = CTX::instance()->mesh.bunin;
+ // optistatus[0] = (ITERB == 1) ?splitFlatQuads(gf, .01, prioritory) : 0;
+ //optistatus[1] =
+ removeTwoQuadsNodes(gf);
+ //optistatus[4] = _defectsRemovalBunin(gf,36);
+ //optistatus[2] =
+ removeDiamonds(gf) ;
+ if(haveParam)laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
+ //optistatus[3] =
+ edgeSwapQuadsForBetterQuality(gf,minqual, prioritory);
+ //optistatus[5] =
+ optiSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
+ // optistatus[5] = untangleInvalidQuads(gf, CTX::instance()->mesh.nbSmoothing);
double bad = printStats(gf, "IN OPTIMIZATION");
- if (bad > .1) break;
- if (ITER == 10){
- ITERB = 1;
- }
+ if (bad > minqual || exbad == bad) break;
+ exbad = bad;
if (ITER > 20) break;
- int nb = 0;
- for (int i=0;i<6;i++) nb += optistatus[i];
- if (!nb && ITERB == 0) ITERB = 1;
- else if (!nb) break;
ITER ++;
}
}
@@ -3656,10 +3656,9 @@ void recombineIntoQuads(GFace *gf,
}
double t2 = Cpu();
Msg::Info("Blossom recombination algorithm completed (%g s)", t2 - t1);
- quadsToTriangles(gf, .01);
+ quadsToTriangles(gf, minqual);
if(haveParam) {
laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
- // optiSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
}
// removeDiamonds(gf);
// removeTwoQuadsNodes(gf);
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index c5682c3e4c..6dde944f99 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -81,7 +81,8 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
void computeEquivalences(GFace *gf,bidimMeshData &DATA);
void recombineIntoQuads(GFace *gf,
bool topologicalOpti = true,
- bool nodeRepositioning = true);
+ bool nodeRepositioning = true,
+ double minqual = 0.1);
//used for meshGFaceRecombine development
int recombineWithBlossom(GFace *gf, double dx, double dy,
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 00b22d068c..29d6c584dc 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1820,7 +1820,7 @@ void sanityCheck1(MTet4 *t)
SortHilbert(v);
double t1 = Cpu();
- // double ta=0,tb=0,tc=0,td=0,T;
+ /// double ta=0,tb=0,tc=0,td=0,T;
for (size_t i=0;i<v.size();i++){
MVertex *pv = v[i];
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 9b26a8a5fb..4748ba728a 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1282,8 +1282,8 @@ void signedDistancesPointsEllipseLine(std::vector<double>&distances,
}
}
-int intersection_segments(SPoint3 &p1, SPoint3 &p2,
- SPoint3 &q1, SPoint3 &q2,
+int intersection_segments(const SPoint3 &p1, const SPoint3 &p2,
+ const SPoint3 &q1, const SPoint3 &q2,
double x[2])
{
double xp_max = std::max(p1.x(), p2.x());
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index a19affad2e..bb059023c3 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -142,8 +142,8 @@ void signedDistancesPointsEllipseLine (std::vector<double>&distances,
const std::vector<SPoint3> &pts,
const SPoint3 &p1, const SPoint3 &p2);
-int intersection_segments (SPoint3 &p1, SPoint3 &p2,
- SPoint3 &q1, SPoint3 &q2,
+int intersection_segments (const SPoint3 &p1, const SPoint3 &p2,
+ const SPoint3 &q1, const SPoint3 &q2,
double x[2]);
//tools for projection onto plane
diff --git a/benchmarks/2d/Square-01.geo b/benchmarks/2d/Square-01.geo
index 42006dad03..6a8d0aee8c 100644
--- a/benchmarks/2d/Square-01.geo
+++ b/benchmarks/2d/Square-01.geo
@@ -1,6 +1,6 @@
fact = 100;
lc = .1 * fact;
-Point(1) = {0.0,0.0,0,lc*1};
+Point(1) = {0.0,0.0,0,lc*2.5e-10};
Point(2) = {1* fact,0.0,0,lc*1};
Point(3) = {1* fact,1* fact,0,lc};
Point(4) = {0,1* fact,0,lc*1};
--
GitLab