From c03d82fbc9cf9de7b60697b497dfaf7bd51a0433 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 25 Oct 2013 12:33:53 +0000
Subject: [PATCH] Fixed surface meshing when singular parametrizations are
presents
TODO : avoid topological unconsistencies for (e.g.) cones when
delaunay meshing is involved.
---
Common/CommandLine.cpp | 2 -
Common/GmshDefines.h | 1 -
Geo/boundaryLayersData.cpp | 182 ++++++++++++++++++-
Geo/boundaryLayersData.h | 2 +-
Mesh/meshGFace.cpp | 6 -
Mesh/meshGFaceDelaunayInsertion.cpp | 262 ++++++++++++++++++++--------
Mesh/meshGFaceDelaunayInsertion.h | 3 -
Mesh/meshGFaceOptimize.cpp | 218 ++++++++++++++++-------
Mesh/meshGFaceOptimize.h | 1 +
benchmarks/2d/Square-01.geo | 6 +-
benchmarks/2d/conge.geo | 4 +-
11 files changed, 531 insertions(+), 156 deletions(-)
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 0cb5506c66..925e484edc 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -792,8 +792,6 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL_QUAD;
else if(!strncmp(argv[i], "pack", 4))
CTX::instance()->mesh.algo2d = ALGO_2D_PACK_PRLGRMS;
- else if(!strncmp(argv[i], "ruppert", 7))
- CTX::instance()->mesh.algo2d = ALGO_2D_RUPPERT;
else if(!strncmp(argv[i], "front2d", 7) || !strncmp(argv[i], "frontal", 7))
CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL;
else if(!strncmp(argv[i], "bamg",4))
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 6226ccd540..6d1cf9eff7 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -229,7 +229,6 @@
#define ALGO_2D_BAMG 7
#define ALGO_2D_FRONTAL_QUAD 8
#define ALGO_2D_PACK_PRLGRMS 9
-#define ALGO_2D_RUPPERT 10
// 3D meshing algorithms (numbers should not be changed)
#define ALGO_3D_DELAUNAY 1
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index 8b1bfd9a80..9adc9e9557 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -15,6 +15,11 @@
#include "MEdge.h"
#include "boundaryLayersData.h"
#include "OS.h"
+#include "BackgroundMesh.h"
+
+#if defined(HAVE_RTREE)
+#include "rtree.h"
+#endif
#if !defined(HAVE_MESH) || !defined(HAVE_ANN)
@@ -777,6 +782,8 @@ bool buildAdditionalPoints2D(GFace *gf)
}
// DEBUG STUFF
+ _columns->filterPoints(gf,0.21);
+
char name[256];
sprintf(name,"points_face_%d.pos",gf->tag());
FILE *f = Fopen (name,"w");
@@ -1087,9 +1094,9 @@ BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
}
}
}
- else
+ else {
createColumnsBetweenFaces (gr,*it,blf,_columns,_allGFaces,_faces,_normals,_treshold);
-
+ }
}
// DEBUG STUFF
@@ -1114,4 +1121,175 @@ BoundaryLayerColumns *buildAdditionalPoints3D(GRegion *gr)
return _columns;
}
+struct blPoint_wrapper
+{
+ bool _tooclose;
+ MVertex *_v;
+ std::map<MVertex*,MVertex*> &_v2v;
+ blPoint_wrapper (MVertex *v, std::map<MVertex*,MVertex*> &v2v)
+ : _tooclose(false), _v(v), _v2v(v2v) {}
+};
+
+struct blPoint_rtree
+{
+ MVertex *_v;
+ double _size;
+ blPoint_rtree (MVertex *v, double size) :
+ _v(v), _size(size) {}
+ bool inExclusionZone (MVertex *v){
+ double d = _v->distance(v);
+ //printf("d = %12.5E\n",d);
+ if (d <= _size) return true;
+ return false;
+ }
+ void minmax (double min[3], double max[3]){
+ min[0] = _v->x() - _size;
+ min[1] = _v->y() - _size;
+ min[2] = _v->z() - _size;
+ max[0] = _v->x() + _size;
+ max[1] = _v->y() + _size;
+ max[2] = _v->z() + _size;
+ }
+};
+
+
+bool rtree_callback(blPoint_rtree *neighbour,void* point){
+ blPoint_wrapper *w = static_cast<blPoint_wrapper*>(point);
+
+ const MVertex *from_1 = w->_v2v[neighbour->_v];
+ const MVertex *from_2 = w->_v2v[w->_v];
+
+ // printf("%p %p\n",from_1,from_2);
+
+ if (from_1 == from_2) {
+ return true;
+ }
+
+ if (neighbour->inExclusionZone(w->_v)){
+ w->_tooclose = true;
+ return false;
+ }
+ return true;
+}
+
+bool inExclusionZone_filter (MVertex* p,
+ std::map <MVertex*, MVertex*> &v2v,
+ RTree< blPoint_rtree *,double,3,double> &rtree){
+ // should assert that the point is inside the domain
+ {
+ double u, v;
+ p->getParameter(0,u);
+ p->getParameter(1,v);
+ if (!backgroundMesh::current()->inDomain(u,v,0)) return true;
+ }
+
+ blPoint_wrapper w (p,v2v);
+ double _min[3] = {p->x()-1.e-1, p->y()-1.e-1,p->z()-1.e-1};
+ double _max[3] = {p->x()+1.e-1, p->y()+1.e-1,p->z()+1.e-1};
+ rtree.Search(_min,_max,rtree_callback,&w);
+
+ return w._tooclose;
+}
+
+
+
+void BoundaryLayerColumns::filterPoints(GEntity *ge, double factor)
+#if defined(HAVE_RTREE)
+{
+ // return;
+ // compute the element sizes
+ std::map<MVertex*,double> sizes;
+ if (ge->dim() == 2){
+ backgroundMesh::set((GFace *)ge);
+ std::list<GEdge*> edges = ge->edges();
+ std::list<GEdge*>::iterator it = edges.begin();
+ for ( ; it != edges.end() ; ++it){
+ GEdge *ged = *it;
+ for (unsigned int i=0;i<ged->lines.size();i++){
+ MLine *e = ged->lines[i];
+ MVertex *v0 = e->getVertex(0);
+ MVertex *v1 = e->getVertex(1);
+ double d = v0->distance(v1);
+ std::map<MVertex*,double>::iterator it0 = sizes.find(v0);
+ if (it0 == sizes.end()) sizes[v0] = d;
+ else it0->second = std::max(d, it0->second);
+ std::map<MVertex*,double>::iterator it1 = sizes.find(v1);
+ if (it1 == sizes.end()) sizes[v1] = d;
+ else it1->second = std::max(d, it1->second);
+ }
+ }
+ }
+ else {
+ Msg::Fatal("code ce truc JF !");
+ }
+
+ // a RTREE data structure that enables to verify if
+ // points are too close
+ RTree<blPoint_rtree*,double,3,double> rtree;
+ // stores the info "where the new vertex comes form"
+ std::map <MVertex*, MVertex*> v2v;
+
+ // compute maximum column size
+ // initialize the RTREE with points on the boundary
+ unsigned int MAXCOLSIZE = 0;
+ BoundaryLayerColumns::iter it = _data.begin();
+ for ( ; it != _data.end() ; ++it) {
+ BoundaryLayerData & d = it->second;
+ MAXCOLSIZE = MAXCOLSIZE > d._column.size() ? MAXCOLSIZE : d._column.size();
+ MVertex * v = it->first;
+ double largeMeshSize = factor*sizes[v];
+ blPoint_rtree *p = new blPoint_rtree(v,largeMeshSize);
+ double _min[3],_max[3];
+ p->minmax (_min,_max);
+ rtree.Insert(_min,_max,p);
+ v2v[v] = v;
+ for (unsigned int k = 0 ; k < d._column.size() ; k++)
+ v2v[d._column[k]] = v;
+ }
+
+ // go layer by layer
+ for (unsigned int LAYER = 0 ; LAYER < MAXCOLSIZE ; LAYER++){
+ // store accepted points that will be inserted in the rtree
+ // afterwards
+ std::set<MVertex*> accepted;
+ it = _data.begin();
+ for ( ; it != _data.end() ; ++it) {
+ MVertex * v = it->first;
+ double largeMeshSize = sizes[v];
+ BoundaryLayerData & d = it->second;
+ // take the point if the number of layers is
+ // large enough
+ if (d._column.size() > LAYER){
+ // check if the vertex in the column at position LAYER
+ // isn't too close to another vertex
+ MVertex *toCheck = d._column[LAYER];
+ if (LAYER){
+ double DD = toCheck->distance ( d._column[LAYER-1] );
+ // do not allow to have elements that are stretched the
+ // other way around !
+ if (DD > largeMeshSize) largeMeshSize *= 100;
+ largeMeshSize = std::max (largeMeshSize, DD);
+ }
+ largeMeshSize *= factor;
+ bool exclude = inExclusionZone_filter (toCheck,v2v,rtree);
+ if (!exclude){
+ v2v [toCheck] = v;
+ blPoint_rtree *p = new blPoint_rtree(toCheck,largeMeshSize);
+ double _min[3],_max[3];
+ p->minmax (_min,_max);
+ rtree.Insert(_min,_max,p);
+ }
+ else {
+ std::vector<MVertex*> newColumn;
+ for (unsigned int k = 0 ; k < LAYER ; k++) newColumn.push_back(d._column[k]);
+ for (unsigned int k = LAYER ; k < d._column.size() ; k++) delete d._column[k];
+ d._column = newColumn;
+ }
+ }
+ }
+ }
+#else
+ Msg::Warning ("Boundary Layer Points cannot be filtered without compiling gmsh with the rtree library");
+#endif
+}
#endif
diff --git a/Geo/boundaryLayersData.h b/Geo/boundaryLayersData.h
index 261ac7b058..ee02ed8516 100644
--- a/Geo/boundaryLayersData.h
+++ b/Geo/boundaryLayersData.h
@@ -205,7 +205,7 @@ public:
static BoundaryLayerData error;
return error;
}
- void filterPoints();
+ void filterPoints(GEntity *ge, double factor);
};
BoundaryLayerField* getBLField(GModel *gm);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 2fad6b8759..91ee53ec0b 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -452,7 +452,6 @@ static bool algoDelaunay2D(GFace *gf)
if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_BAMG ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL ||
- gf->getMeshingAlgo() == ALGO_2D_RUPPERT ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
gf->getMeshingAlgo() == ALGO_2D_BAMG)
@@ -1297,8 +1296,6 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf);
- else if(gf->getMeshingAlgo() == ALGO_2D_RUPPERT)
- gmshRuppert(gf,.5);
else {
bowyerWatson(gf,15000);
meshGFaceBamg(gf);
@@ -2031,8 +2028,6 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf,1000000000, &equivalence, ¶metricCoordinates);
- else if(gf->getMeshingAlgo() == ALGO_2D_RUPPERT)
- gmshRuppert(gf,.4,10000, &equivalence, ¶metricCoordinates);
else
meshGFaceBamg(gf);
if (!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
@@ -2114,7 +2109,6 @@ void meshGFace::operator() (GFace *gf, bool print)
case ALGO_2D_DELAUNAY : algo = "Delaunay"; break;
case ALGO_2D_MESHADAPT_OLD : algo = "MeshAdapt (old)"; break;
case ALGO_2D_BAMG : algo = "Bamg"; break;
- case ALGO_2D_RUPPERT : algo = "Ruppert"; break;
case ALGO_2D_PACK_PRLGRMS : algo = "Square Packing"; break;
case ALGO_2D_AUTO :
algo = (gf->geomType() == GEntity::Plane) ? "Delaunay" : "MeshAdapt";
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 19b59aa852..6103fe15bf 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -42,6 +42,30 @@ static bool isBoundary(MTri3 *t, double limit_, int &active)
return false;
}
*/
+
+struct equivalentTriangle {
+ MTri3 *_t;
+ MVertex *_v[3];
+ equivalentTriangle (MTri3 *t, std::map<MVertex* , MVertex*>* equivalence)
+ :_t(t) {
+ for (int i=0;i<3;i++){
+ MVertex *v = t->tri()->getVertex(i);
+ std::map<MVertex* , MVertex*>::iterator it = equivalence->find(v);
+ if (it == equivalence->end())_v[i] = v;
+ else _v[i] = it->second;
+ }
+ std::sort (_v,_v+3);
+ }
+ bool operator < (const equivalentTriangle &other) const{
+ for (int i=0;i<3;i++){
+ if (other._v[i] > _v[i])return true;
+ if (other._v[i] < _v[i])return false;
+ }
+ return false;
+ }
+};
+
+
template <class ITERATOR>
void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData & data, bool param=true)
{
@@ -227,6 +251,7 @@ void buildMetric(GFace *gf, double *uv, double *metric)
metric[0] = dot(der.first(), der.first());
metric[1] = dot(der.second(), der.first());
metric[2] = dot(der.second(), der.second());
+
}
// m 3x3
@@ -460,6 +485,102 @@ void connectTriangles(std::set<MTri3*, compareTri3Ptr> &l)
connectTris(l.begin(), l.end());
}
+int inCircumCircleTangentPlane(MTriangle *t,
+ SPoint3 &p, SVector3 &t1, SVector3 &t2)
+{
+ double p1[2],p2[2],p3[2], pp[2] = {0,0};
+ MVertex *v1 = t->getVertex(0);
+ MVertex *v2 = t->getVertex(1);
+ MVertex *v3 = t->getVertex(2);
+ SVector3 pmx1 (v1->x()-p.x(),v1->y()-p.y(),v1->z()-p.z());
+ p1[0] = dot(pmx1,t1);p1[1] = dot(pmx1,t2);
+ SVector3 pmx2 (v2->x()-p.x(),v2->y()-p.y(),v2->z()-p.z());
+ p2[0] = dot(pmx2,t1);p2[1] = dot(pmx2,t2);
+ SVector3 pmx3 (v3->x()-p.x(),v3->y()-p.y(),v3->z()-p.z());
+ p3[0] = dot(pmx3,t1);p3[1] = dot(pmx3,t2);
+ 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,
+ MTri3 *t,
+ SPoint3 &p, SVector3 &t1, SVector3 &t2)
+{
+ t->setDeleted(true);
+ 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 = inCircumCircleTangentPlane(neigh->tri(), p, t1, t2);
+ if (circ)
+ recurFindCavityTangentPlane(shell, cavity, neigh, p, t1, t2);
+ else
+ shell.push_back(edgeXface(t, i));
+ }
+ }
+}
+
+///*********************
+/// T E S T : project the cavity on the tangent plane in case of bad mapping
+
+static void computeTangentPlane (GFace *gf, double center[2], SPoint3 &p, SVector3 &t1, SVector3 &t2)
+{
+ GPoint gp = gf->point(SPoint2(center[0],center[1]));
+ p = SPoint3(gp.x(),gp.y(),gp.z());
+ Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(center[0], center[1]));
+ SVector3 n = crossprod(der.first(),der.second());
+ t1 = der.first();
+ t2 = crossprod(n,der.first());
+ t1.normalize();
+ t2.normalize();
+ // printf("%g %g %g :: %g %g %g\n",t1.x(),t1.y(),t1.z(),t2.x(),t2.y(),t2.z());
+}
+
+bool findCavityTangentPlane(GFace *gf, double *center,
+ std::list<edgeXface> &shell,
+ std::list<MTri3*> &cavity,
+ MTri3 *t)
+{
+ return false;
+ SPoint3 p; SVector3 t1,t2;
+ computeTangentPlane (gf,center, p,t1,t2);
+ SVector3 N = crossprod(t2,t1);
+ N.normalize();
+ const double d = -(N.x()*p.x()+N.y()*p.y()+N.z()*p.z());
+ recurFindCavityTangentPlane(shell, cavity, t, p, t1, t2);
+ // double AMIN = 1;
+ double DMAX = 0.0;
+ for (std::list<MTri3*>::iterator i=cavity.begin();i!=cavity.end();i++){
+ t = *i;
+ SPoint3 b = t->tri()->getFace(0).barycenter();
+ // SVector3 n = t->tri()->getFace(0).normal();
+ // double a = fabs(dot(N,n));
+ double dist = fabs(N.x()*b.x()+N.y()*b.y()+N.z()*b.z() + d);
+ DMAX = std::max(DMAX,dist);
+ // AMIN = std::min(a,AMIN);
+ // printf("%g %g %g -- %g %g %g\n",N.x(),N.y(),N.z(),n.x(),n.y(),n.z());
+ }
+ if (DMAX > 30 || cavity.size() < 2){
+ printf("%d elements in the cavity DMAX %g\n",cavity.size(),DMAX);
+ for (std::list<MTri3*>::iterator i=cavity.begin();i!=cavity.end();i++){
+ t = *i;
+ t->setDeleted(false);
+ }
+ cavity.clear();
+ shell.clear();
+ return false;
+ }
+ return true;
+}
+
+
+
void recurFindCavity(std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
double *v, double *param, MTri3 *t, bidimMeshData & data)
{
@@ -585,7 +706,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
double *metric,
MTri3 **oneNewTriangle)
{
- if (shell.size() <= 3 || shell.size() != cavity.size() + 2) return false;
+ if (shell.size() != cavity.size() + 2) return false;
std::list<MTri3*> new_cavity;
@@ -773,15 +894,20 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2], bidimMeshData & data,
return 0;
}
-//double __DT1;
-static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it,
- double center[2], double metric[3], bidimMeshData & data,
+
+///*********************
+
+static bool insertAPoint(GFace *gf,
+ std::set<MTri3*,compareTri3Ptr>::iterator it,
+ double center[2],
+ double metric[3],
+ bidimMeshData & data,
std::set<MTri3*,compareTri3Ptr> &AllTris,
std::set<MTri3*,compareTri3Ptr> *ActiveTris = 0,
- MTri3 *worst = 0, MTri3 **oneNewTriangle = 0)
+ MTri3 *worst = 0,
+ MTri3 **oneNewTriangle = 0)
{
-
if (worst){
it = AllTris.find(worst);
if (worst != *it){
@@ -796,11 +922,11 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
std::list<MTri3*> cavity;
double uv[2];
+ // TEST
// if the point is able to break the bad triangle "worst"
if (inCircumCircleAniso(gf, worst->tri(), center, metric, data)){
- // double t1 = Cpu();
- recurFindCavityAniso(gf, shell, cavity, metric, center, worst, data);
- // __DT1 += Cpu() - t1 ;
+ if (!findCavityTangentPlane(gf,center,shell, cavity, worst))
+ recurFindCavityAniso(gf, shell, cavity, metric, center, worst, data);
for (std::list<MTri3*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
if (invMapUV((*itc)->tri(), center, data, uv, 1.e-8)) {
ptin = *itc;
@@ -808,13 +934,14 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
}
- // else look for it
- else {
- // printf("cocuou\n");
+ else {
ptin = search4Triangle (worst, center, data, AllTris,uv, oneNewTriangle ? true : false);
- if (ptin) {
+ // printf("what's this %g %g\n",center[0],center[1]);
+
+ if (ptin) {
+ if (!findCavityTangentPlane(gf,center,shell, cavity, ptin))
recurFindCavityAniso(gf, shell, cavity, metric, center, ptin, data);
- }
+ }
}
if (ptin) {
@@ -870,55 +997,38 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
-void gmshRuppert(GFace *gf, double minqual, int MAXPNT,
- std::map<MVertex* , MVertex*>* equivalence,
- std::map<MVertex*, SPoint2> * parametricCoordinates){
- MTri3::radiusNorm =3;
-
- std::set<MTri3*,compareTri3Ptr> AllTris;
- bidimMeshData DATA(equivalence,parametricCoordinates);
-
- buildMeshGenerationDataStructures(gf, AllTris, DATA);
-
- int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, DATA);
-
- Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
-
- if(AllTris.empty()){
- Msg::Error("No triangles in initial mesh");
- return;
- }
-
- int NBDELETED = 0;
- while (1){
- MTri3 *worst = *AllTris.begin();
- if (worst->isDeleted()){
- delete worst->tri();
- delete worst;
- AllTris.erase(AllTris.begin());
- NBDELETED ++;
+bool computeEquivalentTriangles (std::set<MTri3*,compareTri3Ptr> &AllTris,
+ std::map<MVertex* , MVertex*>* equivalence) {
+ std::vector<MTri3*> WTF;
+ if (!equivalence)return false;
+ std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
+ std::set<equivalentTriangle> eqTs;
+ int COUNT = 0;
+ for ( ; it!=AllTris.end();++it){
+ if (!(*it)->isDeleted()){
+ COUNT++;
+ equivalentTriangle et ((*it),equivalence);
+ std::set<equivalentTriangle>::iterator iteq = eqTs.find(et);
+ if (iteq == eqTs.end())eqTs.insert(et);
+ else {
+ WTF.push_back(iteq->_t);
+ WTF.push_back(*it);
+ }
}
- else{
- double center[2],metric[3],r2;
- // printf("%12.5E\n",worst->getRadius() );
- if (1./worst->getRadius() > minqual || (int) DATA.vSizes.size() > MAXPNT) break;
- circUV(worst->tri(), DATA, center, gf);
- MTriangle *base = worst->tri();
- int index0 = DATA.getIndex( base->getVertex(0) );
- int index1 = DATA.getIndex( base->getVertex(1) );
- int index2 = DATA.getIndex( base->getVertex(2) );
- double pa[2] = {(DATA.Us[index0] + DATA.Us[index1] + DATA.Us[index2])/ 3.,
- (DATA.Vs[index0] + DATA.Vs[index1] + DATA.Vs[index2])/ 3.};
- buildMetric(gf, pa, metric);
- circumCenterMetric(worst->tri(), metric, DATA, center, r2);
- insertAPoint(gf, AllTris.begin(), center, metric, DATA, AllTris);
+ }
+
+ if (WTF.size()){
+ for (int i=0;i<WTF.size();i++){
+ std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.find(WTF[i]);
+ AllTris.erase(it);
+ WTF[i]->forceRadius(100);
+ AllTris.insert(WTF[i]);
}
+ return true;
}
- MTri3::radiusNorm =2;
- transferDataStructure(gf, AllTris, DATA);
+ return false;
}
-
void bowyerWatson(GFace *gf, int MAXPNT,
std::map<MVertex* , MVertex*>* equivalence,
std::map<MVertex*, SPoint2> * parametricCoordinates)
@@ -937,7 +1047,7 @@ void bowyerWatson(GFace *gf, int MAXPNT,
Msg::Error("No triangles in initial mesh");
return;
}
-
+
int ITER = 0;
int NBDELETED = 0;
// double DT1 = 0 , DT2=0, DT3=0;
@@ -963,8 +1073,12 @@ void bowyerWatson(GFace *gf, int MAXPNT,
DATA.vSizes.size(), worst->getRadius());
// printf("%d %d %d\n",vSizes.size(), AllTris.size(),NBDELETED);
}
+ // VERIFY STOP !!!
+ if (worst->getRadius() < 0.5 * sqrt(2.0) || (int) DATA.vSizes.size() > MAXPNT) {
+ break;
+ }
+
double center[2],metric[3],r2;
- if (worst->getRadius() < 0.5 * sqrt(2.0) || (int) DATA.vSizes.size() > MAXPNT) break;
circUV(worst->tri(), DATA, center, gf);
MTriangle *base = worst->tri();
int index0 = DATA.getIndex( base->getVertex(0) );
@@ -972,14 +1086,15 @@ void bowyerWatson(GFace *gf, int MAXPNT,
int index2 = DATA.getIndex( base->getVertex(2) );
double pa[2] = {(DATA.Us[index0] + DATA.Us[index1] + DATA.Us[index2])/ 3.,
(DATA.Vs[index0] + DATA.Vs[index1] + DATA.Vs[index2])/ 3.};
+
buildMetric(gf, pa, metric);
circumCenterMetric(worst->tri(), metric, DATA, center, r2);
// DT2 += (Cpu() - t2) ;
// double t3 = Cpu() ;
insertAPoint(gf, AllTris.begin(), center, metric, DATA, AllTris);
- // DT3 += (Cpu() - t3) ;
}
}
+ nbSwaps = edgeSwapPass(gf, AllTris, SWCR_QUAL, DATA);
// printf("%12.5E %12.5E %12.5E %12.5E %12.5E\n",DT1,DT2,DT3,__DT1,__DT2);
// printf("%12.5E \n",__DT2);
#if defined(HAVE_ANN)
@@ -994,6 +1109,7 @@ void bowyerWatson(GFace *gf, int MAXPNT,
}
#endif
transferDataStructure(gf, AllTris, DATA);
+ removeThreeTrianglesNodes(gf);
}
/*
@@ -1206,8 +1322,8 @@ void optimalPointFrontalB (GFace *gf,
}
void bowyerWatsonFrontal(GFace *gf,
- std::map<MVertex* , MVertex*>* equivalence,
- std::map<MVertex*, SPoint2> * parametricCoordinates)
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
std::set<MTri3*,compareTri3Ptr> ActiveTris;
@@ -1232,13 +1348,13 @@ void bowyerWatsonFrontal(GFace *gf,
int ITERATION = 0;
while (1){
++ITERATION;
- // 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);
- // sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
- // _printTris (name, ActiveTris.begin(), ActiveTris.end(), DATA,true);
- // }
+ 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);
+ sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ _printTris (name, ActiveTris.begin(), ActiveTris.end(), DATA,true);
+ }
/* if(ITER % 100== 0){
char name[245];
sprintf(name,"delfr2d%d-ITER%d.pos",gf->tag(),ITER);
@@ -1276,7 +1392,9 @@ void bowyerWatsonFrontal(GFace *gf,
// _printTris (name, AllTris, Us, Vs,false);
// sprintf(name,"frontal%d-param.pos", gf->tag());
// _printTris (name, AllTris, Us, Vs,true);
+ nbSwaps = edgeSwapPass(gf, AllTris, SWCR_QUAL, DATA);
transferDataStructure(gf, AllTris, DATA);
+ removeThreeTrianglesNodes(gf);
// in case of boundary layer meshing
#if defined(HAVE_ANN)
{
@@ -1411,7 +1529,7 @@ void buildBackGroundMesh (GFace *gf,
int CurvControl = CTX::instance()->mesh.lcFromCurvature;
CTX::instance()->mesh.lcFromCurvature = 0;
// Do a background mesh
- gmshRuppert(gf,0.3,4000, equivalence,parametricCoordinates);
+ bowyerWatson(gf,4000, equivalence,parametricCoordinates);
// Re-enable curv control if asked
CTX::instance()->mesh.lcFromCurvature = CurvControl;
// apply this to the BGM
@@ -1612,7 +1730,7 @@ void bowyerWatsonParallelograms(GFace *gf,
// buildMetric(gf, newPoint, metrics[i], metric);
buildMetric(gf, newPoint, metric);
- bool success = insertAPoint( gf, AllTris.begin(), newPoint, metric, DATA , AllTris, 0, oneNewTriangle, &oneNewTriangle);
+ bool success = insertAPoint(gf, AllTris.begin(), newPoint, metric, DATA , AllTris, 0, oneNewTriangle, &oneNewTriangle);
if (!success) oneNewTriangle = 0;
// if (!success)printf("success %d %d\n",success,AllTris.size());
i++;
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index ee07a19ad3..96260fc3a5 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -134,9 +134,6 @@ class compareTri3Ptr
void connectTriangles(std::list<MTri3*> &);
void connectTriangles(std::vector<MTri3*> &);
void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
-void gmshRuppert(GFace *gf, double minqual = 0.2, int MAXPNT= 1000000000,
- std::map<MVertex* , MVertex*>* equivalence = 0,
- std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
void bowyerWatson(GFace *gf, int MAXPNT= 1000000000,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index a26fdbbb19..71cfe76607 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -359,6 +359,63 @@ double surfaceTriangleUV(MVertex *v1, MVertex *v2, MVertex *v3,bidimMeshData & d
return 0.5 * fabs (v12[0] * v13[1] - v12[1] * v13[0]);
}
+int _removeThreeTrianglesNodes(GFace *gf)
+{
+ v2t_cont adj;
+ buildVertexToElement(gf->triangles,adj);
+ v2t_cont :: iterator it = adj.begin();
+ int n=0;
+ std::set<MElement*> touched;
+ while (it != adj.end()) {
+ bool skip = false;
+ if(it->second.size() == 3 && it->first->onWhat()->dim() == 2) {
+ const std::vector<MElement*> < = it->second;
+ std::set<MVertex*> vs;
+ for(int i = 0; i < 3; i++) {
+ if(touched.find(lt[i])!=touched.end() || lt[i]->getNumVertices()!=3){
+ skip=true;
+ break;
+ }
+ for(int j = 0; j < 3; j++) {
+ if(lt[i]->getVertex(j) == it->first) {
+ vs.insert(lt[i]->getVertex((j+1)%3));
+ vs.insert(lt[i]->getVertex((j+2)%3));
+ break;
+ }
+ }
+ }
+ if(skip){
+ it++;
+ continue;
+ }
+ std::set<MVertex*>::iterator itt = vs.begin();
+ MVertex *v1 = *itt; ++itt;
+ MVertex *v2 = *itt; ++itt;
+ MVertex *v3 = *itt;
+ MTriangle *newt = new MTriangle(v1,v2,v3);
+ n++;
+ gf->triangles.push_back(newt);
+ for(int i=0;i<3;i++) {
+ touched.insert(lt[i]);
+ }
+ }
+ it++;
+ }
+ std::vector<MTriangle*> triangles2;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ if(touched.find(gf->triangles[i]) == touched.end()){
+ triangles2.push_back(gf->triangles[i]);
+ }
+ else {
+ delete gf->triangles[i];
+ }
+ }
+ gf->triangles = triangles2;
+ Msg::Debug("%i three-triangles vertices removed",n);
+ return n;
+}
+
+
int _removeFourTrianglesNodes(GFace *gf,bool replace_by_quads)
{
v2t_cont adj;
@@ -369,7 +426,7 @@ int _removeFourTrianglesNodes(GFace *gf,bool replace_by_quads)
while (it != adj.end()) {
bool skip = false;
double surfaceRef = 0;
- if(it->second.size() == 4) {
+ if(it->second.size() == 4 && it->first->onWhat()->dim() == 2) {
const std::vector<MElement*> < = it->second;
MVertex* edges[4][2];
for(int i = 0; i < 4; i++) {
@@ -463,11 +520,19 @@ int _removeFourTrianglesNodes(GFace *gf,bool replace_by_quads)
return n;
}
+
+
void removeFourTrianglesNodes(GFace *gf,bool replace_by_quads)
{
while(_removeFourTrianglesNodes(gf,replace_by_quads));
}
+void removeThreeTrianglesNodes(GFace *gf)
+{
+ while(_removeThreeTrianglesNodes(gf));
+}
+
+
/*
+--------+
@@ -2131,10 +2196,6 @@ static int _relocateVertexOpti(GFace *gf, MVertex *ver,
void _relocateVertex(GFace *gf, MVertex *ver,
const std::vector<MElement*> <)
{
- double R;
- SPoint3 c;
- bool isSphere = gf->isSphere(R, c);
-
if(ver->onWhat()->dim() != 2) return;
MFaceVertex *fv = dynamic_cast<MFaceVertex*>(ver);
if(fv && fv->bl_data) return;
@@ -2142,51 +2203,60 @@ void _relocateVertex(GFace *gf, MVertex *ver,
double initu,initv;
ver->getParameter(0, initu);
ver->getParameter(1, initv);
- double cu = 0, cv = 0;
- double XX=0,YY=0,ZZ=0;
- double pu[4], pv[4];
- double fact = 0.0;
+
+ // compute the vertices connected to that one
+ std::map<MVertex*,SPoint2> pts;
for(unsigned int i = 0; i < lt.size(); i++){
- parametricCoordinates(lt[i], gf, pu, pv, ver);
- double XCG=0,YCG=0,ZCG=0;
- for (int j=0;j<lt[i]->getNumVertices();j++){
- XCG += lt[i]->getVertex(j)->x();
- YCG += lt[i]->getVertex(j)->y();
- ZCG += lt[i]->getVertex(j)->z();
- }
- XX += XCG;
- YY += YCG;
- ZZ += ZCG;
-
- // double D = 1./sqrt((XCG-ver->x())*(XCG-ver->x()) +
- // (YCG-ver->y())*(YCG-ver->y()) +
- // (ZCG-ver->z())*(ZCG-ver->z()) );
- double D = 1.0;
- for (int j=0;j<lt[i]->getNumVertices();j++){
- cu += pu[j]*D;
- cv += pv[j]*D;
- }
- fact += lt[i]->getNumVertices() * D;
- }
- SPoint2 newp ;
- if(fact != 0.0){
- SPoint2 before(initu,initv);
- SPoint2 after(cu / fact,cv / fact);
- if (isSphere){
- GPoint gp = gf->closestPoint(SPoint3(XX/fact, YY/fact, ZZ/fact), after);
- after = SPoint2(gp.u(),gp.v());
- }
- bool success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,before,after);
+ for (int j=0;j<lt[i]->getNumEdges();j++){
+ MEdge e = lt[i]->getEdge(j);
+ SPoint2 param0, param1;
+ if (e.getVertex(0) == ver){
+ reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0, param1);
+ pts[e.getVertex(1)] = param1;
+ }
+ else if (e.getVertex(1) == ver){
+ reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0, param1);
+ pts[e.getVertex(0)] = param0;
+ }
+ }
+ }
+
+ SPoint2 before(initu,initv);
+ double metric[3];
+ SPoint2 after(0,0);
+ double COUNT = 0.0;
+ // printf("weights :");
+ for(std::map<MVertex*,SPoint2>::iterator it = pts.begin(); it != pts.end() ; ++it) {
+ SPoint2 adj = it->second;
+ SVector3 d (adj.x()-before.x(),adj.y()-before.y(),0.0);
+ d.normalize();
+ buildMetric (gf,adj,metric);
+ const double F = sqrt(metric[0]*d.x()*d.x() +
+ 2*metric[1]*d.x()*d.y() +
+ metric[2]*d.y()*d.y());
+ // printf("%g ",F);
+ after += adj*F;
+ COUNT += F;
+ // // double RATIO = lt[i]->getVolume()/pow(metric[0]*metric[2]-metric[1]*metric[1],0.5);
+ }
+ // printf("\n");
+ after *= (1.0/COUNT);
+ double FACTOR = 1.0;
+ const int MAXITER = 5;
+ for (int ITER = 0;ITER < MAXITER; ITER ++){
+ SPoint2 trial = after * FACTOR + before * (1.-FACTOR);
+ bool success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,before,trial);
if (success){
- ver->setParameter(0, after.x());
- ver->setParameter(1, after.y());
- GPoint pt = gf->point(after);
+ ver->setParameter(0, trial.x());
+ ver->setParameter(1, trial.y());
+ GPoint pt = gf->point(trial);
if(pt.succeeded()){
- ver->x() = pt.x();
- ver->y() = pt.y();
- ver->z() = pt.z();
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
}
}
+ FACTOR /= 1.4;
}
}
@@ -2651,6 +2721,20 @@ int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*>
return 0;
}
+bool edgeSwapDelProj (MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4){
+ MTriangle t1(v1,v2,v3);
+ MTriangle t2(v2,v1,v4);
+
+ SVector3 n1 = t1.getFace(0).normal();
+ SVector3 n2 = t2.getFace(0).normal();
+ if (dot(n1,n2) <= 0) {
+ // printf("OLA !!!\n");
+ return true;
+ }
+ return false;
+}
+
+
bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
std::vector<MTri3*> &newTris, const swapCriterion &cr, bidimMeshData & data)
{
@@ -2668,22 +2752,26 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
swapquad sq (v1, v2, v3, v4);
if(configs.find(sq) != configs.end()) return false;
configs.insert(sq);
+
+ if (edgeSwapDelProj(v3,v4,v2,v1))return false;
- const double volumeRef = surfaceTriangleUV(v1, v2, v3, data) +
- surfaceTriangleUV(v1, v2, v4, data);
+ // const double volumeRef = surfaceTriangleUV(v1, v2, v3, data) +
+ // surfaceTriangleUV(v1, v2, v4, data);
MTriangle *t1b = new MTriangle(v2, v3, v4);
MTriangle *t2b = new MTriangle(v4, v3, v1);
- const double v1b = surfaceTriangleUV(v2, v3, v4, data);
- const double v2b = surfaceTriangleUV(v4, v3, v1, data);
- const double volume = v1b + v2b;
- if(fabs(volume - volumeRef) > 1.e-10 * (volume + volumeRef) ||
- v1b < 1.e-8 * (volume + volumeRef) ||
- v2b < 1.e-8 * (volume + volumeRef)){
- delete t1b;
- delete t2b;
- return false;
- }
+
+
+ // const double v1b = surfaceTriangleUV(v2, v3, v4, data);
+ // const double v2b = surfaceTriangleUV(v4, v3, v1, data);
+ // const double volume = v1b + v2b;
+ // if(fabs(volume - volumeRef) > 1.e-10 * (volume + volumeRef) ||
+ // v1b < 1.e-8 * (volume + volumeRef) ||
+ // v2b < 1.e-8 * (volume + volumeRef)){
+ // delete t1b;
+ // delete t2b;
+ // return false;
+ // }
switch(cr){
case SWCR_QUAL:
@@ -2692,11 +2780,14 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
qmTriangle(t2->tri(), QMTRI_RHO));
const double triQuality = std::min(qmTriangle(t1b, QMTRI_RHO),
qmTriangle(t2b, QMTRI_RHO));
- if(triQuality < triQualityRef){
- delete t1b;
- delete t2b;
- return false;
+ if (!edgeSwapDelProj(v1,v2,v3,v4)){
+ if(triQuality < triQualityRef){
+ delete t1b;
+ delete t2b;
+ return false;
+ }
}
+ // printf("coucou\n");
break;
}
case SWCR_DEL:
@@ -2884,9 +2975,8 @@ bool buildVertexCavity(MTri3 *t, int iLocalVertex, MVertex **v1,
}
}
-// split one triangle into 3 triangles then apply edge swop (or not)
-// will do better (faster) soon, just to test
-void _triangleSplit (GFace *gf, MElement *t, bool swop = false)
+// split one triangle into 3 triangles
+void _triangleSplit (GFace *gf, MElement *t)
{
MVertex *v1 = t->getVertex(0);
MVertex *v2 = t->getVertex(1);
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 3c6849f570..5d363c1917 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -75,6 +75,7 @@ int edgeSwapPass(GFace *gf,
std::set<MTri3*, compareTri3Ptr> &allTris,
const swapCriterion &cr, bidimMeshData &DATA);
void removeFourTrianglesNodes(GFace *gf, bool replace_by_quads);
+void removeThreeTrianglesNodes(GFace *gf);
void buildMeshGenerationDataStructures(GFace *gf,
std::set<MTri3*, compareTri3Ptr> &AllTris,
bidimMeshData & data);
diff --git a/benchmarks/2d/Square-01.geo b/benchmarks/2d/Square-01.geo
index 5078de8447..42006dad03 100644
--- a/benchmarks/2d/Square-01.geo
+++ b/benchmarks/2d/Square-01.geo
@@ -1,13 +1,13 @@
fact = 100;
lc = .1 * fact;
-Point(1) = {0.0,0.0,0,lc*.3};
+Point(1) = {0.0,0.0,0,lc*1};
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};
+Point(4) = {0,1* fact,0,lc*1};
Line(1) = {3,2};
Line(2) = {2,1};
Line(3) = {1,4};
Line(4) = {4,3};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
-Recombine Surface {6};
+//Recombine Surface {6};
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index bf9543f499..ce88526a7d 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -78,5 +78,5 @@ Physical Line(25) = {12,13,14,15,16,17,18,19,20,10};
Physical Surface(26) = {24,22};
Physical Line(27) = {10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 14, 13, 12, 11};
Physical Line(28) = {17, 16, 20, 19, 18, 15};
-Recombine Surface {24, 22};
-Mesh.Algorithm=8;
+//Recombine Surface {24, 22};
+//Mesh.Algorithm=8;
--
GitLab