diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index ec9bb8a5ce590d0733202f9fd16b65f598f68365..f41a834dd2a1abc93934b3ca1646d6763af441f4 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.82 2007-10-14 17:30:42 remacle Exp $
+// $Id: BDS.cpp,v 1.83 2007-10-14 19:54:16 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -778,15 +778,34 @@ void BDS_Mesh::saturate_edge(BDS_Edge * e, std::vector<BDS_Point *> &mids)
bool BDS_SwapEdgeTestParametric::operator () (BDS_Point *_p1,BDS_Point *_p2,
BDS_Point *_q1,BDS_Point *_q2) const
{
- double p1[2] = {_p1->u,_p1->v};
- double p2[2] = {_p2->u,_p2->v};
- double op1[2] = {_q1->u,_q1->v};
- double op2[2] = {_q2->u,_q2->v};
+ double p1[2] = {_p1->u,_p1->v};
+ double p2[2] = {_p2->u,_p2->v};
+ double op1[2] = {_q1->u,_q1->v};
+ double op2[2] = {_q2->u,_q2->v};
- double ori_t1 = gmsh::orient2d(op1, p1, op2);
- double ori_t2 = gmsh::orient2d(op1,op2, p2);
+ double ori_t1 = gmsh::orient2d(op1, p1, op2);
+ double ori_t2 = gmsh::orient2d(op1,op2, p2);
+
+ return( ori_t1 * ori_t2 > 0 ); // the quadrangle was strictly convex !
+
+// double ori_t1 = gmsh::orient2d(p1, p2, op1);
+// double ori_t2 = gmsh::orient2d(p1, p2, op2);
+
+// double ori_t3 = gmsh::orient2d(op1, op2, p1);
+// double ori_t4 = gmsh::orient2d(op1, op2, p2);
+// return (ori_t1 * ori_t2 < 0 && ori_t3 * ori_t4 < 0);
+
+// double t1 = fabs(surface_triangle_param(_p1,_p2,_q1));
+// double t2 = fabs(surface_triangle_param(_p1,_p2,_q2));
+
+// double t3 = fabs(surface_triangle_param(_q1,_q2,_p1));
+// double t4 = fabs(surface_triangle_param(_q1,_q2,_p2));
+
+// // printf("%12.5E %12.5E %12.5E %12.5E -- %12.5E so %1d\n",t1,t2,t3,t4,fabs(t1+t2-t3-t4),fabs(t1+t2-t3-t4) > 1.e-15 * (t1+t2+t3+t4));
+
+// if (fabs(t1+t2-t3-t4) > 1.e-13*(t1+t2+t3+t4))return false;
+// return true;
- return (ori_t1 * ori_t2 > 0);
}
bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
@@ -1174,25 +1193,25 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
double U = 0;
double V = 0;
- // double tot_length = 0;
+ double tot_length = 0;
double LC = 0;
std::list < BDS_Edge * >::iterator eit = p->edges.begin();
while(eit != p->edges.end()) {
if ((*eit)->numfaces() == 1) return false;
BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
- // const double l_e = (*eit)->length();
- U += op->u;
- V += op->v;
- // tot_length += l_e;
+ const double l_e = (*eit)->length();
+ U += op->u*l_e;
+ V += op->v*l_e;
+ tot_length += l_e;
LC += op->lc();
++eit;
}
- U /= (p->edges.size());
- V /= (p->edges.size());
- // U /= tot_length;
- // V /= tot_length;
+ //U /= (p->edges.size());
+ //V /= (p->edges.size());
+ U /= tot_length;
+ V /= tot_length;
LC /= p->edges.size();
std::list < BDS_Face * >ts;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index b9dce76f47b9b9a30072287adef7b9e066cbc3be..296dce3b0e1f4e9d572e528da23db2a2e5fd4852 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.97 2007-10-14 17:30:42 remacle Exp $
+// $Id: meshGFace.cpp,v 1.98 2007-10-14 19:54:16 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -326,10 +326,6 @@ 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);
@@ -393,7 +389,7 @@ void OptimizeMesh(GFace *gf, BDS_Mesh &m, const int NIT)
{
int result = edgeSwapTestQuality(*it,5);
if (result >= 0)
- if(edgeSwapTestDelaunay(*it,gf) || result > 0)
+ if(edgeSwapTestDelaunay(*it,gf))
m.swap_edge ( *it , BDS_SwapEdgeTestParametric());
}
++it;
@@ -544,8 +540,7 @@ void smoothVertexPass ( GFace *gf, BDS_Mesh &m, int &nb_smooth)
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
- {
-
+ {
if(m.smooth_point_centroid(*itp,gf))
nb_smooth ++;
++itp;
@@ -627,15 +622,15 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
splitEdgePass ( gf, m, maxE, nb_split);
//saturateEdgePass ( gf, m, maxE, nb_split);
clock_t t2 = clock();
- swapEdgePass ( gf, m, nb_swap);
+ // swapEdgePass ( gf, m, nb_swap);
clock_t t3 = clock();
collapseEdgePass ( gf, m, minE , MAXNP, nb_collaps);
clock_t t4 = clock();
- swapEdgePass ( gf, m, nb_swap);
+ // swapEdgePass ( gf, m, nb_swap);
clock_t t5 = clock();
smoothVertexPass ( gf, m, nb_smooth);
clock_t t6 = clock();
- swapEdgePass ( gf, m, nb_swap);
+ // swapEdgePass ( gf, m, nb_swap);
clock_t t7 = clock();
// clean up the mesh
@@ -653,6 +648,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
if (nb_split==0 && nb_collaps == 0)break;
}
+
double t_total = t_spl + t_sw + t_col + t_sm;
Msg(DEBUG1," ---------------------------------------");
Msg(DEBUG1," CPU Report ");
@@ -1143,9 +1139,9 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
if (!AlgoDelaunay2D ( gf ))
{
RefineMesh (gf,*m, CTX.mesh.refine_steps);
- OptimizeMesh(gf, *m, 2);
- RefineMesh (gf,*m, -CTX.mesh.refine_steps);
- OptimizeMesh(gf, *m, 2);
+// OptimizeMesh(gf, *m, 2);
+// RefineMesh (gf,*m, -CTX.mesh.refine_steps);
+// OptimizeMesh(gf, *m, 2);
if (gf->meshAttributes.recombine)
{
m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
@@ -1824,7 +1820,7 @@ void meshGFace::operator() (GFace *gf)
// temp fix until we create MEdgeLoops in gmshFace
Msg(DEBUG1, "Generating the mesh");
if(gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
- gmsh2DMeshGenerator(gf,0, false);
+ gmsh2DMeshGenerator(gf,0, true);
}
else{
if(!gmsh2DMeshGeneratorPeriodic(gf,false))