From 54befb5842d893c1f7dd225e5eaa0bea9da09524 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Thu, 17 Jan 2008 17:48:39 +0000
Subject: [PATCH] Aniso delaunay surface mesher is now fast : the initial mesh
is now "delaunized" using swaps using a "fast" algorithm
---
Geo/OCCFace.cpp | 3 +-
Mesh/BDS.cpp | 6 +--
Mesh/meshGFace.cpp | 6 +--
Mesh/meshGFaceDelaunayInsertion.cpp | 10 ++--
Mesh/meshGFaceOptimize.cpp | 80 +++++++++++++++++++++--------
5 files changed, 72 insertions(+), 33 deletions(-)
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 7b05b04e28..2e3182fc2e 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.24 2007-11-11 19:53:57 remacle Exp $
+// $Id: OCCFace.cpp,v 1.25 2008-01-17 17:48:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -52,6 +52,7 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
GEdge *e = m->edgeByTag(index);
if(!e) throw;
l_wire.push_back(e);
+ Msg(DEBUG2,"Edge %d",e->tag());
e->addFace(this);
if(!e->is3D()){
OCCEdge *occe = (OCCEdge*)e;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 1d4fc36ac9..00131e7940 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.88 2008-01-15 19:50:58 remacle Exp $
+// $Id: BDS.cpp,v 1.89 2008-01-17 17:48:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -44,7 +44,7 @@ void outputScalarField(std::list < BDS_Face * >t, const char *iii, int param)
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
pts[0]->u, pts[0]->v, 0.0,
pts[1]->u, pts[1]->v, 0.0,
- pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->lc(),(double)pts[1]->lc(),(double)pts[2]->lc());
+ pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
else
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
@@ -1326,7 +1326,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf)
double sold = fabs(surface_triangle_param(n[0],n[1],n[2]));
s2 += sold;
// printf("%22.15E %22.15E\n",snew,sold);
- // if ( snew < .1 * sold) return false;
+ if ( snew < .1 * sold) return false;
// if (!test_move_point_parametric_triangle ( p, U, V, t)){
// return false;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 74e28e80d9..9a4670f893 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.105 2008-01-15 19:50:58 remacle Exp $
+// $Id: meshGFace.cpp,v 1.106 2008-01-17 17:48:38 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -2253,8 +2253,8 @@ void meshGFace::operator() (GFace *gf)
{
Msg(DEBUG1, "Generating the mesh");
if(noseam (gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
- //gmsh2DMeshGenerator(gf,0, true);
- gmsh2DMeshGenerator(gf,0, false);
+ gmsh2DMeshGenerator(gf,0, true);
+ //gmsh2DMeshGenerator(gf,0, false);
}
else{
if(!gmsh2DMeshGeneratorPeriodic(gf,false))
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 42cd02777d..c8624ed87b 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.7 2008-01-14 21:29:14 remacle Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.8 2008-01-17 17:48:39 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -584,15 +584,13 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
gf->triangles.clear();
connectTris ( AllTris.begin(), AllTris.end() );
+ Msg(DEBUG,"All %d tris were connected",AllTris.size());
// _printTris ("before.pos", AllTris, Us,Vs);
- // this should be MUCH faster !
- for (int i=0;i<1200;i++)
- if(!edgeSwapPass(gf,AllTris,SWCR_DEL,Us,Vs,vSizes,vSizesBGM))break;
-
+ int nbSwaps = edgeSwapPass(gf,AllTris,SWCR_DEL,Us,Vs,vSizes,vSizesBGM);
// _printTris ("after2.pos", AllTris, Us,Vs);
+ Msg(DEBUG,"Delaunization of the initial mesh done (%d swaps)",nbSwaps);
- Msg(DEBUG,"All %d tris were connected",AllTris.size());
// here the classification should be done
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index f460af3b7c..d20d091e6a 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -127,9 +127,29 @@ double surfaceTriangleUV (MVertex *v1, MVertex *v2, MVertex *v3,
return 0.5*fabs (v12[0]*v13[1]-v12[1]*v13[0]);
}
+struct swapquad{
+ int v[4];
+ bool operator < (const swapquad &o) const{
+ if (v[0] < o.v[0])return true;
+ if (v[0] > o.v[0])return false;
+ if (v[1] < o.v[1])return true;
+ if (v[1] > o.v[1])return false;
+ if (v[2] < o.v[2])return true;
+ if (v[2] > o.v[2])return false;
+ if (v[3] < o.v[3])return true;
+ return false;
+ }
+ swapquad(MVertex *v1,MVertex *v2,MVertex *v3,MVertex *v4){
+ v[0] = v1->getNum();
+ v[1] = v2->getNum();
+ v[2] = v3->getNum();
+ v[3] = v4->getNum();
+ std::sort(v,v+4);
+ }
+};
-
-bool gmshEdgeSwap(MTri3 *t1,
+bool gmshEdgeSwap(std::set<swapquad> & configs,
+ MTri3 *t1,
GFace *gf,
int iLocalEdge,
std::vector<MTri3*> &newTris,
@@ -137,7 +157,7 @@ bool gmshEdgeSwap(MTri3 *t1,
const std::vector<double> & Us,
const std::vector<double> & Vs,
const std::vector<double> & vSizes ,
- const std::vector<double> & vSizesBGM ){
+ const std::vector<double> & vSizesBGM){
MTri3 *t2 = t1->getNeigh(iLocalEdge);
if (!t2) return false;
@@ -149,6 +169,16 @@ bool gmshEdgeSwap(MTri3 *t1,
for (int i=0;i<3;i++)
if (t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
v4 = t2->tri()->getVertex(i);
+
+ // printf("%d %d %d %d %d\n",Us.size(),v1->getNum(),v2->getNum(),v3->getNum(),v4->getNum());
+
+ // printf("%d %d %d %d\n",tv1,tv2,tv3,tv4);
+
+ swapquad sq (v1,v2,v3,v4);
+ if (configs.find(sq) != configs.end())return false;
+ configs.insert(sq);
+
+ // if (tv1 != 0 && tv1 == tv2 && tv1 == tv3 && tv1 == tv4)return false;
const double volumeRef = surfaceTriangleUV (v1,v2,v3,Us,Vs) + surfaceTriangleUV (v1,v2,v4,Us,Vs);
@@ -182,8 +212,8 @@ bool gmshEdgeSwap(MTri3 *t1,
}
case SWCR_DEL:
{
- double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()])*.5,
- (Vs[v1->getNum()]+Vs[v2->getNum()])*.5};
+ double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()]+Us[v3->getNum()]+Us[v4->getNum()])*.25,
+ (Vs[v1->getNum()]+Vs[v2->getNum()]+Vs[v3->getNum()]+Vs[v4->getNum()])*.25};
double uv4[2] ={Us[v4->getNum()],Vs[v4->getNum()]};
double metric[3];
buildMetric ( gf , edgeCenter , metric);
@@ -255,6 +285,7 @@ bool gmshEdgeSwap(MTri3 *t1,
connectTriangles ( cavity );
newTris.push_back(t1b3);
newTris.push_back(t2b3);
+
return true;
}
@@ -572,6 +603,8 @@ bool gmshVertexCollapse(const double lMin,
return true;
}
+
+
int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
const gmshSwapCriterion &cr,
const std::vector<double> & Us ,
@@ -580,27 +613,34 @@ int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
const std::vector<double> & vSizesBGM)
{
typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
- std::vector<MTri3*> newTris;
- int nbSwap = 0;
-
- for (CONTAINER::iterator it = allTris.begin();it!=allTris.end();++it){
- if (!(*it)->isDeleted()){
- for (int i=0;i<3;i++){
- if (gmshEdgeSwap(*it,gf,i,newTris,cr,Us,Vs,vSizes,vSizesBGM)) {nbSwap++;break;}
+ int nbSwapTot=0;
+ std::set<swapquad> configs;
+ for (int iter=0;iter<1200;iter++){
+ int nbSwap = 0;
+ std::vector<MTri3*> newTris;
+ for (CONTAINER::iterator it = allTris.begin();it!=allTris.end();++it){
+ if (!(*it)->isDeleted()){
+ for (int i=0;i<3;i++){
+ if (gmshEdgeSwap(configs,*it,gf,i,newTris,cr,Us,Vs,vSizes,vSizesBGM)) {nbSwap++;break;}
+ }
+ }
+ else{
+ delete *it;
+ CONTAINER::iterator itb = it;
+ ++it;
+ allTris.erase(itb);
}
}
- else{
- CONTAINER::iterator itb = it;
- ++it;
- delete *itb;
- allTris.erase(itb);
- }
+ allTris.insert(newTris.begin(),newTris.end());
+ // printf("iter %d nbswam %d\n",iter,nbSwap);
+ nbSwapTot+=nbSwap;
+ if (nbSwap == 0)break;
}
+
// printf("B %d %d tris ",allTris.size(),newTris.size());
- allTris.insert(newTris.begin(),newTris.end());
// printf("A %d %d tris\n",allTris.size(),newTris.size());
- return nbSwap;
+ return nbSwapTot;
}
int edgeSplitPass (double maxLC,
--
GitLab