diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 310b1b7d48248a6b64efc66afde28c876e819439..133d06e0608b966583999cac498a0304e9f4a71f 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.57 2006-08-19 20:51:44 geuzaine Exp $
+// $Id: BDS.cpp,v 1.58 2006-08-21 13:32:42 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -470,7 +470,7 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
rhs[1] = y3 - y1;
if(!sys2x2(mat, rhs, x))
return 0;
- if(x[0] > 0.0 && x[0] < 1.0 && x[1] > 0.0 && x[1] < 1.0)
+ if(x[0] >= 0.0 && x[0] <= 1.0 && x[1] >= 0.0 && x[1] <= 1.0)
return 1;
return 0;
}
@@ -496,6 +496,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2)
while (it!=edges.end())
{
e = (*it);
+ // if (e->p1->iD >= 0 && e->p2->iD >= 0)Msg(INFO," %d %d %g %g - %g %g",e->p1->iD,e->p2->iD,e->p1->u,e->p1->v,e->p2->u,e->p2->v);
if (!e->deleted && e->p1 != p1 && e->p1 != p2 && e->p2 != p1 && e->p2 != p2)
if(Intersect_Edges_2d(e->p1->X, e->p1->Y,
e->p2->X, e->p2->Y,
@@ -1586,6 +1587,12 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
triangles.push_back(t2);
triangles.push_back(t3);
triangles.push_back(t4);
+
+ // config has changed
+ p1->config_modified = true;
+ p2->config_modified = true;
+ op[0]->config_modified = true;
+ op[1]->config_modified = true;
}
bool BDS_SwapEdgeTestNonPlanar::operator () (BDS_Edge *e,
@@ -1672,7 +1679,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
// return false;
if(e->deleted)
return false;
-
+
int nbFaces = e->numfaces();
if(nbFaces != 2)
return false;
@@ -1750,11 +1757,17 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
triangles.push_back(t1);
triangles.push_back(t2);
+
+ p1->config_modified = true;
+ p2->config_modified = true;
+ op[0]->config_modified = true;
+ op[1]->config_modified = true;
+
return true;
}
-bool BDS_Mesh::collapse_edge(BDS_Edge * e, BDS_Point * p, const double eps)
+bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
{
if(e->numfaces() != 2)
@@ -1827,6 +1840,7 @@ bool BDS_Mesh::collapse_edge(BDS_Edge * e, BDS_Point * p, const double eps)
std::list < BDS_Edge * >::iterator eit = edges.begin();
std::list < BDS_Edge * >::iterator eite = edges.end();
while(eit != eite) {
+ (*eit)->p1->config_modified = (*eit)->p2->config_modified = true;
ept[0][kk] = ((*eit)->p1 == p) ? o->iD : (*eit)->p1->iD;
ept[1][kk] = ((*eit)->p2 == p) ? o->iD : (*eit)->p2->iD;
egs[kk++] = (*eit)->g;
@@ -1994,6 +2008,8 @@ void BDS_Mesh::snap_point(BDS_Point * p, BDS_Mesh * geom_mesh)
bool BDS_Mesh::smooth_point(BDS_Point * p, BDS_Mesh * geom_mesh)
{
+
+
if(p->g && p->g->classif_degree <= 1)
return false;
// if(!p->g->surf)
@@ -2048,11 +2064,13 @@ bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS
bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
{
+ if (!p->config_modified)return false;
if(p->g && p->g->classif_degree <= 1)
return false;
double U = 0;
double V = 0;
+ double LC = 0;
std::list < BDS_Edge * >::iterator eit = p->edges.begin();
while(eit != p->edges.end()) {
@@ -2060,11 +2078,13 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
U += op->u;
V += op->v;
+ LC += op->lc();
++eit;
}
U /= p->edges.size();
V /= p->edges.size();
+ LC /= p->edges.size();
std::list < BDS_Triangle * >ts;
p->getTriangles(ts);
@@ -2081,6 +2101,7 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
GPoint gp = gf->point(U,V);
p->u = U;
p->v = V;
+ p->lc() = LC;
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
@@ -2198,23 +2219,6 @@ void BDS_Mesh::compute_metric_edge_lengths(const BDS_Metric & metric)
// printf("end computation of metric lengths\n");
}
-
-void BDS_Mesh::applyOptimizationPatterns()
-{
- std::set<BDS_Point*, PointLessThan>::iterator it = points.begin();
- std::set<BDS_Point*, PointLessThan>::iterator ite = points.end();
- while (it != ite)
- {
- if ((*it)->edges.size() == 3 && (*it)->g->classif_degree == 2)
- {
- BDS_Edge *e = *((*it)->edges.begin());
- if (!e->deleted)collapse_edge ( e , *it , .1);
- }
- ++it;
- }
-}
-
-
int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
BDS_Mesh * geom_mesh)
{
@@ -2265,16 +2269,16 @@ int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
{
double length = (*it)->length();
if (!(*it)->deleted && length < (*it)->target_length * 0.7 ){
- if (nb_modif %2 )
- {
- if (collapse_edge (*it, (*it)->p1, 0.1 ))
- nb_modif++;
- }
- else
- {
- if (collapse_edge (*it, (*it)->p2, 0.1 ))
- nb_modif++;
- }
+// if (nb_modif %2 )
+// {
+// if (collapse_edge (*it, (*it)->p1, 0.1 ))
+// nb_modif++;
+// }
+// else
+// {
+// if (collapse_edge (*it, (*it)->p2, 0.1 ))
+// nb_modif++;
+// }
}
++it;
}
@@ -2332,7 +2336,6 @@ int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
Msg(INFO,"%d snaps have succeeded , %d have failed\n",SNAP_SUCCESS,SNAP_FAILURE);
// outputScalarField (triangles,"b.pos");
- applyOptimizationPatterns(); // FIXME: this is buggy
cleanup();
return nb_modif;
}
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index e331480e049ce3c2d704d044bdd8f9a66badcdd3..4d1a462381ebc8c8279fed650b444e908374e18c 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -227,6 +227,7 @@ class BDS_Point
public:
double X,Y,Z; // Real COORDINATES
double u,v; // Parametric COORDINATES
+ bool config_modified;
int iD;
BDS_GeomEntity *g;
std::list<BDS_Edge*> edges;
@@ -257,7 +258,7 @@ public:
void getTriangles(std::list<BDS_Triangle *> &t) const;
void compute_curvature();
BDS_Point(int id, double x=0, double y=0, double z=0)
- : radius_of_curvature(1.e22),X(x),Y(y),Z(z),u(0),v(0),iD(id),g(0)
+ : radius_of_curvature(1.e22),X(x),Y(y),Z(z),u(0),v(0),config_modified(true),iD(id),g(0)
{
}
};
@@ -702,13 +703,12 @@ public:
// 2D operators
BDS_Edge *recover_edge(int p1, int p2);
bool swap_edge(BDS_Edge *, const BDS_SwapEdgeTest &theTest);
- bool collapse_edge(BDS_Edge *, BDS_Point*, const double eps);
+ bool collapse_edge_parametric(BDS_Edge *, BDS_Point*);
void snap_point(BDS_Point* , BDS_Mesh *geom = 0);
bool smooth_point(BDS_Point* , BDS_Mesh *geom = 0);
bool smooth_point_parametric(BDS_Point * p, GFace *gf);
bool move_point(BDS_Point *p , double X, double Y, double Z);
void split_edge(BDS_Edge *, BDS_Point *);
- bool delaunay_insertion ( int num , double x, double y,double z );
bool edge_constraint ( BDS_Point *p1, BDS_Point *p2 );
// Global operators
void classify(double angle, int nb = -1);
@@ -724,6 +724,7 @@ public:
void applyOptimizationPatterns();
};
+void outputScalarField(std::list < BDS_Triangle * >t, const char *fn);
void recur_tag(BDS_Triangle * t, BDS_GeomEntity * g);
bool project_point_on_a_list_of_triangles(BDS_Point *p , const std::list<BDS_Triangle*> &t,
double &X, double &Y, double &Z);
diff --git a/Mesh/Utils.cpp b/Mesh/Utils.cpp
index 477957c576e609d339a640837d132446ebc96412..b0cef3eeb6b0fcacf925af5bbdfb7129d5dd3474 100644
--- a/Mesh/Utils.cpp
+++ b/Mesh/Utils.cpp
@@ -1,4 +1,4 @@
-// $Id: Utils.cpp,v 1.32 2006-07-25 12:08:24 remacle Exp $
+// $Id: Utils.cpp,v 1.33 2006-08-21 13:32:42 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -418,13 +418,17 @@ void XYZtoUV(Surface * s, double X, double Y, double Z, double *U, double *V,
(jac[0][1] * (X - P.Pos.X) + jac[1][1] * (Y - P.Pos.Y) +
jac[2][1] * (Z - P.Pos.Z));
- err = DSQR(Unew - *U) + DSQR(Vnew - *V);
+ // err = DSQR(Unew - *U) + DSQR(Vnew - *V);
+ // A BETTER TEST !! (JFR/AUG 2006)
+ err = DSQR(X - P.Pos.X) + DSQR(Y - P.Pos.Y) + + DSQR(Z - P.Pos.Z);
iter++;
*U = Unew;
*V = Vnew;
}
+
+
if(iter < MaxIter && err <= Precision &&
Unew <= umax && Vnew <= vmax &&
Unew >= umin && Vnew >= vmin){
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 5d78a080718badfc2243bfbf4bcaab0641105384..173618737aaedebfbed164b31d78f1c117842164 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.12 2006-08-19 08:26:47 remacle Exp $
+// $Id: meshGEdge.cpp,v 1.13 2006-08-21 13:32:42 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -49,12 +49,12 @@ double F_Lc_bis(double t)
const double fact = (t-t_begin)/(t_end-t_begin);
double lc_here = lc_begin + fact * (lc_end-lc_begin);
SVector3 der = _myGEdge -> firstDer(t) ;
- const double d = norm(der);
+ const double d = norm(der)/CTX.mesh.lc_factor;
// TESTTTT
GPoint points = _myGEdge -> point(t) ;
- double l_bgm = F_LC_ANALY(points.x(),points.y(),points.z()) ;
- lc_here = l_bgm;
+ // double l_bgm = F_LC_ANALY(points.x(),points.y(),points.z()) ;
+ // lc_here = l_bgm;
if(CTX.mesh.bgmesh_type == ONFILE) {
GPoint point = _myGEdge -> point(t) ;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 5e10deb6de5655792484362090ae257aa61f253f..f8e14e3ec0ea1bb60067591e335197ab084023fc 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.6 2006-08-19 08:26:47 remacle Exp $
+// $Id: meshGFace.cpp,v 1.7 2006-08-21 13:32:42 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -219,6 +219,7 @@ extern double F_LC_ANALY (double xx, double yy, double zz);
double NewGetLc ( BDS_Point *p )
{
+ return p->lc();
return F_LC_ANALY(p->X,p->Y,p->Z);
if(BGMExists())
{
@@ -249,6 +250,8 @@ bool edgeSwapTest(BDS_Edge *e)
{
BDS_Point *op[2];
+ if (!e->p1->config_modified && ! e->p2->config_modified) return false;
+
if (e->numfaces() != 2)return false;
e->oppositeof (op);
@@ -287,32 +290,94 @@ bool edgeSwapTest(BDS_Edge *e)
}
-void RefineMesh ( GFace *gf, BDS_Mesh &m )
+void OptimizeMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
+{
+ // optimize
+ for(int i = 0 ; i < NIT ; i++){
+ {
+ std::set<BDS_Point*,PointLessThan> PTS (m.points);
+ std::set<BDS_Point*,PointLessThan>::iterator itp = PTS.begin();
+ while (itp != PTS.end())
+ {
+ std::list < BDS_Triangle * >t;
+ (*itp)->getTriangles(t);
+ if ((t.size()==3 && (*itp)->edges.size() == 3)||
+ (t.size()==4 && (*itp)->edges.size() == 4))
+ m.collapse_edge_parametric ( *(*itp)->edges.begin(), (*itp));
+ else
+ m.smooth_point_parametric(*itp,gf);
+ ++itp;
+ }
+ }
+ {
+ // swap edges that provide a better configuration
+ std::list<BDS_Edge *> temp (m.edges);
+ std::list<BDS_Edge*>::iterator it = temp.begin();
+ while (it != temp.end())
+ {
+ if (!(*it)->deleted && (*it)->numfaces() == 2)
+ {
+ BDS_Point *op[2];
+ (*it)->oppositeof(op);
+
+ bool c1 = (op[0]-> edges.size() == 5 && op[1]-> edges.size() == 5 );
+ bool c2 = ((*it)->p1-> edges.size() == 7 && (*it)->p2-> edges.size() == 7);
+ if ( c1 && c2 )
+ m.swap_edge ( *it , BDS_SwapEdgeTestParametric());
+ }
+ ++it;
+ }
+ }
+ m.cleanup();
+ }
+}
+
+void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
{
int NUMP = 0;
int IT =0;
int MAXNP = m.MAXPOINTNUMBER;
+ clock_t t1 = clock();
+
+ // computecharacteristic lengths using mesh edge spacing
+
+ std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
+ while (itp != m.points.end())
+ {
+ std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
+ std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
+ double L = 1.e245;
+ while(it!=ite){
+ double l = (*it)->length();
+ if (l<L && (*it)->g && (*it)->g->classif_degree == 1)L=l;
+ ++it;
+ }
+ (*itp)->lc() = L;
+ ++itp;
+ }
+
+
while (1)
{
- double stot = 0;
- std::list<BDS_Triangle *>::iterator ittt = m.triangles.begin();
- while (ittt!= m.triangles.end())
- {
- if (!(*ittt)->deleted)
- {
- BDS_Point *pts[3];
- (*ittt)->getNodes(pts);
- stot += fabs( surface_triangle_param(pts[0], pts[1], pts[2]));
- }
- ++ittt;
- }
+ // double stot = 0;
+// std::list<BDS_Triangle *>::iterator ittt = m.triangles.begin();
+// while (ittt!= m.triangles.end())
+// {
+// if (!(*ittt)->deleted)
+// {
+// BDS_Point *pts[3];
+// (*ittt)->getNodes(pts);
+// stot += fabs( surface_triangle_param(pts[0], pts[1], pts[2]));
+// }
+// ++ittt;
+// }
std::list<BDS_Edge *> temp (m.edges);
std::list<BDS_Edge*>::iterator it = temp.begin();
int nb_split=0;
- int nb_short=0;
+ int nb_smooth=0;
int nb_collaps=0;
int nb_swap=0;
@@ -323,6 +388,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m )
{
if (!(*it)->deleted)
{
+ (*it)->p1->config_modified = false;
+ (*it)->p2->config_modified = false;
double lone = NewGetLc ( *it);
maxL = std::max(maxL,lone);
minL = std::min(minL,lone);
@@ -330,8 +397,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m )
++it;
}
- Msg(INFO,"stot %22.15E minL %g maxL %g",stot,minL,maxL);
- if ((minL > 0.2 && maxL < 1.3) || IT > 7)break;
+ Msg(STATUS2," %d triangles : conv %g -> %g (%g sec)",m.triangles.size(),maxL,1.5,(double)(clock()-t1)/CLOCKS_PER_SEC);
+ if ((minL > 0.2 && maxL < 1.5) || IT > NIT)break;
it = temp.begin();
@@ -340,13 +407,14 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m )
if (!(*it)->deleted)
{
double lone = NewGetLc ( *it);
- if ((*it)->numfaces() == 2 && (lone > 1.4))
+ if ((*it)->numfaces() == 2 && (lone > 1.3))
{
BDS_Point *mid ;
double coord = 0.5;
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);
+ mid->lc() = 0.5 * ( (*it)->p1->lc() + (*it)->p2->lc() );
m.split_edge ( *it, mid );
nb_split++;
}
@@ -354,24 +422,32 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m )
++it;
}
+ // swap edges that provide a better configuration
+ temp = m.edges;
+ it = temp.begin();
+ while (it != temp.end())
+ {
+ if (!(*it)->deleted && edgeSwapTest(*it))
+ if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
+ nb_swap++;
+ ++it;
+ }
+
// collapse short edges
temp = m.edges;
it = temp.begin();
while (it != temp.end())
{
double lone = NewGetLc ( *it);
- if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < 0.7 )
+ if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < 0.6 )
{
bool res = false;
if ( (*it)->p1->iD > MAXNP )
- res =m.collapse_edge ( *it, (*it)->p1, 0);
+ res =m.collapse_edge_parametric ( *it, (*it)->p1);
else if ( (*it)->p2->iD > MAXNP )
- res =m.collapse_edge ( *it, (*it)->p2, 0);
-
+ res =m.collapse_edge_parametric ( *it, (*it)->p2);
if (res)
nb_collaps ++;
- else
- nb_short++;
}
++it;
}
@@ -387,35 +463,21 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m )
++it;
}
// smooth resulting mesh
- // if (IT % 4 == 0 )
+ //if (IT % 4 == 0 )
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
{
- m.smooth_point_parametric(*itp,gf);
+
+ if(m.smooth_point_parametric(*itp,gf))
+ nb_smooth ++;
++itp;
}
}
- IT++;
- m.cleanup();
- }
-
- // optimize
- {
- std::set<BDS_Point*,PointLessThan> PTS (m.points);
- std::set<BDS_Point*,PointLessThan>::iterator itp = PTS.begin();
- while (itp != PTS.end())
- {
- std::list < BDS_Triangle * >t;
- (*itp)->getTriangles(t);
- if ((t.size()==3 && (*itp)->edges.size() == 3)||
- (t.size()==4 && (*itp)->edges.size() == 4))
- m.collapse_edge ( *(*itp)->edges.begin(), (*itp), 0);
- else
- m.smooth_point_parametric(*itp,gf);
- ++itp;
- }
- }
+ IT++;
+ // Msg(INFO," %d split %d swap %d collapse %d smooth",nb_split,nb_swap,nb_collaps,nb_smooth);
+ m.cleanup();
+ }
}
void OldMeshGenerator ( GFace *gf,
@@ -542,7 +604,7 @@ void OldMeshGenerator ( GFace *gf,
}
-void recover_medge ( BDS_Mesh *m, GEdge *ge)
+bool recover_medge ( BDS_Mesh *m, GEdge *ge)
{
m->add_geom (ge->tag(), 1);
@@ -550,10 +612,12 @@ void recover_medge ( BDS_Mesh *m, GEdge *ge)
if (ge->mesh_vertices.size() == 0)
{
- MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
- MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
+
+ MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
+ MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
BDS_Point *pstart = m->find_point(vstart->getNum());
- BDS_Point *pend = m->find_point(vend->getNum());
+ BDS_Point *pend = m->find_point(vend->getNum());
+
if(!pstart->g)
{
m->add_geom (vstart->getNum(), 0);
@@ -568,8 +632,12 @@ void recover_medge ( BDS_Mesh *m, GEdge *ge)
}
BDS_Edge * e = m->recover_edge ( vstart->getNum(), vend->getNum());
if (e)e->g = g;
- else throw;
- return;
+ else {
+ Msg(GERROR,"Unable to recover an edge %g %g && %g %g (size 0)",vstart->x(),vstart->y(), vend->x(),vend->y());
+ Msg(GERROR,"Mesh with %d triangles",m->triangles.size());
+ return false;
+ }
+ return true;
}
@@ -594,13 +662,19 @@ void recover_medge ( BDS_Mesh *m, GEdge *ge)
vend = ge->mesh_vertices[i];
e = m->recover_edge ( vstart->getNum(), vend->getNum());
if (e)e->g = g;
- else throw;
+ else {
+ Msg(GERROR,"Unable to recover an edge %g %g && %g %g (%d/*d)",vstart->x(),vstart->y(), vend->x(),vend->y(),i,ge->mesh_vertices.size());
+ return false;
+ }
}
vstart = vend;
vend = *(ge->getEndVertex()->mesh_vertices.begin());
e = m->recover_edge ( vstart->getNum(), vend->getNum());
if (e)e->g = g;
- else throw;
+ else {
+ Msg(GERROR,"Unable to recover an edge %g %g && %g %g",vstart->x(),vstart->y(), vend->x(),vend->y());
+ return false;
+ }
BDS_Point *pend = m->find_point(vend->getNum());
if(!pend->g)
@@ -609,7 +683,7 @@ void recover_medge ( BDS_Mesh *m, GEdge *ge)
BDS_GeomEntity *g0 = m->get_geom(vend->getNum(), 0);
pend->g = g0;
}
-
+ return true;
}
@@ -620,6 +694,8 @@ void recover_medge ( BDS_Mesh *m, GEdge *ge)
void gmsh2DMeshGenerator ( GFace *gf )
{
+ // if(gf->tag() == 138|| gf->tag() == 196|| gf->tag() == 200|| gf->tag() == 262) return;
+
std::set<MVertex*>all_vertices;
std::map<int, MVertex*>numbered_vertices;
std::list<GEdge*> edges = gf->edges();
@@ -689,7 +765,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
}
/// Increase the size of the bounding box by 20 %
- bbox *= 1.2;
+ bbox *= 1.5;
/// add 4 points than encloses the domain
/// Use negative number to distinguish thos fake vertices
MVertex *bb[4];
@@ -734,7 +810,10 @@ void gmsh2DMeshGenerator ( GFace *gf )
for ( int ip = 0 ; ip<4 ; ip++ ) delete bb[ip];
// Recover the boundary edges
+ // and compute characteristic lenghts using mesh edge spacing
+ BDS_GeomEntity CLASS_F (1,2);
+
it = edges.begin();
while(it != edges.end())
{
@@ -748,14 +827,12 @@ void gmsh2DMeshGenerator ( GFace *gf )
++it;
}
- Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());
-
+ // Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());
// Look for an edge that is on the boundary for which one of the
// two neighbors has a negative number node. The other triangle
// is inside the domain and, because all edges were recovered,
// triangles inside the domain can be recovered using a simple
// recursive algorithm
- BDS_GeomEntity CLASS_F (1,2);
{
std::list<BDS_Edge*>::iterator ite = m->edges.begin();
while (ite != m->edges.end())
@@ -814,9 +891,11 @@ void gmsh2DMeshGenerator ( GFace *gf )
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
- // start mesh generation
+ // goto hhh;
+ // start mesh generation
- RefineMesh (gf,*m);
+ RefineMesh (gf,*m,27);
+ OptimizeMesh (gf,*m,2);
// fill the small gmsh structures
@@ -827,7 +906,8 @@ void gmsh2DMeshGenerator ( GFace *gf )
BDS_Point *p = *itp;
if (numbered_vertices.find(p->iD) == numbered_vertices.end())
{
- MVertex *v = new MFaceVertex (p->X,p->Y,p->Z,gf,p->u,p->v);
+ //MVertex *v = new MFaceVertex (p->X,p->Y,p->Z,gf,p->u,p->v);
+ MVertex *v = new MFaceVertex (p->u,p->v,0,gf,p->u,p->v);
numbered_vertices[p->iD]=v;
gf->mesh_vertices.push_back(v);
}
@@ -854,13 +934,13 @@ void gmsh2DMeshGenerator ( GFace *gf )
// delete the mesh
+ // outputScalarField(m->triangles, "lc.pos");
delete m;
+
fromParametricToCartesian p2c ( gf );
std::for_each(all_vertices.begin(),all_vertices.end(),p2c);
- // std::for_each(gf->mesh_vertices.begin(),gf->mesh_vertices.end(),p2c);
-
-
+ std::for_each(gf->mesh_vertices.begin(),gf->mesh_vertices.end(),p2c);
}
@@ -881,7 +961,7 @@ void meshGFace :: operator() (GFace *gf)
if(gf->geomType() == GEntity::DiscreteSurface) return;
// Send a messsage to the GMSH environment
- Msg(INFO, "Meshing surface %d", gf->tag());
+ Msg(STATUS2, "Meshing surface %d", gf->tag());
// destroy the mesh if it exists
deMeshGFace dem;
@@ -902,10 +982,13 @@ void meshGFace :: operator() (GFace *gf)
Msg(DEBUG1, "points were put in parametric coords ...");
-
- gmsh2DMeshGenerator ( gf ) ;
-
-
+ try
+ {
+ gmsh2DMeshGenerator ( gf ) ;
+ }
+ catch (...)
+ {
+ }
Msg(DEBUG1, "type %d %d triangles generated, %d internal vertices",
gf->geomType(),gf->triangles.size(),gf->mesh_vertices.size());
diff --git a/benchmarks/2d/Square-Emb.geo b/benchmarks/2d/Square-Emb.geo
index 602285056ad97992448ed5fcdc8db063d8d01e2f..74aef27ddc5cd0e266f853374cc3463a6e8951b6 100644
--- a/benchmarks/2d/Square-Emb.geo
+++ b/benchmarks/2d/Square-Emb.geo
@@ -1,16 +1,17 @@
/******************************
Square uniformly meshed
******************************/
-lc = .49999;
+lc = .1;
+lc2 = .01;
Point(1) = {0.0,0.0,0,lc};
Point(2) = {1,0.0,0,lc};
Point(3) = {1,1,0,lc};
Point(4) = {0,1,0,lc};
-Point(11) = {0.2,0.2,0,lc};
-Point(22) = {0.3,0.8,0,lc};
-Point(33) = {0.8,0.8,0,lc};
-Point(44) = {0.9,0.2,0,lc};
+Point(11) = {0.2,0.2,0,lc2};
+Point(22) = {0.3,0.8,0,lc2};
+Point(33) = {0.8,0.8,0,lc2};
+Point(44) = {0.9,0.2,0,lc2};
Line(1) = {3,2};
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index dea5672ff84b7e2b976ef5280d2cdde9c4578a28..43358fce04fd05b4c2090d7337cd935d1d77b610 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -1,4 +1,6 @@
-lc = 0.05;
+lc = 0.001;
+lc2 = 0.2;
+lc3 = 0.05;
Point(1) = {1.000000e+00,0.000000e+00,0.000000e+00,lc};
Point(2) = {9.997533e-01,0.000000e+00,-3.498543e-05,lc};
Point(3) = {9.990134e-01,0.000000e+00,-1.398841e-04,lc};
@@ -48,8 +50,8 @@ Point(46) = {5.782179e-01,0.000000e+00,-4.655984e-02,lc};
Point(47) = {5.626673e-01,0.000000e+00,-4.777199e-02,lc};
Point(48) = {5.470549e-01,0.000000e+00,-4.894463e-02,lc};
Point(49) = {5.313960e-01,0.000000e+00,-5.007425e-02,lc};
-Point(50) = {5.157061e-01,0.000000e+00,-5.115728e-02,lc};
-Point(51) = {5.000008e-01,0.000000e+00,-5.219014e-02,lc};
+Point(50) = {5.157061e-01,0.000000e+00,-5.115728e-02,lc3};
+Point(51) = {5.000008e-01,0.000000e+00,-5.219014e-02,lc3};
Point(52) = {4.842954e-01,0.000000e+00,-5.316926e-02,lc};
Point(53) = {4.686055e-01,0.000000e+00,-5.409108e-02,lc};
Point(54) = {4.529467e-01,0.000000e+00,-5.495201e-02,lc};
@@ -149,7 +151,7 @@ Point(147) = {4.373342e-01,0.000000e+00,5.574857e-02,lc};
Point(148) = {4.529467e-01,0.000000e+00,5.495201e-02,lc};
Point(149) = {4.686055e-01,0.000000e+00,5.409108e-02,lc};
Point(150) = {4.842954e-01,0.000000e+00,5.316926e-02,lc};
-Point(151) = {5.000008e-01,0.000000e+00,5.219014e-02,lc};
+Point(151) = {5.000008e-01,0.000000e+00,5.219014e-02,lc3};
Point(152) = {5.157061e-01,0.000000e+00,5.115728e-02,lc};
Point(153) = {5.313960e-01,0.000000e+00,5.007425e-02,lc};
Point(154) = {5.470549e-01,0.000000e+00,4.894463e-02,lc};
@@ -208,10 +210,10 @@ Spline(4) = { 151 ... 200, 1};
Rotate { {1,0,0},{0,0,0},Pi/2 } { Line{1,2,3,4}; }
Translate {-0.5,0,0} { Line{1,2,3,4}; }
-Point(1000) = {2,2,0,lc};
-Point(1001) = {-2,2,0,lc};
-Point(1002) = {-2,-2,0,lc};
-Point(1003) = {2,-2,0,lc};
+Point(1000) = {2,2,0,lc2};
+Point(1001) = {-2,2,0,lc2};
+Point(1002) = {-2,-2,0,lc2};
+Point(1003) = {2,-2,0,lc2};
Line(5) = {1000,1001};
Line(6) = {1001,1002};
Line(7) = {1002,1003};