diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 8ab4a9d14013793caecda8acf16e8e58941073f2..b439812fa74153a3a8a14841412d978bcc2a2aed 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -49,6 +49,9 @@ class GEdge : public GEntity {
   virtual bool continuous(int dim=0) const = 0;
   virtual void setVisibility(char val, bool recursive=false);
 
+  /** True if the edge is a seam for the given face. */
+  virtual int isSeam(GFace *face) const {return 0;}
+
   // The bounding box
   SBoundingBox3d bounds() const;
 
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 7c930cdcc7e4c4818e5c828a9246b01eb137c368..8589767f4a2a54b4f7d0a19654962c6df572a9a3 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -105,6 +105,7 @@ GEdgeLoop::GEdgeLoop ( const std::list<GEdge*> & cwire )
   while (wire.size())
     {
       ges = nextOne ( prevOne , wire );
+      
       prevOne = &ges;
       ges.print();
       loop.push_back(ges);
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index d7e4d8625d0cde55ee8b8e6b0c75878861198ed4..9ae7a5e18e14544ee61f88f1034245d0f9a88c17 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.20 2006-11-16 18:32:41 remacle Exp $
+// $Id: GFace.cpp,v 1.21 2006-11-20 12:44:09 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -475,10 +475,3 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p) const
   return SPoint2(U,V);
 }
 
-void GFace::parFromPoint(const SPoint3 &p, std::list<double> &u, std::list<double> &v ) const
-{
-  double U,V;  
-  XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0);
-  u.push_back(U);
-  v.push_back(V);
-}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 3916b146f4b18c7fa327356afdc1976e3813a290..b25055dd9e08b2afb70739b70b49b89154422ba4 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -44,7 +44,6 @@ class GFace : public GEntity
 {
  protected: 
   // edge loops, will replace what follows
-  std::list<GEdgeLoop> edgeLoops;
   // list of al the edges of the face
   std::list<GEdge *> l_edges;
   std::list<int> l_dirs;
@@ -96,7 +95,6 @@ class GFace : public GEntity
   // Return the parmater location on the face given a point in space
   // that is on the face.
   virtual SPoint2            parFromPoint(const SPoint3 &) const;
-  virtual void               parFromPoint(const SPoint3 &, std::list<double> &u, std::list<double> &v ) const;
 
   // True if the parameter value is interior to the face.
   virtual int containsParam(const SPoint2 &pt) const = 0;
@@ -143,6 +141,7 @@ class GFace : public GEntity
 
   std::vector<MTriangle*> triangles;
   std::vector<MQuadrangle*> quadrangles;
+  std::list<GEdgeLoop> edgeLoops;
 
   struct {
     // do we recombine the triangles of the mesh ?
@@ -158,6 +157,7 @@ class GFace : public GEntity
     int transfiniteArrangement;
     // the extrusion parameters (if any)
     ExtrudeParams *extrude;
+    // edge loops
   } meshAttributes ;
 
 };
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index cf847ffb918cf9d7d9455e990b520a752eccbe3d..b792a7a65f1796d5e385478df6949973c2617628 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: GVertex.cpp,v 1.7 2006-11-16 21:14:10 remacle Exp $
+// $Id: GVertex.cpp,v 1.8 2006-11-20 12:44:09 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -45,7 +45,7 @@ void GVertex::delEdge(GEdge *e)
 }
 
 
-SPoint2 GVertex::reparamOnFace ( GFace *gf ) const
+SPoint2 GVertex::reparamOnFace ( GFace *gf , int) const
 {
   return gf->parFromPoint ( SPoint3(x(),y(),z() ));
 }
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index 91b7a5c069480da08486bfb86b8660f406b2e310..ffa6687378eaa59266e9184a0d4ba98fb091272f 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -44,7 +44,7 @@ class GVertex : public GEntity
   virtual GeomType geomType() const {return Point;}
   virtual double prescribedMeshSizeAtVertex() const {return 0;}
   virtual SBoundingBox3d bounds(){ return SBoundingBox3d(SPoint3(x(), y(), z())); }
-  virtual SPoint2 reparamOnFace ( GFace *gf ) const;
+  virtual SPoint2 reparamOnFace ( GFace *gf , int) const;
   virtual std::string getAdditionalInfoString();
 };
 
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index f9b6ce720fbb1dd2517e3b0158d9dcaea1a92432..72880a4232fb42d4a93a745d1092ad464df00597 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -100,7 +100,7 @@ class MVertex{
 };
 
 class MEdgeVertex : public MVertex{
- private:
+ protected:
   double _u;
  public :
   MEdgeVertex(double x, double y, double z, GEntity *ge, double u) 
@@ -112,6 +112,19 @@ class MEdgeVertex : public MVertex{
   virtual bool setParameter(int i, double par){ _u = par; return true; }
 };
 
+class MSeamEdgeVertex : public MEdgeVertex{
+ private:
+  MEdgeVertex *_ev;
+ public :
+  MSeamEdgeVertex(MEdgeVertex* ev) 
+    : MEdgeVertex(ev->x(), ev->y(), ev->z(), ev->onWhat(),0), _ev(ev)
+    {
+      _ev->getParameter(0, _u);
+    }
+  ~MSeamEdgeVertex(){}
+};
+
+
 class MFaceVertex : public MVertex{
  protected:
   double _u, _v;
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 22724d5af3f0bb9c140690de2e10744d418deffb..1fe1b794da907bb2eda4e0c983876ccabff0ddb7 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.7 2006-11-16 21:14:10 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.8 2006-11-20 12:44:09 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -52,17 +52,57 @@ void OCCEdge::setTrimmed (OCCFace *f)
 
 SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
 {
-  double t0,t1;
   const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
+  double t0,t1;
   Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface(c, *s, t0, t1);
   if (c2d.IsNull())
     return GEdge::reparamOnFace(face, epar,dir);
-
   double u,v;
   c2d->Value(epar).Coord(u,v);
-  return SPoint2 (u,v);
+  if (! isSeam ( face ) )
+    {
+      return SPoint2 (u,v);
+    }
+  else
+    {
+      BRepAdaptor_Surface surface( *s );
+
+//       printf ("surface %d (%d %d) firstu %g lastu %g firstv %g lastv %g\n",
+// 	      surface.IsUPeriodic() ,
+// 	      surface.IsVPeriodic() ,
+// 	      face->tag(),
+// 	      surface.FirstUParameter(),
+// 	      surface.LastUParameter(),
+// 	      surface.FirstVParameter(),
+// 	      surface.LastVParameter());
+
+      if ( surface.IsUPeriodic() )
+	{
+	  if (dir == -1) 
+	    return SPoint2(surface.FirstUParameter(), v);
+	  else
+	    return SPoint2(surface.LastUParameter(), v);
+	}
+      else {
+	if (dir == -1) 
+	  return SPoint2(u , surface.FirstVParameter());
+	else
+	  return SPoint2(u , surface.LastVParameter());
+      }
+    }
 }
 
+/** True if the edge is a seam for the given face. */
+int OCCEdge::isSeam(GFace *face) const
+{
+  const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
+  BRepAdaptor_Surface surface( *s );
+  if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
+    {
+      return BRep_Tool::IsClosed( c, *s );
+    }
+  return 0;
+}
 
 
 GPoint OCCEdge::point(double par) const
diff --git a/Geo/OCCEdge.h b/Geo/OCCEdge.h
index 6a045ad6c14836c156ef820666e62d9f00e66c99..f3a909e69c61c524cf84eab2eae49167cdabe2ae 100644
--- a/Geo/OCCEdge.h
+++ b/Geo/OCCEdge.h
@@ -60,6 +60,7 @@ class OCCEdge : public GEdge {
   virtual int minimumDrawSegments () const;
   bool is3D() const {return !curve.IsNull();}
   void setTrimmed (OCCFace *);
+  int isSeam (GFace *) const;
 };
 
 #endif
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index da730f051b3e90652672c428f71790586c1c07a0..19a2be842d127fee5dd73dc2c54270ae7e40d9b4 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.8 2006-11-16 21:14:10 remacle Exp $
+// $Id: OCCFace.cpp,v 1.9 2006-11-20 12:44:09 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -25,8 +25,12 @@
 #include "OCCEdge.h"
 #include "OCCFace.h"
 #include "Message.h"
+#include "Utils.h"
 
 #if defined(HAVE_OCC)
+#include "Geom_CylindricalSurface.hxx"
+#include "Geom_Plane.hxx"
+#include "gp_Pln.hxx"
 
 OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape &emap)
   : GFace(m, num), s(_s),_periodic(false)
@@ -156,39 +160,17 @@ SPoint2 OCCFace::parFromPoint(const SPoint3 &qp) const
       Msg(GERROR,"OCC Project Point on Surface FAIL");
       return GFace::parFromPoint(qp);
     }   
-//   if (proj.NbPoints() != 1)
-//     {
-//       Msg(WARNING,"More than one points match the projection (%d) ", proj.NbPoints());
-//     }
-
   double U,V;
   proj.LowerDistanceParameters (U, V);
-  //  proj.Parameters (1,U,V);
   return SPoint2(U,V);
 }
 
-void OCCFace::parFromPoint(const SPoint3 &qp, std::list<double> &u, std::list<double> &v ) const
-{
-  gp_Pnt pnt(qp.x(),qp.y(),qp.z());
-  GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
-  if (!proj.NbPoints())
-    {
-      Msg(GERROR,"OCC Project Point on Surface FAIL");
-      GFace::parFromPoint(qp,u,v);
-      return;
-    }   
-
-  for (int i = 1;i <= proj.NbPoints();i++)
-    {
-      double U,V;
-      proj.Parameters (i,U,V);
-      u.push_back(U);
-      v.push_back(V);
-    }
-}
-
 GEntity::GeomType OCCFace::geomType() const
 {
+  if (occface->DynamicType() == STANDARD_TYPE(Geom_Plane))
+    return Plane;
+  else if (occface->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface))
+    return Cylinder;
   return Unknown;
 }
 
@@ -208,7 +190,45 @@ double OCCFace::curvature (const SPoint2 &param) const
 
 int OCCFace::containsPoint(const SPoint3 &pt) const
 { 
-  Msg(GERROR,"Not Done Yet ...");
+  if(geomType() == Plane)
+    {
+      gp_Pln pl = Handle(Geom_Plane)::DownCast(occface)->Pln();
+
+      // OK to use the normal from the mean plane here: we compensate
+      // for the (possibly wrong) orientation at the end
+      double n[3] , c;
+      pl.Coefficients ( n[0], n[1], n[2], c );
+      norme(n);
+      double angle = 0.;
+      Vertex v,v1,v2;
+      v.Pos.X = pt.x();
+      v.Pos.Y = pt.y();
+      v.Pos.Z = pt.z();
+      for(std::list<GEdge*>::const_iterator  it = l_edges.begin(); it != l_edges.end(); it++) {
+	GEdge *c = *it;
+	int N=10;
+	Range<double> range = c->parBounds(0);
+	for(int j = 0; j < N ; j++) {
+	  double u1 = (double)j / (double)N;
+	  double u2 = (double)(j + 1) / (double)N;
+	  GPoint pp1 = c->point(range.low() + u1 * (range.high() - range.low() ));
+	  GPoint pp2 = c->point(range.low() + u2 * (range.high() - range.low() ));
+	  v1.Pos.X = pp1.x();
+	  v1.Pos.Y = pp1.y();
+	  v1.Pos.Z = pp1.z();
+	  v2.Pos.X = pp2.x();
+	  v2.Pos.Y = pp2.y();
+	  v2.Pos.Z = pp2.z();	  
+	  angle += angle_plan(&v, &v1, &v2, n);
+	}
+      }
+      // we're inside if angle equals 2 * pi
+      if(fabs(angle) > 2 * M_PI - 0.5 && fabs(angle) < 2 * M_PI + 0.5) 
+	return true;
+      return false;
+    }
+  else
+    Msg(GERROR,"Not Done Yet ...");
 }
 // void OCCFace::buildVisTriangulation ();
 // {
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index 803843f865c523bf8b692188f2f87a008db1c6f4..25a84fbe4165789d2bb525fedfff0230a5296d6f 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -63,7 +63,6 @@ class OCCFace : public GFace {
   void * getNativePtr() const { return (void*)&s; }
   virtual bool surfPeriodic(int dim) const {return _periodic;}
   virtual SPoint2 parFromPoint(const SPoint3 &) const;
-  virtual void parFromPoint(const SPoint3 &, std::list<double> &u, std::list<double> &v ) const;
   virtual double curvature(const SPoint2 &param) const;
 };
 
diff --git a/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index 8ac63a82ded3115fd37405e574df3334c845a741..0268d7210c3f822d1245d6465628b2cc57d6773a 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCVertex.cpp,v 1.3 2006-11-16 21:14:10 remacle Exp $
+// $Id: OCCVertex.cpp,v 1.4 2006-11-20 12:44:09 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -33,7 +33,7 @@ double max_surf_curvature ( const GVertex *gv, double x, double y, double z , co
   double curv = 0;
   while (it != faces.end())
     {
-      SPoint2 par = gv->reparamOnFace((*it));
+      SPoint2 par = gv->reparamOnFace((*it),1);
       double cc = (*it)->curvature ( par );
       if (cc < 1.e2) curv = std::max(curv, cc );					
       ++it;
@@ -55,13 +55,32 @@ double OCCVertex::max_curvature_of_surfaces() const
   return max_curvature;
 }
 
-SPoint2 OCCVertex::reparamOnFace ( GFace *gf ) const
+SPoint2 OCCVertex::reparamOnFace ( GFace *gf , int dir) const
 {
   std::list<GEdge*>::const_iterator it = l_edges.begin();
   while (it != l_edges.end())
     {
       std::list<GEdge*> l_edges = gf->edges();
-      
+      if (std::find(l_edges.begin(),l_edges.end(),*it) != l_edges.end())
+	{
+	  if ((*it)->isSeam(gf))
+	    {
+	      const TopoDS_Face *s = (TopoDS_Face*) gf->getNativePtr();
+	      const TopoDS_Edge *c = (TopoDS_Edge*) (*it)->getNativePtr();
+	      double s1,s0;
+	      Handle(Geom2d_Curve) curve2d = BRep_Tool::CurveOnSurface(*c, *s, s0, s1);
+	      if ((*it)->getBeginVertex() == this)
+		return (*it)->reparamOnFace(gf,s0,dir);
+	      else if ((*it)->getEndVertex() == this)
+		return (*it)->reparamOnFace(gf,s1,dir);
+	    }
+	}
+      ++it;
+    }  
+  it = l_edges.begin();
+  while (it != l_edges.end())
+    {
+      std::list<GEdge*> l_edges = gf->edges();
       if (std::find(l_edges.begin(),l_edges.end(),*it) != l_edges.end())
 	{
 	  const TopoDS_Face *s = (TopoDS_Face*) gf->getNativePtr();
@@ -69,9 +88,9 @@ SPoint2 OCCVertex::reparamOnFace ( GFace *gf ) const
 	  double s1,s0;
 	  Handle(Geom2d_Curve) curve2d = BRep_Tool::CurveOnSurface(*c, *s, s0, s1);
 	  if ((*it)->getBeginVertex() == this)
-	    return (*it)->reparamOnFace(gf,s0,1);
+	    return (*it)->reparamOnFace(gf,s0,dir);
 	  else if ((*it)->getEndVertex() == this)
-	    return (*it)->reparamOnFace(gf,s1,1);
+	    return (*it)->reparamOnFace(gf,s1,dir);
 	}
       ++it;
     }
diff --git a/Geo/OCCVertex.h b/Geo/OCCVertex.h
index 6fbd53d906850cfdfcbacc52f4b6c83fb131977b..a34f2b2eaca0ef6f65998115b811c269a815bcfc 100644
--- a/Geo/OCCVertex.h
+++ b/Geo/OCCVertex.h
@@ -69,7 +69,7 @@ class OCCVertex : public GVertex {
       lc = std::min (lc,6.28/(CTX.mesh.min_circ_points*maxc));
     return lc;
   }
-  virtual SPoint2 reparamOnFace ( GFace *gf ) const;
+  virtual SPoint2 reparamOnFace ( GFace *gf , int) const;
 };
 
 #endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 0e04a234b8beac3b0b363b65fe7ab00fb816380b..8da54f41a6d6ba32ded8f890158a4b1efd7f65e7 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.19 2006-11-16 21:14:10 remacle Exp $
+// $Id: meshGFace.cpp,v 1.20 2006-11-20 12:44:09 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -107,7 +107,7 @@ public :
     if (ge->dim() == 0)
       {
 	GVertex *gve = (GVertex*)ge;
-	SPoint2 param = gve->reparamOnFace(gf);
+	SPoint2 param = gve->reparamOnFace(gf,1);
 	v->x() = param.x();  
 	v->y() = param.y();
 	v->z() = 0.0;
@@ -400,7 +400,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 	    {
 	      (*it)->p1->config_modified = false;
 	      (*it)->p2->config_modified = false;
-	      double lone = NewGetLc ( *it);
+	      double lone = NewGetLc ( *it);	      
 	      maxL = std::max(maxL,lone);
 	      minL = std::min(minL,lone);
 	    }
@@ -722,29 +722,27 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
 // that respects the boundaries of the
 // domain, including embedded points 
 // and surfaces
+
 void gmsh2DMeshGenerator ( GFace *gf )
 {
 
-  std::set<MVertex*>all_vertices;
+  typedef std::set<MVertex*> v_container ;
+  v_container all_vertices;
   std::map<int, MVertex*>numbered_vertices;
   std::list<GEdge*> edges = gf->edges();
   std::list<GEdge*> emb_edges = gf->emb_edges();
   std::list<GEdge*>::iterator it = edges.begin();
-  
+
   // build a set with all points of the boundaries
+  it = edges.begin();
   while(it != edges.end())
     {
       all_vertices.insert ( (*it)->mesh_vertices.begin() , (*it)->mesh_vertices.end() );
       all_vertices.insert ( (*it)->getBeginVertex()->mesh_vertices.begin() , (*it)->getBeginVertex()->mesh_vertices.end() );
-//       GVertex *v1 = (*it)->getBeginVertex();
-//       GVertex *v2 = (*it)->getEndVertex();
-//       MVertex *mv1 = *(v1->mesh_vertices.begin());
-//       MVertex *mv2 = *(v2->mesh_vertices.begin());
-//       Msg(INFO,"%d -> %d = %g %g %g -- %g %g %g",v1->tag(),v2->tag(),v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z());
-//       Msg(INFO,"%g %g %g -- %g %g %g",mv1->x(),mv1->y(),mv1->z(),mv2->x(),mv2->y(),mv2->z());
       all_vertices.insert ( (*it)->getEndVertex()->mesh_vertices.begin() , (*it)->getEndVertex()->mesh_vertices.end() );
       ++it;
     }
+
   it = emb_edges.begin();
   while(it != emb_edges.end())
     {
@@ -754,12 +752,11 @@ void gmsh2DMeshGenerator ( GFace *gf )
       ++it;
     }
 
-
   double * X_ = new double [all_vertices.size()];
   double * Y_ = new double [all_vertices.size()];
   double * Z_ = new double [all_vertices.size()];
 
-  std::set<MVertex*>::iterator itv = all_vertices.begin();
+  v_container::iterator itv = all_vertices.begin();
 
   //  FILE *fdeb = fopen("debug.dat","w");
   //  fprintf(fdeb,"surface %d\n" ,gf->tag());
@@ -1047,6 +1044,398 @@ void gmsh2DMeshGenerator ( GFace *gf )
   //  std::for_each(gf->mesh_vertices.begin(),gf->mesh_vertices.end(),p2c);    
 }
 
+void buildConsecutiveListOfVertices (  GFace *gf,
+				       GEdgeLoop  &gel , 
+				       std::vector<BDS_Point*> &result,
+				       SBoundingBox3d &bbox,
+				       BDS_Mesh *m,
+				       std::map<BDS_Point*,MVertex*> &recover_map, 
+				       int &count)
+{
+
+  std::list<GEdgeSigned > temp;
+  temp.insert( temp.begin(), gel.begin() ,gel.end() );
+  
+  GEdgeLoop::iter it = temp.begin();
+  GEdgeSigned & ges = *it;
+  if (ges.ge->getBeginVertex() == ges.ge->getEndVertex())
+    {
+      // do not put any undecidable loop at front
+      temp.erase(it);
+      temp.push_back(ges);
+    }
+
+  // when we get the seed 
+  // first we consider the minimum parameter for the seed.
+
+  it = temp.begin();
+
+  int _SIGN = -1;
+  int nbSeam = 0;
+  while (it != temp.end())
+    {
+      ges = *it;
+      
+      //      printf("GEdge %d %d \n",ges.getBeginVertex()->tag(),ges.getEndVertex()->tag());
+      
+      bool seam = (ges.ge->isSeam(gf));
+      if (seam)nbSeam++;
+      if (nbSeam == 2)_SIGN = 1;
+
+      std::vector<MVertex *> edgeLoop;
+      std::vector<BDS_Point *> edgeLoop_BDS;
+      
+      // we decide to reverse the ordering in case of closed edge connected to
+      // a seam edge
+      if (it != temp.begin() && (ges.ge->getBeginVertex() == ges.ge->getEndVertex()) && (nbSeam >= 1))
+	{
+	  ges._sign = -1;
+	}
+
+      if ( ges._sign == 1)
+	{
+	  edgeLoop.push_back(ges.ge->getBeginVertex()->mesh_vertices[0]);
+	  for (int i=0;i<ges.ge->mesh_vertices.size();i++)
+	    edgeLoop.push_back(ges.ge->mesh_vertices[i]);
+	}
+      else
+	{
+	  edgeLoop.push_back(ges.ge->getEndVertex()->mesh_vertices[0]);
+	  for (int i=ges.ge->mesh_vertices.size()-1;i>=0;i--)	    
+	    edgeLoop.push_back(ges.ge->mesh_vertices[i]);
+	}
+      
+      for (int i=0;i<edgeLoop.size();i++)	    
+	{
+	  MVertex *here     = edgeLoop[i];
+	  GEntity *ge       = here->onWhat();    
+	  BDS_Point *pp;
+	  double U,V;
+	  SPoint2 param;
+	  switch (ge->dim())
+	    {
+	    case 0:
+	      {
+		GVertex *gve = (GVertex*)ge;
+		//		SPoint2 param1 = gve->reparamOnFace(gf,1);
+		SPoint2 param = gve->reparamOnFace(gf,_SIGN);
+		// 		printf("BOUM %g %g -- %g %g\n",param1.x(),param1.y(),param2.x(),param2.y());
+// 		SPoint2 dif1 = param1-param;
+// 		SPoint2 dif2 = param2-param;
+// 		double d1 = dif1.x()*dif1.x() + dif1.y()*dif1.y();
+// 		double d2 = dif2.x()*dif2.x() + dif2.y()*dif2.y();
+// 		if (d1 > d2)param = param2;
+// 		else param = param1;
+		U = param.x();  
+		V = param.y();
+		pp = m->add_point ( count, U,V,gf );
+		m->add_geom (gve->tag(), 0);
+		BDS_GeomEntity *g = m->get_geom(gve->tag(),0);
+		pp->g = g;
+	      }
+	      break;
+	    case 1:
+	      {
+		GEdge *ged = (GEdge*)ge;
+		if (ge != ges.ge)
+		  printf("argh\n");	     
+		double u;
+		here->getParameter(0,u);
+// 		if (seam)
+// 		  {
+// 		    SPoint2  param1 = ged->reparamOnFace(gf,u,1);
+// 		    SPoint2  param2 = ged->reparamOnFace(gf,u,-1);
+// 		    //		    printf("%g %g -- %g %g\n",param1.x(),param1.y(),param2.x(),param2.y());
+// 		    SPoint2 dif1 = param1-param;
+// 		    SPoint2 dif2 = param2-param;
+// 		    double d1 = dif1.x()*dif1.x() + dif1.y()*dif1.y();
+// 		    double d2 = dif2.x()*dif2.x() + dif2.y()*dif2.y();
+// 		    if (d1 > d2)param = param2;
+// 		    else param = param1;
+// 		  }
+// 		else
+		  {
+		    param = ged->reparamOnFace(gf,u,_SIGN);
+		  }
+		U = param.x();
+		V = param.y();
+		pp = m->add_point ( count, U,V,gf );
+		m->add_geom (ged->tag(), 1);
+		BDS_GeomEntity *g = m->get_geom(ged->tag(),1);
+		pp->g = g;
+	      }
+	      break;
+	    default:
+	      {
+		throw;
+	      }
+	    }	    
+	  bbox += SPoint3(U,V,0);	  
+	  edgeLoop_BDS.push_back(pp);
+	  recover_map[pp] = here;
+	  count++;
+	}
+
+      
+      result.insert(result.end(),edgeLoop_BDS.begin(),edgeLoop_BDS.end());	    
+
+      ++it;
+    }
+
+//   for (int i=0;i<result.size();i++)
+//     {
+//       //      printf("point %d (%g %g) (%d,%d)\n",i,result[i]->u,result[i]->v,result[i]->g->classif_tag,result[i]->g->classif_degree);
+//     }
+
+}
+
+
+void gmsh2DMeshGeneratorPeriodic ( GFace *gf )
+{
+
+  std::map<BDS_Point*,MVertex*> recover_map;
+
+  Range<double> rangeU = gf->parBounds(0);
+  Range<double> rangeV = gf->parBounds(1);
+  
+  const double du = rangeU.high() -rangeU.low();
+  const double dv = rangeV.high() -rangeV.low();
+
+  const double LC2D = sqrt ( du*du + dv*dv ); 
+
+  //  printf("LC2D %g (%g,%g), (%g,%g)\n",LC2D,rangeU.high(),rangeU.low(),rangeV.high(),rangeV.low());
+
+  // Buid a BDS_Mesh structure that is convenient for doing the actual meshing procedure    
+  BDS_Mesh *m = new BDS_Mesh;
+  std::vector< std::vector<BDS_Point* > > edgeLoops_BDS;
+  SBoundingBox3d bbox;
+  int nbPointsTotal = 0;
+  {
+    for (std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin() ; it != gf->edgeLoops.end() ; it++)
+      {
+	std::vector<BDS_Point* > edgeLoop_BDS;
+	buildConsecutiveListOfVertices ( gf, *it , edgeLoop_BDS, bbox, m, recover_map , nbPointsTotal);
+	edgeLoops_BDS.push_back(edgeLoop_BDS);
+      }
+  }
+
+  // build a point record structure to create the initial mesh
+  DocRecord doc;  
+  doc.points =  (PointRecord*)malloc((nbPointsTotal+4) * sizeof(PointRecord));
+  int count = 0;
+  for ( int i=0;i<edgeLoops_BDS.size();i++)
+    {
+      std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
+      for ( int j=0;j<edgeLoop_BDS.size();j++)
+	{
+	  BDS_Point *pp = edgeLoop_BDS[j];
+	  const double U = pp->u;
+	  const double V = pp->v;
+	  double XX = CTX.mesh.rand_factor * LC2D *
+	    (double)rand() / (double)RAND_MAX;
+	  double YY = CTX.mesh.rand_factor * LC2D *
+	    (double)rand() / (double)RAND_MAX;
+	  doc.points[count].where.h = U + XX;
+	  doc.points[count].where.v = V + YY;
+	  doc.points[count].adjacent = NULL;
+	  doc.points[count].data = pp;      
+	  count++;	  
+	}
+    }
+  /// Increase the size of the bounding box by 20 %
+  /// add 4 points than encloses the domain
+  /// Use negative number to distinguish thos fake vertices
+  bbox *= 1.5;
+
+  MVertex *bb[4];
+  bb[0] = new MVertex ( bbox.min().x(), bbox.min().y(), 0,0,-1);
+  bb[1] = new MVertex ( bbox.min().x(), bbox.max().y(), 0,0,-2);
+  bb[2] = new MVertex ( bbox.max().x(), bbox.min().y(), 0,0,-3);
+  bb[3] = new MVertex ( bbox.max().x(), bbox.max().y(), 0,0,-4);
+  
+  for ( int ip = 0 ; ip<4 ; ip++ )
+    {
+      BDS_Point *pp = m->add_point ( -ip-1, bb[ip]->x(),bb[ip]->y(), gf);
+      m->add_geom (gf->tag(), 2);
+      BDS_GeomEntity *g = m->get_geom(gf->tag(),2);
+      pp->g = g;
+      doc.points[nbPointsTotal+ip].where.h  = bb[ip]->x();
+      doc.points[nbPointsTotal+ip].where.v  = bb[ip]->y();
+      doc.points[nbPointsTotal+ip].adjacent = 0;
+      doc.points[nbPointsTotal+ip].data = pp;
+    }
+
+  for ( int ip = 0 ; ip<4 ; ip++ ) delete bb[ip];
+  
+  /// Use "fast" inhouse recursive algo to generate the triangulation
+  /// At this stage the triangulation is not what we need
+  ///   -) It does not necessary recover the boundaries 
+  ///   -) It contains triangles outside the domain (the first edge loop is the outer one)
+  Msg(DEBUG1,"Meshing of the convex hull (%d points)",nbPointsTotal);
+  Make_Mesh_With_Points(&doc,doc.points,nbPointsTotal+4);
+  for(int i = 0; i < doc.numTriangles; i++) 
+    {
+      BDS_Point *p1 = (BDS_Point*)doc.points[doc.delaunay[i].t.a].data;
+      BDS_Point *p2 = (BDS_Point*)doc.points[doc.delaunay[i].t.b].data;
+      BDS_Point *p3 = (BDS_Point*)doc.points[doc.delaunay[i].t.c].data;
+      BDS_Face *t = m->add_triangle ( p1->iD,p2->iD,p3->iD);
+    }  
+  // Free stuff
+  free (doc.points);
+  free (doc.delaunay);
+ 
+  // Recover the boundary edges
+  // and compute characteristic lenghts using mesh edge spacing
+    
+  BDS_GeomEntity CLASS_F (1,2);
+  BDS_GeomEntity CLASS_E (1,1);
+  
+  for ( int i=0;i<edgeLoops_BDS.size();i++)
+    {
+      std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
+      for ( int j=0;j<edgeLoop_BDS.size();j++)
+	{
+	  BDS_Edge * e = m->recover_edge ( edgeLoop_BDS[j]->iD,edgeLoop_BDS[(j+1)%edgeLoop_BDS.size()]->iD);	  
+	  if (!e)Msg(GERROR,"impossible to recover the edge %d %d\n",edgeLoop_BDS[j]->iD,edgeLoop_BDS[(j+1)%edgeLoop_BDS.size()]->iD);
+	  else e->g = &CLASS_E;
+	}
+    }	  
+  
+  //  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
+  {
+    std::list<BDS_Edge*>::iterator ite = m->edges.begin();
+    while (ite != m->edges.end())
+      {
+	BDS_Edge *e = *ite;
+	if ( e->g  && e->numfaces () == 2)
+	  {
+	    BDS_Point *oface[2];
+	    e->oppositeof(oface);
+	    if (oface[0]->iD < 0) 
+	      {
+		recur_tag ( e->faces(1) , &CLASS_F); 
+		break;
+	      }
+	    else if (oface[1]->iD < 0) 
+	      {
+		recur_tag ( e->faces(0) , &CLASS_F); 
+		break;
+	      }
+	    }
+	++ite;
+      }
+  }
+
+  // delete useless stuff
+  {
+    std::list<BDS_Face*>::iterator itt = m->triangles.begin();
+    while (itt != m->triangles.end())
+      {
+	BDS_Face *t = *itt;
+	if (!t->g)
+	  {
+	    m->del_face (t);
+	  }
+	++itt;
+      }
+  }
+
+  m->cleanup();
+
+ 
+
+
+  {
+    std::list<BDS_Edge*>::iterator ite = m->edges.begin();
+    while (ite != m->edges.end())
+      {
+	BDS_Edge *e = *ite;
+	if (e->numfaces() == 0)
+	  m->del_edge(e);
+	else
+	  { 
+	    if (!e->g)
+	      e->g = &CLASS_F;
+	    if (!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree)e->p1->g = e->g;
+	    if (!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree)e->p2->g = e->g;
+	  }
+	++ite;
+      }
+  }
+  m->cleanup();
+  m->del_point(m->find_point(-1));
+  m->del_point(m->find_point(-2));
+  m->del_point(m->find_point(-3));
+  m->del_point(m->find_point(-4));
+    
+  // goto hhh;
+  // start mesh generation
+  
+  RefineMesh (gf,*m,15);
+  //  OptimizeMesh (gf,*m,2);
+
+ //    char name[245];
+//    sprintf(name,"s%d.pos",gf->tag());
+//    outputScalarField(m->triangles, name);
+  if (gf->meshAttributes.recombine)
+    {
+      m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
+    }
+  
+  // fill the small gmsh structures
+  
+  {
+    std::set<BDS_Point*,PointLessThan>::iterator itp =  m->points.begin(); 
+    while (itp != m->points.end())
+      {
+	BDS_Point *p = *itp;
+	if (recover_map.find(p) == recover_map.end())
+	  {
+	    MVertex *v = new MFaceVertex (p->X,p->Y,p->Z,gf,p->u,p->v);
+	    recover_map[p] = v;
+	    gf->mesh_vertices.push_back(v);
+	  }
+	++itp;
+      }
+  }
+  
+  {
+    std::list<BDS_Face*>::iterator itt = m->triangles.begin();
+    while (itt != m->triangles.end())
+      {
+	BDS_Face *t = *itt;
+	if (!t->deleted)
+	  {
+	    BDS_Point *n[4];
+	    t->getNodes(n);
+	    MVertex *v1 = recover_map[n[0]];
+	    MVertex *v2 = recover_map[n[1]];
+	    MVertex *v3 = recover_map[n[2]];
+	    if (!n[3])
+	      gf->triangles.push_back(new MTriangle (v1,v2,v3) );	
+	    else
+	      {
+		MVertex *v4 = recover_map[n[3]];
+		gf->quadrangles.push_back(new MQuadrangle (v1,v2,v3,v4) );	
+	      }
+	  }
+	++itt;
+      }
+  }
+  
+  
+  // delete the mesh
+  
+  delete m; 
+}
+
+
+
 void deMeshGFace :: operator() (GFace *gf) 
 {
   for (unsigned int i=0;i<gf->mesh_vertices.size();i++) delete gf->mesh_vertices[i];
@@ -1066,7 +1455,7 @@ void meshGFace :: operator() (GFace *gf)
   Msg(STATUS2, "Meshing surface %d", gf->tag());
 
   // TEST TEST 
-  if (gf->surfPeriodic(2)) return;
+  //  if (gf->surfPeriodic(2)) return;
 
   // destroy the mesh if it exists
   deMeshGFace dem;
@@ -1094,7 +1483,11 @@ void meshGFace :: operator() (GFace *gf)
 
   Msg(DEBUG1, "Generating the mesh");
   // mesh the face
-  gmsh2DMeshGenerator ( gf ) ;
+  //    gmsh2DMeshGeneratorPeriodic ( gf ) ;
+    //  else
+  if(!gf->surfPeriodic(0))
+    gmsh2DMeshGenerator ( gf ) ;
+
 
   Msg(DEBUG1, "type %d %d triangles generated, %d internal vertices",
       gf->geomType(),gf->triangles.size(),gf->mesh_vertices.size());