From 4fbec4e23ec897b69b78cc7b6eac7381a8c5d603 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Sun, 14 Oct 2007 17:30:42 +0000
Subject: [PATCH] *** empty log message ***
---
Mesh/BDS.cpp | 129 ++++++++++++++++++++++++++++++++++++++++++++-
Mesh/BDS.h | 1 +
Mesh/meshGFace.cpp | 87 +++++++++++++++++++++++++++---
3 files changed, 208 insertions(+), 9 deletions(-)
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index a682301007..ec9bb8a5ce 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.81 2007-10-11 14:34:04 remacle Exp $
+// $Id: BDS.cpp,v 1.82 2007-10-14 17:30:42 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -645,6 +645,131 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
}
+void BDS_Mesh::saturate_edge(BDS_Edge * e, std::vector<BDS_Point *> &mids)
+{
+ BDS_Point *op[2];
+ BDS_Point *p1 = e->p1;
+ BDS_Point *p2 = e->p2;
+
+ BDS_Point *pts1[4];
+ e->faces(0)->getNodes(pts1);
+
+ int orientation = 0;
+ for(int i = 0; i < 3; i++) {
+ if(pts1[i] == p1) {
+ if(pts1[(i + 1) % 3] == p2)
+ orientation = 1;
+ else
+ orientation = -1;
+
+ break;
+ }
+ }
+
+ // printf("sat %d\n",mids.size());
+
+ // we should project
+ e->oppositeof(op);
+ BDS_GeomEntity *g1 = 0, *g2 = 0, *ge = e->g;
+
+ BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
+ BDS_Edge *op1_p2 = find_edge(op[0], p2, e->faces(0));
+ BDS_Edge *p1_op2 = find_edge(p1, op[1], e->faces(1));
+ BDS_Edge *op2_p2 = find_edge(op[1], p2, e->faces(1));
+
+ if(e->faces(0)) {
+ g1 = e->faces(0)->g;
+ del_face(e->faces(0));
+ }
+ // not a bug !!!
+ if(e->faces(0)) {
+ g2 = e->faces(0)->g;
+ del_face(e->faces(0));
+ }
+
+ // double t_l = e->target_length;
+
+ del_edge(e);
+
+ // create all the sub-edges of e
+ std::vector<BDS_Edge*> subs;
+ for (int i=0; i<mids.size()+1; i++)
+ {
+ BDS_Point *a,*b;
+ if (i == 0)a = p1;
+ else a = mids[i-1];
+ if (i == mids.size())b = p2;
+ else b = mids[i];
+ BDS_Edge *sub = new BDS_Edge(a, b);
+ // printf("%g %g -> %g %g\n",a->X,a->Y,b->X,b->Y);
+ edges.push_back(sub);
+ subs.push_back(sub);
+ sub->g = ge;
+ }
+ // create edges that connect op1 and op2
+ std::vector<BDS_Edge*> conn1,conn2;
+ for (int i=0;i<mids.size();i++)
+ {
+ BDS_Edge *c1 = new BDS_Edge(mids[i], op[0]);
+ edges.push_back(c1);
+ conn1.push_back(c1);
+ BDS_Edge *c2 = new BDS_Edge(mids[i], op[1]);
+ edges.push_back(c2);
+ conn2.push_back(c2);
+ c1->g = g1;
+ c2->g = g2;
+ mids[i]->g = ge;
+ }
+
+ // create the triangles
+
+ for (int i=0;i<mids.size()+1;i++)
+ {
+ BDS_Edge *e1,*e2,*e3=subs[i];
+
+ // printf("--> %d %g %g ->%d %g %g\n",e3->p1->iD,e3->p1->X,e3->p1->Y,e3->p2->iD,e3->p2->X,e3->p2->Y);
+
+
+ if (i==0)e1 = p1_op1;
+ else e1 = conn1[i-1];
+ if (i==mids.size())e2 = op1_p2;
+ else e2 = conn1[i];
+
+// printf("--> %d %g %g ->%d %g %g\n",e1->p1->iD,e1->p1->X,e1->p1->Y,e1->p2->iD,e1->p2->X,e1->p2->Y);
+// printf("--> %d %g %g ->%d %g %g\n",e2->p1->iD,e2->p1->X,e2->p1->Y,e2->p2->iD,e2->p2->X,e2->p2->Y);
+
+ BDS_Face *t1;
+ if(orientation == 1)
+ t1 = new BDS_Face(e3, e2, e1);
+ else
+ t1 = new BDS_Face(e3, e1, e2);
+
+
+ if (i==0)e1 = p1_op2;
+ else e1 = conn2[i-1];
+ if (i==mids.size())e2 = op2_p2;
+ else e2 = conn2[i];
+
+ BDS_Face *t2;
+ if(orientation != 1)
+ t2 = new BDS_Face(e3, e2, e1);
+ else
+ t2 = new BDS_Face(e3, e1, e2);
+ t1->g = g1;
+ t2->g = g2;
+
+ triangles.push_back(t1);
+ triangles.push_back(t2);
+ }
+ // config has changed
+ p1->config_modified = true;
+ p2->config_modified = true;
+ op[0]->config_modified = true;
+ op[1]->config_modified = true;
+}
+
+
+
// This function does actually the swap without taking into account
// the feasability of the operation. Those conditions have to be
@@ -682,6 +807,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
// we test if the edge is deleted
// return false;
+
if(e->deleted)
return false;
@@ -692,6 +818,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
if(e->g && e->g->classif_degree == 1)
return false;
+
BDS_Point *op[2];
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index eabd18d6f8..d3b7c74310 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -457,6 +457,7 @@ public:
bool smooth_point_centroid(BDS_Point * p, GFace *gf);
bool move_point(BDS_Point *p , double X, double Y, double Z);
void split_edge(BDS_Edge *, BDS_Point *);
+ void saturate_edge(BDS_Edge *, std::vector<BDS_Point *>&);
bool edge_constraint ( BDS_Point *p1, BDS_Point *p2 );
bool recombine_edge ( BDS_Edge *e );
// Global operators
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e819f854b4..b9dce76f47 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.96 2007-10-11 16:08:23 remacle Exp $
+// $Id: meshGFace.cpp,v 1.97 2007-10-14 17:30:42 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -255,6 +255,14 @@ double NewGetLc(BDS_Edge *e, GFace *f)
return 2*linearLength / (l1 + l2);
}
+double NewGetLc(BDS_Point *p1,BDS_Point *p2, GFace *f)
+{
+ double linearLength = computeEdgeLinearLength(p1,p2,f);
+ double l1 = NewGetLc(p1);
+ double l2 = NewGetLc(p2);
+ return 2*linearLength / (l1 + l2);
+}
+
bool edgeSwapTest(BDS_Edge *e,GFace *gf)
{
BDS_Point *op[2];
@@ -318,9 +326,15 @@ bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
double p2x[3] = {e->p2->X,e->p2->Y,e->p2->Z};
double op1x[3] = {op[0]->X,op[0]->Y,op[0]->Z};
double op2x[3] = {op[1]->X,op[1]->Y,op[1]->Z};
+ double p1u[2] = {e->p1->u,e->p1->v};
+ double p2u[2] = {e->p2->u,e->p2->v};
+ double op1u[2] = {op[0]->u,op[0]->v};
+ double op2u[2] = {op[1]->u,op[1]->v};
double fourth[3];
fourthPoint(p1x,p2x,op1x,fourth);
double result = gmsh::insphere(p1x, p2x, op1x, fourth, op2x) * gmsh::orient3d(p1x, p2x, op1x, fourth);
+ //double result = gmsh::incircle(p1u, p2u, op1u, op2u) * gmsh::orient2d(p1u, p2u, op1u);
+ // printf("result = a%12.5E\n",result);
return result > 0.;
}
@@ -404,11 +418,13 @@ void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
{
int result = edgeSwapTestQuality(*it,5);
if (result >= 0)
- if(edgeSwapTestDelaunay(*it,gf) || result > 0)
- if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
- nb_swap++;
- ++it;
+ {
+ if(edgeSwapTestDelaunay(*it,gf) || result > 0)
+ if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
+ nb_swap++;
+ }
}
+ ++it;
}
}
@@ -425,13 +441,18 @@ void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
double lone = NewGetLc ( *it,gf);
if ((*it)->numfaces() == 2 && (lone > MAXE_))
{
+
+// BDS_Point *op[2];
+// (*it)->oppositeof(op);
+// double lone1 = NewGetLc ( op[0],mid,gf);
+// double lone2 = NewGetLc ( op[1],mid,gf);
+
const double coord = 0.5;
BDS_Point *mid ;
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v,gf);
- double l1;
- // if (BGMExists())
+
mid->lcBGM() = BGM_MeshSize(gf,
(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
(coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
@@ -440,6 +461,52 @@ void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
mid->lc() = 0.5 * ( (*it)->p1->lc() + (*it)->p2->lc() );
m.split_edge ( *it, mid );
nb_split++;
+ }
+ }
+ ++it;
+ }
+}
+
+void saturateEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
+{
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1)
+ {
+ if (NN2++ >= NN1)break;
+ if (!(*it)->deleted)
+ {
+ double lone = NewGetLc ( *it,gf);
+ if ((*it)->numfaces() == 2 && (lone > MAXE_))
+ {
+ int nbSub = (int) (lone/MAXE_) ;
+ // nbSub = std::min(nbSub,2);
+ // printf("%d %g\n",nbSub,lone/MAXE_);
+ std::vector<BDS_Point*> mids;
+ for (int i=0;i<nbSub;i++)
+ {
+ const double coord = (double)(i+1)/(nbSub+1);
+ BDS_Point *mid ;
+ mid = m.add_point(++m.MAXPOINTNUMBER,
+ (1.-coord) * (*it)->p1->u + coord * (*it)->p2->u,
+ (1.-coord) * (*it)->p1->v + coord * (*it)->p2->v,gf);
+ double l1;
+ // if (BGMExists())
+ mid->lcBGM() = BGM_MeshSize(gf,
+ ((1.-coord) * (*it)->p1->u + (coord) * (*it)->p2->u)*m.scalingU,
+ ((1.-coord) * (*it)->p1->v + (coord) * (*it)->p2->v)*m.scalingV,
+ mid->X,mid->Y,mid->Z);
+ //mid->lc() = 2./ ( 1./(*it)->p1->lc() + 1./(*it)->p2->lc() );
+ mid->lc() = ( (1.-coord)*(*it)->p1->lc() + coord*(*it)->p2->lc() );
+ mids.push_back(mid);
+ // printf("new point %g %g lc %g\n",mid->X,mid->Y,mid->lc());
+ }
+ // printf("saturating an edge with %d points %d triangles\n",mids.size(),m.triangles.size());
+ if(nbSub>0)m.saturate_edge ( *it, mids );
+ // printf("-> %d triangles\n",m.triangles.size());
+ nb_split++;
+ // if (nb_split == )break;
}
}
++it;
@@ -518,7 +585,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
double t_spl=0, t_sw=0,t_col=0,t_sm=0;
- const double MINE_ = 0.7, MAXE_=1.4;
+ const double MINE_ = 0.67, MAXE_=1.4;
while (1)
{
// we count the number of local mesh modifs.
@@ -558,6 +625,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
double minE = MINE_;//std::min(MINE_,minL * 1.2);
clock_t t1 = clock();
splitEdgePass ( gf, m, maxE, nb_split);
+ //saturateEdgePass ( gf, m, maxE, nb_split);
clock_t t2 = clock();
swapEdgePass ( gf, m, nb_swap);
clock_t t3 = clock();
@@ -567,12 +635,15 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
clock_t t5 = clock();
smoothVertexPass ( gf, m, nb_smooth);
clock_t t6 = clock();
+ swapEdgePass ( gf, m, nb_swap);
+ clock_t t7 = clock();
// clean up the mesh
t_spl += (double)(t2-t1)/CLOCKS_PER_SEC;
t_sw += (double)(t3-t2)/CLOCKS_PER_SEC;
t_sw += (double)(t5-t4)/CLOCKS_PER_SEC;
+ t_sw += (double)(t7-t6)/CLOCKS_PER_SEC;
t_col += (double)(t4-t3)/CLOCKS_PER_SEC;
t_sm += (double)(t6-t5)/CLOCKS_PER_SEC;
--
GitLab