From 806e031649a6dadc827dc064656cc89b01ec5a60 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 30 Jan 2008 15:27:41 +0000
Subject: [PATCH] *** empty log message ***

---
 Common/CommandLine.cpp  |  21 +++++-
 Common/Context.h        |   2 +-
 Common/DefaultOptions.h |   2 +
 Common/Options.cpp      |  13 +++-
 Common/Options.h        |   1 +
 Fltk/Callbacks.cpp      |   3 +-
 Fltk/GUI.cpp            |  10 ++-
 Geo/GEdgeLoop.cpp       |  13 +++-
 Geo/GEdgeLoop.h         |   1 +
 Geo/OCCEdge.cpp         |   3 +-
 Geo/OCCFace.cpp         |   7 +-
 Mesh/BDS.cpp            |  43 +++++++++++-
 Mesh/BackgroundMesh.cpp |   3 +-
 Mesh/meshGEdge.cpp      |   4 +-
 Mesh/meshGFace.cpp      | 148 ++++++++++++++++++++++++++++++----------
 Mesh/meshGFaceBDS.cpp   |  33 +++------
 16 files changed, 233 insertions(+), 74 deletions(-)

diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 749fc1418d..15a3428091 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.113 2008-01-25 21:37:08 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.114 2008-01-30 15:27:40 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -88,6 +88,8 @@ void Print_Usage(char *name){
   Msg(DIRECT, "  -order int            Set mesh order (1, ..., 5)");
   Msg(DIRECT, "  -optimize_hom         Optimize higher order meshes (in 2D)");
   Msg(DIRECT, "  -clscale float        Set characteristic length scaling factor");
+  Msg(DIRECT, "  -clmin   float        Set minimum characteristic length");
+  Msg(DIRECT, "  -clmax   float        Set maximum characteristic length");
   Msg(DIRECT, "  -clcurv               Compute characteristic lengths from curvatures");
   Msg(DIRECT, "  -epslc1d              Set the accuracy of the evaluation of the LCFIELD for 1D mesh");
   Msg(DIRECT, "  -swapangle            Set the treshold angle (in degree) between two adjacent faces below which a swap is allowed");
@@ -362,7 +364,7 @@ void Get_Options(int argc, char *argv[])
         i++;
         if(argv[i] != NULL) {
           CTX.mesh.lc_min = atof(argv[i++]);
-          if(CTX.mesh.lc_factor <= 0.0) {
+          if(CTX.mesh.lc_min <= 0.0) {
             fprintf(stderr, ERROR_STR
                     "Minimum length size must be > 0\n");
             exit(1);
@@ -373,6 +375,21 @@ void Get_Options(int argc, char *argv[])
           exit(1);
         }
       }
+      else if(!strcmp(argv[i] + 1, "clmax")) {
+        i++;
+        if(argv[i] != NULL) {
+          CTX.mesh.lc_max = atof(argv[i++]);
+          if(CTX.mesh.lc_max <= 0.0) {
+            fprintf(stderr, ERROR_STR
+                    "Maximum length size must be > 0\n");
+            exit(1);
+          }
+        }
+        else {
+          fprintf(stderr, ERROR_STR "Missing number\n");
+          exit(1);
+        }
+      }
       else if(!strcmp(argv[i] + 1, "epslc1d")) {
         i++;
         if(argv[i] != NULL) {
diff --git a/Common/Context.h b/Common/Context.h
index dd6e8a60a8..c987ac96f3 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -181,7 +181,7 @@ public :
     int optimize, optimize_netgen, refine_steps;
     int quality_type, label_type;
     double quality_inf, quality_sup, radius_inf, radius_sup;
-    double scaling_factor, lc_factor, rand_factor, lc_integration_precision,lc_min;
+    double scaling_factor, lc_factor, rand_factor, lc_integration_precision,lc_min,lc_max;
     int dual, draw_skin_only;
     int light, light_two_side, light_lines;
     int format, nb_smoothing, algo2d, algo3d, algo_recombine;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 2d9d8e6fd6..279bad8944 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -897,6 +897,8 @@ StringXNumber MeshOptions_Number[] = {
     "Factor applied to all characteristic lengths" },
   { F|O, "CharacteristicLengthMin" , opt_mesh_lc_min, 0.0 ,
     "Minimum mesh size" },
+  { F|O, "CharacteristicLengthMax" , opt_mesh_lc_max, 1.e22 ,
+    "Maximum mesh size" },
   { F|O, "CharacteristicLengthFromCurvature" , opt_mesh_lc_from_curvature , 0. ,
     "Compute characteritic lenghts automatically from curvatures" },
   { F|O, "ColorCarousel" , opt_mesh_color_carousel , 1. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 43b54c8115..d67e17c82f 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.379 2008-01-28 09:59:52 geuzaine Exp $
+// $Id: Options.cpp,v 1.380 2008-01-30 15:27:40 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -4305,6 +4305,17 @@ double opt_mesh_lc_min(OPT_ARGS_NUM)
   return CTX.mesh.lc_min;
 }
 
+double opt_mesh_lc_max(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.mesh.lc_max = val;
+#if defined(HAVE_FLTK)
+  if(WID && (action & GMSH_GUI))
+    WID->mesh_value[26]->value(CTX.mesh.lc_max);
+#endif
+  return CTX.mesh.lc_max;
+}
+
 double opt_mesh_lc_from_curvature(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index a507131e2b..1d33b98811 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -434,6 +434,7 @@ double opt_mesh_tangents(OPT_ARGS_NUM);
 double opt_mesh_explode(OPT_ARGS_NUM);
 double opt_mesh_scaling_factor(OPT_ARGS_NUM);
 double opt_mesh_lc_min(OPT_ARGS_NUM);
+double opt_mesh_lc_max(OPT_ARGS_NUM);
 double opt_mesh_lc_factor(OPT_ARGS_NUM);
 double opt_mesh_lc_from_curvature(OPT_ARGS_NUM);
 double opt_mesh_rand_factor(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index ef35df0d65..685118ebeb 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.560 2008-01-19 22:06:00 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.561 2008-01-30 15:27:40 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1134,6 +1134,7 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
   opt_mesh_nb_smoothing(0, GMSH_SET, WID->mesh_value[0]->value());
   opt_mesh_lc_factor(0, GMSH_SET, WID->mesh_value[2]->value());
   opt_mesh_lc_min(0, GMSH_SET, WID->mesh_value[25]->value());
+  opt_mesh_lc_max(0, GMSH_SET, WID->mesh_value[26]->value());
   opt_mesh_quality_inf(0, GMSH_SET, WID->mesh_value[4]->value());
   opt_mesh_quality_sup(0, GMSH_SET, WID->mesh_value[5]->value());
   opt_mesh_radius_inf(0, GMSH_SET, WID->mesh_value[6]->value());
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 16d6c54c33..4853a2883a 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.650 2008-01-19 22:06:00 geuzaine Exp $
+// $Id: GUI.cpp,v 1.651 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -2436,7 +2436,11 @@ void GUI::create_option_window()
       mesh_value[25]->align(FL_ALIGN_RIGHT);
       mesh_value[25]->callback(mesh_options_ok_cb);
 
-      mesh_value[3] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 7 * BH, IW, BH, "Element order");
+      mesh_value[26] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 7 * BH, IW, BH, "Maximum mesh size");
+      mesh_value[26]->align(FL_ALIGN_RIGHT);
+      mesh_value[26]->callback(mesh_options_ok_cb);
+
+      mesh_value[3] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 8 * BH, IW, BH, "Element order");
       mesh_value[3]->minimum(1);
       // FIXME: this makes it possible to set > 2 by hand, but not by
       // dragging (>2 is too buggy for general use)
@@ -2445,7 +2449,7 @@ void GUI::create_option_window()
       mesh_value[3]->align(FL_ALIGN_RIGHT);
       mesh_value[3]->callback(mesh_options_ok_cb);
 
-      mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 8 * BH, BW, BH, "Use incomplete high order elements (8-node quads, etc.)");
+      mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 9 * BH, BW, BH, "Use incomplete high order elements (8-node quads, etc.)");
       mesh_butt[4]->type(FL_TOGGLE_BUTTON);
       mesh_butt[4]->callback(mesh_options_ok_cb);
 
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 54b997f9eb..39135111ed 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdgeLoop.cpp,v 1.10 2008-01-22 16:57:36 geuzaine Exp $
+// $Id: GEdgeLoop.cpp,v 1.11 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -136,3 +136,14 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
     loop.push_back(ges);
   }
 }
+
+
+GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire, const std::list<int> &dir)
+{
+  std::list<GEdge*>::const_iterator it = cwire.begin();
+  std::list<int>::const_iterator itdir = dir.begin();
+  for ( ; it != cwire.end() ; ++it,++itdir){
+    loop.push_back(GEdgeSigned(*itdir,*it));
+  }
+}
+
diff --git a/Geo/GEdgeLoop.h b/Geo/GEdgeLoop.h
index 2c1a1da0e2..ac2d5f15c7 100644
--- a/Geo/GEdgeLoop.h
+++ b/Geo/GEdgeLoop.h
@@ -45,6 +45,7 @@ public:
   typedef std::list<GEdgeSigned>::iterator iter;
   typedef std::list<GEdgeSigned>::const_iterator citer;
   GEdgeLoop(const std::list<GEdge*> &);
+  GEdgeLoop(const std::list<GEdge*> &, const std::list<int> &dir);
   inline iter begin() { return loop.begin(); }
   inline iter end() { return loop.end(); }
   inline citer begin() const { return loop.begin(); }
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index eac220ac94..087aa55adf 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.31 2008-01-21 19:22:50 geuzaine Exp $
+// $Id: OCCEdge.cpp,v 1.32 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -62,7 +62,6 @@ void OCCEdge::setTrimmed (OCCFace *f)
   }
 }
 
-// Plugin Netgen Salome !
 SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar, int dir) const
 {
   const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 1e859579dd..668c24a34c 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.26 2008-01-20 10:10:41 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.27 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -46,13 +46,15 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
     TopoDS_Shape wire = exp2.Current();
     Msg(DEBUG2,"OCC Face %d - New Wire",num);
     std::list<GEdge*> l_wire;
+    std::list<int> l_oris;
     for(exp3.Init(wire, TopAbs_EDGE); exp3.More(); exp3.Next()){	  
       TopoDS_Edge edge = TopoDS::Edge(exp3.Current());
       int index = emap.FindIndex(edge);
       GEdge *e = m->edgeByTag(index);
       if(!e) throw;
       l_wire.push_back(e);
-      Msg(DEBUG2,"Edge %d",e->tag());
+      l_oris.push_back(edge.Orientation());
+      Msg(DEBUG2,"Edge %d ori %d",e->tag(),edge.Orientation());
       e->addFace(this);
       if(!e->is3D()){
 	OCCEdge *occe = (OCCEdge*)e;
@@ -60,6 +62,7 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
       }
     }      
     
+    //    GEdgeLoop el(l_wire,l_oris);
     GEdgeLoop el(l_wire);
     
     for(GEdgeLoop::citer it = el.begin() ; it != el.end() ; ++it){
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 0864c22859..603ecf9e1c 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.95 2008-01-26 17:47:58 remacle Exp $
+// $Id: BDS.cpp,v 1.96 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -509,6 +509,47 @@ BDS_Mesh ::~ BDS_Mesh ()
 }
 
 
+bool BDS_Mesh::split_face(BDS_Face * f, BDS_Point *mid)
+{
+  BDS_Point *p1 = f->e1->commonvertex (f->e2); 
+  BDS_Point *p2 = f->e3->commonvertex (f->e2); 
+  BDS_Point *p3 = f->e1->commonvertex (f->e3); 
+  BDS_Edge *p1_mid = new BDS_Edge(p1, mid);
+  edges.push_back(p1_mid);
+  BDS_Edge *p2_mid = new BDS_Edge(p2, mid);
+  edges.push_back(p2_mid);
+  BDS_Edge *p3_mid = new BDS_Edge(p3, mid);
+  edges.push_back(p2_mid);
+  BDS_Face *t1, *t2, *t3;
+  t1 = new BDS_Face(f->e1, p1_mid,p3_mid);
+  t2 = new BDS_Face(f->e2, p2_mid,p1_mid);
+  t3 = new BDS_Face(f->e3, p3_mid,p2_mid);
+
+  t1->g = f->g;
+  t2->g = f->g;
+  t3->g = f->g;
+
+  p1_mid->g = f->g;
+  p2_mid->g = f->g;
+  p3_mid->g = f->g;
+
+  mid->g = f->g;
+
+  triangles.push_back(t1);
+  triangles.push_back(t2);
+  triangles.push_back(t3);
+
+  // config has changed
+  p1->config_modified = true;
+  p2->config_modified = true;
+  p3->config_modified = true;
+
+  // delete the face
+  del_face(f);
+
+  return true;
+}
+
 bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
 {
   /*
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 11324b4440..c9f60b1312 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.32 2008-01-19 22:06:03 geuzaine Exp $
+// $Id: BackgroundMesh.cpp,v 1.33 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -197,6 +197,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
     lc = std::min (lc,LC_MVertex_CURV(ge, U, V));
 
   lc = std::max(lc,CTX.mesh.lc_min*CTX.mesh.lc_factor);
+  lc = std::min(lc,CTX.mesh.lc_max*CTX.mesh.lc_factor);
 
   if(lc <= 0.){
     Msg(GERROR, "Incorrect char. length lc = %g: using default instead", lc);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 899a05f4db..5d618b1be9 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.51 2008-01-21 19:22:50 geuzaine Exp $
+// $Id: meshGEdge.cpp,v 1.52 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -303,7 +303,7 @@ void meshGEdge::operator() (GEdge *ge)
   // Send a messsage to the GMSH environment
 
   if (length == 0.0)
-    Msg(GERROR,"Curve %d has a zero length", ge->tag());
+    Msg(DEBUG2,"Curve %d has a zero length", ge->tag());
 
   List_Reset(Points);
     
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 02dd505bd1..13e8c9fa4a 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.113 2008-01-26 17:47:58 remacle Exp $
+// $Id: meshGFace.cpp,v 1.114 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -824,15 +824,63 @@ inline double dist2 (const SPoint2 &p1,const SPoint2 &p2)
   return dx*dx+dy*dy;
 }
 
-
+// bool buildConsecutiveListOfVertices_b (  GFace *gf,
+// 					 GEdgeLoop  &gel , 
+// 					 std::vector<BDS_Point*> &result,
+// 					 SBoundingBox3d &bbox,
+// 					 BDS_Mesh *m,
+// 					 std::map<BDS_Point*,MVertex*> &recover_map, 
+// 					 int &count, double tol){
+//   std::map<GEntity*,std::vector<SPoint2> > meshes;
+//   std::map<GEntity*,std::vector<SPoint2> > meshes_seam;  
+
+//   result.clear();
+  
+//   GEdgeLoop::iter it  = gel.begin();  
+
+//   while (it != gel.end())   
+//    {
+//      // I get the signed edge
+//      GEdgeSigned ges = *it ;      
+//      std::vector<SPoint2> mesh1d;
+//      // I look if it is a seam
+//      bool seam = ges.ge->isSeam(gf);
+//      // I get parameter bounds
+//      Range<double> range = ges.ge->parBounds(0);
+//      // I Get the first vertex
+//      MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0];
+//      if ( seam && ges._sign == 0 )mesh1d.push_back(ges.ge->reparamOnFace(gf,range.low(),-1));
+//      else mesh1d.push_back(ges.ge->reparamOnFace(gf,range.low(),1));
+//      for (unsigned int i=0;i<ges.ge->mesh_vertices.size();i++)
+//        {
+// 	 double u;
+// 	 here = ges.ge->mesh_vertices[i];
+// 	 here->getParameter(0,u);
+// 	 if ( seam && ges._sign == 0) mesh1d.push_back(ges.ge->reparamOnFace(gf,u,-1));
+// 	 else mesh1d.push_back(ges.ge->reparamOnFace(gf,u,1));
+//        }
+//      here = ges.ge->getEndVertex()->mesh_vertices[0];
+//      if ( seam && ges._sign == 0) mesh1d.push_back(ges.ge->reparamOnFace(gf,range.high(),-1));
+//      else mesh1d.push_back(ges.ge->reparamOnFace(gf,range.high(),1));
+//      it++;
+//    }
+// }
+
+
+static void printMesh1d (int iEdge, int seam, std::vector<SPoint2> &m){
+  printf("Mesh1D for edge %d seam %d\n",iEdge,seam);
+  for (int i=0;i<m.size();i++){
+    printf("%12.5E %12.5E\n",m[i].x(),m[i].y());
+  }
+}
 
 bool buildConsecutiveListOfVertices (  GFace *gf,
 				       GEdgeLoop  &gel , 
 				       std::vector<BDS_Point*> &result,
 				       SBoundingBox3d &bbox,
 				       BDS_Mesh *m,
-				       std::map<BDS_Point*,MVertex*> &recover_map, 
-				       int &count, double tol)
+				       std::map<BDS_Point*,MVertex*> &recover_map_global, 
+				       int &count, int countTot, double tol, bool seam_the_first = false)
 {
 
   // for each edge, we build a list of points that 
@@ -840,18 +888,19 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
   // for seams, we build the list for every side
   // for closed loops, we build it on both senses
 
-  //  printf("new edge loop\n");
   std::map<GEntity*,std::vector<SPoint2> > meshes;
   std::map<GEntity*,std::vector<SPoint2> > meshes_seam;  
 
   const int _DEBUG = false;
 
+  std::map<BDS_Point*,MVertex*> recover_map; 
+
   result.clear();
+  count = 0;
   
   GEdgeLoop::iter it  = gel.begin();  
 
-
-  if (_DEBUG)printf("face %d with %d edges\n",gf->tag(), (int)gf->edges().size());
+  if (_DEBUG)printf("face %d with %d edges case %d\n",gf->tag(), (int)gf->edges().size(),seam_the_first);
 
   while (it != gel.end())   
    {
@@ -862,8 +911,6 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 
      bool seam = ges.ge->isSeam(gf);
      
-     if (_DEBUG)printf("face %d edge %d seam %d (%d %d)\n",gf->tag(),ges.ge->tag(),seam,ges.ge->getBeginVertex()->tag(),ges.ge->getEndVertex()->tag());
-     
      Range<double> range = ges.ge->parBounds(0);
 
      MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0];
@@ -882,6 +929,10 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
      if ( seam ) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf,range.high(),-1));
      meshes.insert(std::pair<GEntity*,std::vector<SPoint2> > (ges.ge,mesh1d) );
      if(seam)meshes_seam.insert(std::pair<GEntity*,std::vector<SPoint2> > (ges.ge,mesh1d_seam) );     
+
+     //     printMesh1d (ges.ge->tag(), seam, mesh1d);
+     //     if (seam)printMesh1d (ges.ge->tag(), seam, mesh1d_seam);
+
      it++;
    }
 
@@ -896,6 +947,7 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 
   while (unordered.size())
    {
+     if (_DEBUG)printf("unordered.size() = %d\n",unordered.size());
      std::list<GEdgeSigned>::iterator it = unordered.begin();     
      std::vector<SPoint2>  coords;
 
@@ -914,18 +966,29 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 	 if (!counter)
 	   {
 	     counter++;
-	     coords = ((*it)._sign == 1)?mesh1d:mesh1d_reversed;	 
-	     found = (*it);
+	     if (seam && seam_the_first){
+	       coords = ((*it)._sign == 1)?mesh1d_seam:mesh1d_seam_reversed;	 
+	       found = (*it);
+	       Msg(INFO,"This test case would have failed in Previous Gmsh Version ;-)");
+	     }
+	     else{
+	       coords = ((*it)._sign == 1)?mesh1d:mesh1d_reversed;	 
+	       found = (*it);
+	     }
 	     unordered.erase(it);
+	     if (_DEBUG)printf("Starting with edge = %d seam %d\n",(*it).ge->tag(),seam);
 	     break;
 	   }
 	 else
 	   {
+	     if (_DEBUG)printf("Followed by edge = %d\n",(*it).ge->tag());
 	     SPoint2 first_coord         = mesh1d[0];
-	     double d = dist2(last_coord,first_coord);
+	     double d=-1,d_reversed=-1,d_seam=-1,d_seam_reversed=-1;
+	     d = dist2(last_coord,first_coord);
+	     if (_DEBUG)printf("%g %g dist = %12.5E\n",first_coord.x(),first_coord.y(),d);
 	     if (d < tol) 
 	       {
-		 if (_DEBUG)printf("d = %12.5E %d\n",d, (int)coords.size());
+		 //		 if (_DEBUG)printf("d = %12.5E %d\n",d, (int)coords.size());
 		 coords.clear();
 		 coords = mesh1d;
 		 found = GEdgeSigned(1,ge);
@@ -933,11 +996,11 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 		 goto Finalize;
 	       }
 	     SPoint2 first_coord_reversed = mesh1d_reversed[0];
-	     double d_reversed = dist2(last_coord,first_coord_reversed);
-	     if (_DEBUG)printf("d_r = %12.5E\n",d_reversed);
+	     d_reversed = dist2(last_coord,first_coord_reversed);
+	     if (_DEBUG)printf("%g %g dist_reversed = %12.5E\n",first_coord_reversed.x(),first_coord_reversed.y(),d_reversed);
 	     if (d_reversed < tol) 
 	       {
-		 if (_DEBUG)printf("d_r = %12.5E\n",d_reversed);
+		 //		 if (_DEBUG)printf("d_r = %12.5E\n",d_reversed);
 		 coords.clear();
 		 coords = mesh1d_reversed;
 		 found = (GEdgeSigned(-1,ge));
@@ -948,20 +1011,20 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 	       {
 		 SPoint2 first_coord_seam         = mesh1d_seam[0];
 		 SPoint2 first_coord_seam_reversed = mesh1d_seam_reversed[0];
-		 double d_seam = dist2(last_coord,first_coord_seam);
+		 d_seam = dist2(last_coord,first_coord_seam);
+		 if (_DEBUG)printf("dist_seam = %12.5E\n",d_seam);
 		 if (d_seam < tol)
 		   {
-		     if (_DEBUG)printf("d_seam = %12.5E\n",d_seam);
 		     coords.clear();
 		     coords = mesh1d_seam;
 		     found = (GEdgeSigned(1,ge));
 		     unordered.erase(it);
 		     goto Finalize;
 		   }
-		 double d_seam_reversed = dist2(last_coord,first_coord_seam_reversed);
+		 d_seam_reversed = dist2(last_coord,first_coord_seam_reversed);
+		 if (_DEBUG)printf("dist_seam_reversed = %12.5E\n",d_seam_reversed);
 		 if (d_seam_reversed < tol)
 		   {
-		     if (_DEBUG)printf("d_seam_reversed = %12.5E\n",d_seam_reversed);
 		     coords.clear();
 		     coords = mesh1d_seam_reversed;
 		     found = (GEdgeSigned(-1,ge));
@@ -975,8 +1038,17 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
        }
    Finalize:
 
-     if (_DEBUG)printf("Finalize\n");
-     if (coords.size() == 0)return false;
+     if (_DEBUG)printf("Finalize, found %d points\n",coords.size());
+     if (coords.size() == 0){
+       // It has not worked : either tolerance is wrong or the first seam edge
+       // has to be taken with the other parametric coordinates (because it is
+       // only present once in the closure of the domain). 
+       for (std::map<BDS_Point*,MVertex*> :: iterator it = recover_map.begin();
+	    it != recover_map.end(); ++it){
+	 m->del_point(it->first);
+       }
+       return false;
+     }
      
      std::vector<MVertex*>    edgeLoop;
      if ( found._sign == 1)
@@ -1003,7 +1075,7 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 	 SPoint2 param = coords [i];
 	 U = param.x() / m->scalingU ;
 	 V = param.y() / m->scalingV;	
-	 BDS_Point *pp = m->add_point ( count, U,V,gf );
+	 BDS_Point *pp = m->add_point ( count+countTot, U,V,gf );
  	 if(ge->dim() == 0)
  	   {
 	     pp->lcBGM() = BGM_MeshSize(ge,0,0,here->x(),here->y(),here->z());
@@ -1027,7 +1099,7 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 	 m->add_geom (ge->tag(), ge->dim());
 	 BDS_GeomEntity *g = m->get_geom(ge->tag(),ge->dim());
 	 pp->g = g;
-	 if (_DEBUG)printf("point %3d (%8.5f %8.5f) (%2d,%2d)\n",count,pp->u,pp->v,pp->g->classif_tag,pp->g->classif_degree);
+	 if (_DEBUG)printf("point %3d (%8.5f %8.5f : %8.5f %8.5f) (%2d,%2d)\n",count,pp->u,pp->v,param.x(),param.y(),pp->g->classif_tag,pp->g->classif_degree);
 	 bbox += SPoint3(U,V,0);	  
 	 edgeLoop_BDS.push_back(pp);
 	 recover_map[pp] = here;	 
@@ -1048,6 +1120,8 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
 //       {
 //         printf("point %3d (%8.5f %8.5f) (%2d,%2d)\n",i,result[i]->u,result[i]->v,result[i]->g->classif_tag,result[i]->g->classif_degree);
 //       }
+  // It has worked, so we add all the points to the recover map
+  recover_map_global.insert(recover_map.begin(),recover_map.end());
 
   return true;
 }
@@ -1082,15 +1156,19 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
     for (std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin() ; it != gf->edgeLoops.end() ; it++)
       {
 	std::vector<BDS_Point* > edgeLoop_BDS;
-	if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-7*LC2D))
-	  if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-5*LC2D))
-	    if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal, 1.e-3*LC2D))
-	      {
-		gf->meshStatistics.status = GFace::FAILED;
-		Msg(GERROR,"The 1D Mesh seems not to be forming a closed loop");
-		m->scalingU = m->scalingV = 1.0;
-		return false;
-	      }
+	int nbPointsLocal;
+	if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsLocal, nbPointsTotal, 1.e-7*LC2D))
+	  if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsLocal, nbPointsTotal, 1.e-7*LC2D,true))
+	    if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsLocal, nbPointsTotal, 1.e-5*LC2D))
+	      if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsLocal, nbPointsTotal, 1.e-5*LC2D,true))
+		if(!buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsLocal, nbPointsTotal, 1.e-3*LC2D))
+		  {
+		    gf->meshStatistics.status = GFace::FAILED;
+		    Msg(GERROR,"The 1D Mesh seems not to be forming a closed loop");
+		    m->scalingU = m->scalingV = 1.0;
+		    return false;
+		  }
+	nbPointsTotal += nbPointsLocal;
 	edgeLoops_BDS.push_back(edgeLoop_BDS);
       }
   }
@@ -1114,7 +1192,7 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
 	  doc.points[count].where.h = U + XX;
 	  doc.points[count].where.v = V + YY;
 	  doc.points[count].adjacent = NULL;
-	  doc.points[count].data = pp;      
+	  doc.points[count].data = pp; 
 	  count++;	  
 	}
     }
@@ -1420,7 +1498,7 @@ void meshGFace::operator() (GFace *gf)
   computeEdgeLoops(gf, points, indices);
 
   // temp fix until we create MEdgeLoops in gmshFace
-  if (1 || gf->tag() == 6)
+  if (1 || gf->tag() == 46)
     {
       Msg(DEBUG1, "Generating the mesh");
       if(noseam (gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index dbb4cfe6db..13de94d874 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceBDS.cpp,v 1.2 2008-01-26 17:47:58 remacle Exp $
+// $Id: meshGFaceBDS.cpp,v 1.3 2008-01-30 15:27:41 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -713,25 +713,14 @@ void gmshOptimizeMeshBDS(GFace *gf,
 
 
 
-// void delaunayPointInsertionBDS ( GFace *gf, BDS_Mesh &m, BDS_Vertex *v, BDS_Face *f){
-//   const double p[2] = {v->u,v->v};
-
-//   BDS_Edge *e1 = f->e1;
-//   BDS_Edge *e2 = f->e2;
-//   BDS_Edge *e3 = f->e3;
-
-//   m.split_face ( f , v );
+void delaunayPointInsertionBDS ( GFace *gf, BDS_Mesh &m, BDS_Point *v, BDS_Face *f){
+  const double p[2] = {v->u,v->v};
   
-//   recur_delaunay_swap ( e1 );
-//   recur_delaunay_swap ( e2 );
-//   recur_delaunay_swap ( e3 );
-
-//   return;
-
-//   BDS_Point *n2[4];
-//   f1->getNodes(n1);
-
-//   const double a[2] = {v->u,v->v};
-//   const double b[2] = {v->u,v->v};
-//   const double c[2] = {v->u,v->v};
-// } 
+  BDS_Edge *e1 = f->e1;
+  BDS_Edge *e2 = f->e2;
+  BDS_Edge *e3 = f->e3;
+  // face is splitted, 
+  m.split_face ( f , v );
+  int nb_swap = 0;
+  gmshDelaunayizeBDS ( gf, m, nb_swap );
+}
-- 
GitLab