diff --git a/Geo/GEdgeLoop.h b/Geo/GEdgeLoop.h
index 972c5ffc9220107e686aac7e08fa1acc7d386633..04fc678242a27ded00dc5e2180e87974043d3111 100644
--- a/Geo/GEdgeLoop.h
+++ b/Geo/GEdgeLoop.h
@@ -52,6 +52,7 @@ public:
   inline citer begin () const {return loop.begin();}
   inline citer end   () const {return loop.end();}
   int  count (GEdge*) const;
+  int  count () const {return loop.size();}
 private:
   std::list< GEdgeSigned > loop ;
 };
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 8308ab9e1de80657a195909f6b033248d4b0d678..18270f149e5bceb5a0057e16cb018656c5aedecc 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-  // $Id: GModelIO_OCC.cpp,v 1.8 2006-11-16 21:14:10 remacle Exp $
+  // $Id: GModelIO_OCC.cpp,v 1.9 2006-11-21 23:52:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -502,9 +502,9 @@ void OCC_Internals :: loadSTEP (const char *fn)
   Standard_Integer nb = reader.NbRootsForTransfer();
   reader.TransferRoots (); 
   shape = reader.OneShape();  
-  BRepTools::Clean (shape);
-  buildLists();
-  HealGeometry();
+//    BRepTools::Clean (shape);
+//    buildLists();
+//    HealGeometry();
   BRepTools::Clean (shape);
 }
 
@@ -516,9 +516,9 @@ void OCC_Internals :: loadIGES (const char *fn)
   Standard_Integer nb = reader.NbRootsForTransfer();
   reader.TransferRoots (); 
   shape = reader.OneShape();  
-  BRepTools::Clean (shape);
-  buildLists();
-  HealGeometry();
+//   BRepTools::Clean (shape);
+//   buildLists();
+//   HealGeometry();
   BRepTools::Clean (shape);
 }
 
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 1fe1b794da907bb2eda4e0c983876ccabff0ddb7..3ed9a5c3d7fa72abd5ae714b5ad03c92cd12bc7b 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.8 2006-11-20 12:44:09 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.9 2006-11-21 23:52:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -25,6 +25,9 @@
 #include "OCCFace.h"
 
 #if defined(HAVE_OCC)
+#include "Geom_Ellipse.hxx"
+#include "Geom_Circle.hxx"
+#include "Geom_Line.hxx"
 
 OCCEdge::OCCEdge(GModel *model, TopoDS_Edge edge, int num, GVertex *v1, GVertex *v2)
   : GEdge(model, num, v1, v2), trimmed(0),c(edge)
@@ -153,19 +156,54 @@ double OCCEdge::parFromPoint(const SPoint3 &pt) const
 
 GEntity::GeomType OCCEdge::geomType() const
 {
-  return Unknown;
+  if (curve.IsNull())
+    {
+      if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Circle))
+	return Circle;
+      else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Line))
+	return Line;
+      else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Ellipse))
+	return Ellipse;
+      //   else if (occface->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
+      //     return Cone;
+      return Unknown;
+    }
+  else
+    {
+      if (curve->DynamicType() == STANDARD_TYPE(Geom_Circle))
+	return Circle;
+      else if (curve->DynamicType() == STANDARD_TYPE(Geom_Line))
+	return Line;
+      else if (curve->DynamicType() == STANDARD_TYPE(Geom_Ellipse))
+	return Ellipse;
+      //   else if (occface->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
+      //     return Cone;
+      return Unknown;
+    }
 }
 
 int OCCEdge::minimumMeshSegments () const
 {
-  return 2 ;
+  if(geomType() == Circle || geomType() == Ellipse)
+    return (int)(fabs(s1 - s0) *
+		 (double)CTX.mesh.min_circ_points / Pi) - 1;
+  else
+    return GEdge::minimumMeshSegments () ;
 }
 
 int OCCEdge::minimumDrawSegments () const
 {
-  return CTX.geom.circle_points;
+  int n = GEdge::minimumDrawSegments();
+
+  if(geomType() == Line)
+    return n;
+  else if(geomType() == Circle || geomType() == Ellipse)
+    return CTX.geom.circle_points;
+  else
+    return 10 * n;
 }
 
+
 double OCCEdge::curvature(double par) const 
 {
   return GEdge::curvature(par);
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 19a2be842d127fee5dd73dc2c54270ae7e40d9b4..f2a618ef3dee8c1ac09799b96c0b105935fa75a2 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.9 2006-11-20 12:44:09 remacle Exp $
+// $Id: OCCFace.cpp,v 1.10 2006-11-21 23:52:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -33,7 +33,7 @@
 #include "gp_Pln.hxx"
 
 OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape &emap)
-  : GFace(m, num), s(_s),_periodic(false)
+  : GFace(m, num), s(_s)
 {
   TopExp_Explorer exp0, exp01, exp1, exp2, exp3;
   for (exp2.Init (s, TopAbs_WIRE); exp2.More(); exp2.Next())
@@ -60,14 +60,16 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
 
       for (GEdgeLoop::citer it = el.begin() ; it != el.end() ; ++it)
 	{
-	  if(!it->ge->is3D())_periodic = true;
-	  if (el.count (it->ge) > 1)_periodic = true;
 	  l_edges.push_back(it->ge);
 	  l_dirs.push_back(it->_sign);
 	}
       
       edgeLoops.push_back(el);
     }
+  BRepAdaptor_Surface surface( s );
+  _periodic[0] = surface.IsUPeriodic();
+  _periodic[1] = surface.IsVPeriodic();
+// 	      surface.IsVPeriodic()
 
   Msg(INFO,"OCC Face %d with %d edges",num,l_edges.size());
   ShapeAnalysis::GetFaceUVBounds (s, umin, umax, vmin, vmax);
@@ -171,9 +173,12 @@ GEntity::GeomType OCCFace::geomType() const
     return Plane;
   else if (occface->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface))
     return Cylinder;
+//   else if (occface->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
+//     return Cone;
   return Unknown;
 }
 
+
 double OCCFace::curvature (const SPoint2 &param) const
 {
   BRepAdaptor_Surface sf(s, Standard_True);
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index 25a84fbe4165789d2bb525fedfff0230a5296d6f..5819e34bad90dee2431e4996381d2aeea24832df 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -33,7 +33,7 @@ class OCCFace : public GFace {
   TopoDS_Face s;
   Handle(Geom_Surface) occface;
   double umin, umax, vmin, vmax;
-  bool _periodic;
+  bool _periodic[2];
  public:
   OCCFace(GModel *m, TopoDS_Face s, int num, TopTools_IndexedMapOfShape &emap);
 
@@ -58,10 +58,13 @@ class OCCFace : public GFace {
   
   virtual bool continuous(int dim) const { return true; }
   virtual bool degenerate(int dim) const { return false; }
-  virtual double period(int dir) const {throw;}
+  virtual double period(int dir) const 
+  {
+    return (dir)? vmax-vmin:umax-umin;
+  }
   ModelType getNativeType() const { return OpenCascadeModel; }
   void * getNativePtr() const { return (void*)&s; }
-  virtual bool surfPeriodic(int dim) const {return _periodic;}
+  virtual bool surfPeriodic(int dim) const {return _periodic[dim];}
   virtual SPoint2 parFromPoint(const SPoint3 &) const;
   virtual double curvature(const SPoint2 &param) const;
 };
diff --git a/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index 0268d7210c3f822d1245d6465628b2cc57d6773a..ca054bebd18a8a5f8364d5bf53496297e41312de 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCVertex.cpp,v 1.4 2006-11-20 12:44:09 remacle Exp $
+// $Id: OCCVertex.cpp,v 1.5 2006-11-21 23:52:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -35,7 +35,7 @@ double max_surf_curvature ( const GVertex *gv, double x, double y, double z , co
     {
       SPoint2 par = gv->reparamOnFace((*it),1);
       double cc = (*it)->curvature ( par );
-      if (cc < 1.e2) curv = std::max(curv, cc );					
+      curv = std::max(curv, cc );					
       ++it;
     }  
   return curv;
@@ -96,4 +96,15 @@ SPoint2 OCCVertex::reparamOnFace ( GFace *gf , int dir) const
     }
 }
 
+double OCCVertex::prescribedMeshSizeAtVertex() const { 
+  SBoundingBox3d b = model()->bounds();
+  double lc     = 0.1 * norm( SVector3 ( b.max() , b.min() ) ) * CTX.mesh.lc_factor;
+  double lc_min = 0.004 * norm(SVector3 ( b.max() , b.min() ) ) * CTX.mesh.lc_factor;
+  double maxc = max_curvature_of_surfaces();
+  if (maxc !=0)       
+    lc = std::max(lc_min,std::min (lc,6.28/(CTX.mesh.min_circ_points*maxc)));
+  return lc;
+}
+
+
 #endif
diff --git a/Geo/OCCVertex.h b/Geo/OCCVertex.h
index a34f2b2eaca0ef6f65998115b811c269a815bcfc..2783057a96847e876fa32d10d12c9ea30f7d25e7 100644
--- a/Geo/OCCVertex.h
+++ b/Geo/OCCVertex.h
@@ -61,14 +61,7 @@ class OCCVertex : public GVertex {
   }
   ModelType getNativeType() const { return OpenCascadeModel; }
   void * getNativePtr() const { return (void*) &v; }
-  virtual double prescribedMeshSizeAtVertex() const { 
-    SBoundingBox3d b = model()->bounds();
-    double lc = 0.1*norm ( SVector3 ( b.max() , b.min() ) ) * CTX.mesh.lc_factor;
-    double maxc = max_curvature_of_surfaces();
-    if (maxc !=0)       
-      lc = std::min (lc,6.28/(CTX.mesh.min_circ_points*maxc));
-    return lc;
-  }
+  virtual double prescribedMeshSizeAtVertex() const;
   virtual SPoint2 reparamOnFace ( GFace *gf , int) const;
 };
 
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 323ad8afc0739776b64954165af53bd9605f84d7..0d1613138f79ef390ea5ed466b35411fbfd68e24 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.64 2006-11-16 21:14:10 remacle Exp $
+// $Id: BDS.cpp,v 1.65 2006-11-21 23:52:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -29,7 +29,7 @@
 
 bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Face *t);
 
-void outputScalarField(std::list < BDS_Face * >t, const char *iii)
+void outputScalarField(std::list < BDS_Face * >t, const char *iii, int param)
 {
   FILE *f = fopen(iii, "w");
   fprintf(f, "View \"scalar\" {\n");
@@ -38,10 +38,16 @@ void outputScalarField(std::list < BDS_Face * >t, const char *iii)
   while(tit != tite) {
     BDS_Point *pts[4];
     (*tit)->getNodes(pts);
-    fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
-            pts[0]->u, pts[0]->v, 0.0,
-            pts[1]->u, pts[1]->v, 0.0,
-            pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
+    if (param)
+      fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+	      pts[0]->u, pts[0]->v, 0.0,
+	      pts[1]->u, pts[1]->v, 0.0,
+	      pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
+    else
+      fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+	      pts[0]->X, pts[0]->Y, pts[0]->Z,
+	      pts[1]->X, pts[1]->Y, pts[1]->Z,
+	      pts[2]->X, pts[2]->Y, pts[2]->Z,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
     ++tit;
   }
   fprintf(f, "};\n");
@@ -236,7 +242,8 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2)
 	  BDS_Edge *eee = find_edge (num1, num2);
 	  if (!eee)
 	    {
-	      outputScalarField(triangles, "debug.pos");
+	      outputScalarField(triangles, "debugp.pos",1);
+	      outputScalarField(triangles, "debugr.pos",0);
 	      return 0;
 	    }
 	  return eee;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index fb165472f7b1bb2492dec2ba9236fd27224e485b..17a5cd38d0f4550a84957913953a549edaa77e9d 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -424,5 +424,5 @@ public:
   bool import_view(Post_View *view, const double tolerance);
 };
 
-void outputScalarField(std::list < BDS_Face * >t, const char *fn);
+void outputScalarField(std::list < BDS_Face * >t, const char *fn, int param);
 void recur_tag(BDS_Face * t, BDS_GeomEntity * g);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 95f58c506e47b4c3e68ee67291dc7f4e56b4d117..e805dd4fbb620dfbf71df7b3129974b8603231d0 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.16 2006-11-16 18:32:41 remacle Exp $
+// $Id: meshGEdge.cpp,v 1.17 2006-11-21 23:52:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -139,7 +139,6 @@ void deMeshGEdge :: operator() (GEdge *ge)
 void meshGEdge :: operator() (GEdge *ge) 
 {  
   if(ge->geomType() == GEntity::DiscreteCurve) return;
-  if(!ge->is3D()) return;
 
   // Send a messsage to the GMSH environment
   Msg(INFO, "Meshing curve %d", ge->tag());
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 8da54f41a6d6ba32ded8f890158a4b1efd7f65e7..66618dc2cef935cd561997c2bb6641fbc1f2811c 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.20 2006-11-20 12:44:09 remacle Exp $
+// $Id: meshGFace.cpp,v 1.21 2006-11-21 23:52:59 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -223,6 +223,13 @@ inline double computeEdgeLinearLength ( BDS_Point *p1, BDS_Point *p2)
   const double l = sqrt (dx*dx+dy*dy+dz*dz);
   return l;
 }
+inline double computeParametricEdgeLength ( BDS_Point *p1, BDS_Point *p2)
+{
+  const double dx = p1->u-p2->u;
+  const double dy = p1->v-p2->v;
+  const double l = sqrt (dx*dx+dy*dy);
+  return l;
+}
 
 double NewGetLc ( BDS_Edge *e  )
 {
@@ -420,6 +427,9 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 	  if (!(*it)->deleted)
 	    {
 	      double lone = NewGetLc ( *it);
+
+	      if (lone < 1.e-10 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-4) lone = 2;
+
 	      if ((*it)->numfaces() == 2 && (lone >  1.3))
 		{
 		  BDS_Point *mid ;
@@ -464,6 +474,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 	{
 	  if (NN2++ >= NN1)break;
 	  double lone = NewGetLc ( *it);
+	  if (lone < 1.e-10 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-4) lone = 2;
 	  if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < 0.6 )
 	    {
 	      bool res = false;
@@ -1044,6 +1055,19 @@ void gmsh2DMeshGenerator ( GFace *gf )
   //  std::for_each(gf->mesh_vertices.begin(),gf->mesh_vertices.end(),p2c);    
 }
 
+// this function buils a list of vertices (BDS) that 
+// are consecutive in one given edge loop. We take
+// care of periodic surfaces. In the case of periodicity, some
+// curves are present 2 times in the wire (seams). Those
+// must be meshed separately
+
+inline double dist2 (const SPoint2 &p1,const SPoint2 &p2)
+{
+  const double dx = p1.x() - p2.x(); 
+  const double dy = p1.y() - p2.y(); 
+  return dx*dx+dy*dy;
+}
+
 void buildConsecutiveListOfVertices (  GFace *gf,
 				       GEdgeLoop  &gel , 
 				       std::vector<BDS_Point*> &result,
@@ -1053,139 +1077,181 @@ void buildConsecutiveListOfVertices (  GFace *gf,
 				       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);
-    }
+  // for each edge, we build a list of points that 
+  // are the mapping of the edge points on the face
+  // for seams, we build the list for every side
+  // for closed loops, we build it on both senses
 
-  // 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++;
-	}
+  //  printf("new edge loop\n");
+  std::map<GEntity*,std::vector<SPoint2> > meshes;
+  std::map<GEntity*,std::vector<SPoint2> > meshes_seam;
+  
+  GEdgeLoop::iter it  = gel.begin();  
 
+  while (it != gel.end())   
+   {
+     GEdgeSigned ges = *it ;
       
-      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);
-//     }
+     std::vector<SPoint2> mesh1d;
+     std::vector<SPoint2> mesh1d_seam;
+
+     bool seam = ges.ge->isSeam(gf);
+     Range<double> range = ges.ge->parBounds(0);
+
+     MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0];
+     mesh1d.push_back(ges.ge->reparamOnFace(gf,range.low(),1));
+     if ( seam ) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf,range.low(),-1));
+     for (int i=0;i<ges.ge->mesh_vertices.size();i++)
+       {
+	 double u;
+	 here = ges.ge->mesh_vertices[i];
+	 here->getParameter(0,u);
+	 mesh1d.push_back(ges.ge->reparamOnFace(gf,u,1));
+	 if ( seam ) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf,u,-1));
+       }
+     here = ges.ge->getEndVertex()->mesh_vertices[0];
+     mesh1d.push_back(ges.ge->reparamOnFace(gf,range.high(),1));
+     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) );
+
+     
+     it++;
+   }
+
+
+  std::list<GEdgeSigned> unordered;
+  unordered.insert(unordered.begin(),gel.begin(),gel.end());
+  
+  GEdgeSigned found(0,0);
+
+  SPoint2 last_coord;  
+  int counter = 0;
+
+  while (unordered.size())
+   {
+     std::list<GEdgeSigned>::iterator it = unordered.begin();     
+     std::vector<SPoint2>  coords;
+
+     while (it != unordered.end())   
+       {
+	 std::vector<SPoint2>  mesh1d;
+	 std::vector<SPoint2>  mesh1d_seam;
+	 std::vector<SPoint2>  mesh1d_reversed;
+	 std::vector<SPoint2>  mesh1d_seam_reversed;
+	 GEdge *ge = (*it).ge;
+	 bool seam = ge->isSeam(gf);
+	 mesh1d = meshes[ge];
+	 if (seam){mesh1d_seam = meshes_seam[ge];}
+	 mesh1d_reversed.insert(mesh1d_reversed.begin(),mesh1d.rbegin(),mesh1d.rend());
+	 if (seam)mesh1d_seam_reversed.insert(mesh1d_seam_reversed.begin(),mesh1d_seam.rbegin(),mesh1d_seam.rend());
+	 if (!counter)
+	   {
+	     counter++;
+	     coords = ((*it)._sign == 1)?mesh1d:mesh1d_reversed;	 
+	     found = (*it);
+	     unordered.erase(it);
+	     break;
+	   }
+	 else
+	   {
+	     SPoint2 first_coord         = mesh1d[0];
+	     double d = dist2(last_coord,first_coord);
+	     //	     printf("d = %12.5E %d\n",d, coords.size());
+	     if (d < 1.e-8) 
+	       {
+		 coords.clear();
+		 coords = mesh1d;
+		 found = GEdgeSigned(1,ge);
+		 unordered.erase(it);
+		 break;
+	       }
+	     SPoint2 first_coord_reversed = mesh1d_reversed[0];
+	     double d_reversed = dist2(last_coord,first_coord_reversed);
+	     //	     printf("d_r = %12.5E\n",d_reversed);
+	     if (d_reversed < 1.e-8) 
+	       {
+		 coords.clear();
+		 coords = mesh1d_reversed;
+		 found = (GEdgeSigned(-1,ge));
+		 unordered.erase(it);
+		 break;
+	       }
+	     if (seam)
+	       {
+		 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);
+		 //		 printf("d_seam = %12.5E\n",d_seam);
+		 if (d_seam < 1.e-8)
+		   {
+		     coords.clear();
+		     coords = mesh1d_seam;
+		     found = (GEdgeSigned(1,ge));
+		     unordered.erase(it);
+		     break;
+		   }
+		 double d_seam_reversed = dist2(last_coord,first_coord_seam_reversed);
+		 //		 printf("d_seam_reversed = %12.5E\n",d_seam_reversed);
+		 if (d_seam_reversed < 1.e-8)
+		   {
+		     coords.clear();
+		     coords = mesh1d_seam_reversed;
+		     found = (GEdgeSigned(-1,ge));
+		     unordered.erase(it);
+		     break;
+		   }
+	       }
+	   }
+	 ++it;
+       }
+
+     std::vector<MVertex*>    edgeLoop;
+     if ( found._sign == 1)
+       {
+	 edgeLoop.push_back(found.ge->getBeginVertex()->mesh_vertices[0]);
+	 for (int i=0;i<found.ge->mesh_vertices.size();i++)
+	   edgeLoop.push_back(found.ge->mesh_vertices[i]);
+       }
+     else
+       {
+	 edgeLoop.push_back(found.ge->getEndVertex()->mesh_vertices[0]);
+	 for (int i=found.ge->mesh_vertices.size()-1;i>=0;i--)	    
+	   edgeLoop.push_back(found.ge->mesh_vertices[i]);
+       }
+     
+     //     printf("edge %d size %d size %d\n",found.ge->tag(),edgeLoop.size(), coords.size());
+     
+     std::vector<BDS_Point*>  edgeLoop_BDS;
+     for (int i=0;i<edgeLoop.size();i++)	    
+       {
+	 MVertex *here     = edgeLoop[i];
+	 GEntity *ge       = here->onWhat();    
+	 double U,V;
+	 SPoint2 param = coords [i];
+	 U = param.x();
+	 V = param.y();
+	 BDS_Point *pp;
+	 pp = m->add_point ( count, U,V,gf );
+	 m->add_geom (ge->tag(), ge->dim());
+	 BDS_GeomEntity *g = m->get_geom(ge->tag(),ge->dim());
+	 pp->g = g;
+	 //	 printf("point %3d (%8.5f %8.5f) (%2d,%2d)\n",count,pp->u,pp->v,pp->g->classif_tag,pp->g->classif_degree);
+	 bbox += SPoint3(U,V,0);	  
+	 edgeLoop_BDS.push_back(pp);
+	 recover_map[pp] = here;	 
+	 count++;
+       }     
+     last_coord = coords[coords.size()-1];
+     //     printf("last coord %g %g\n",last_coord.x(),last_coord.y());
+     result.insert(result.end(),edgeLoop_BDS.begin(),edgeLoop_BDS.end());	         
+   }
+
+//    for (int i=0;i<result.size();i++)
+//      {
+//        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);
+//      }
 
 }
 
@@ -1290,6 +1356,7 @@ void gmsh2DMeshGeneratorPeriodic ( GFace *gf )
   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];
@@ -1379,9 +1446,8 @@ void gmsh2DMeshGeneratorPeriodic ( GFace *gf )
   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);
@@ -1431,6 +1497,11 @@ void gmsh2DMeshGeneratorPeriodic ( GFace *gf )
   
   // delete the mesh
   
+//   char name[245];
+//   sprintf(name,"param%d.pos",gf->tag());
+//   outputScalarField(m->triangles, name,1);
+//   sprintf(name,"real%d.pos",gf->tag());
+//   outputScalarField(m->triangles, name,0);
   delete m; 
 }
 
@@ -1470,7 +1541,7 @@ void meshGFace :: operator() (GFace *gf)
   std::vector<MVertex*> points;
   std::vector<int> indices;
   // compute loops on the fly
-  // indices indicate start and end points of a loop
+// indices indicate start and end points of a loop
   // loops are not yet oriented
   Msg(DEBUG1, "Computing edge loops");
   computeEdgeLoops(gf, points, indices);
@@ -1483,10 +1554,9 @@ void meshGFace :: operator() (GFace *gf)
 
   Msg(DEBUG1, "Generating the mesh");
   // mesh the face
-  //    gmsh2DMeshGeneratorPeriodic ( gf ) ;
-    //  else
-  if(!gf->surfPeriodic(0))
-    gmsh2DMeshGenerator ( gf ) ;
+  gmsh2DMeshGeneratorPeriodic ( gf ) ;
+      //  else
+  //gmsh2DMeshGenerator ( gf ) ;
 
 
   Msg(DEBUG1, "type %d %d triangles generated, %d internal vertices",