diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 5bc122a1ed944f673161f4b275d314e6819beb3f..a70fb278b68cc26e6fe3ba66fd021989c184938c 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -243,6 +243,7 @@ class GEntity {
 
   // Vertex arrays to draw the mesh efficiently
   VertexArray *va_lines, *va_triangles;
+
 };
 
 class GEntityLessThan {
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index ec5ee8fc70f0b4c74c50ef739fb9428639be324d..4782cf9f194913bad71acd9269e5df1391ba049c 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.42 2008-01-22 16:47:10 geuzaine Exp $
+// $Id: GFace.cpp,v 1.43 2008-02-05 14:40:29 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -23,6 +23,7 @@
 #include "GFace.h"
 #include "GEdge.h"
 #include "MElement.h"
+#include "VertexArray.h"
 
 #if defined(HAVE_GMSH_EMBEDDED)
 #  include "GmshEmbedded.h"
@@ -42,7 +43,7 @@ void dsvdcmp(double **a, int m, int n, double w[], double **v);
 
 extern Context_T CTX;
 
-GFace::GFace(GModel *model, int tag) : GEntity(model, tag), r1(0), r2(0) 
+GFace::GFace(GModel *model, int tag) : GEntity(model, tag), r1(0), r2(0), va_geom_triangles(0)
 {
   meshStatistics.status = GFace::PENDING;
   resetMeshAttributes();
@@ -65,6 +66,9 @@ GFace::~GFace()
 
   for(unsigned int i = 0; i < quadrangles.size(); i++) 
     delete quadrangles[i];
+
+  if (va_geom_triangles)
+    delete va_geom_triangles;
 }
 
 void GFace::resetMeshAttributes()
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 325c723a51129a64891d783e9b394317fcd2a5dc..34399563e027fcbdb4e65a1ae1d961d198b121d8 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -73,7 +73,9 @@ class GFace : public GEntity
   // i.e. closed edge loops.  the first wire is the one that is the
   // outer contour of the face.
   void resolveWires();
-
+  // builds a STL triangulation and fill the vertex array va_geom_triangles
+  // by default, do nothing
+  virtual bool buildSTLTriangulation () {return false;}
  public:
   GFace(GModel *model, int tag);
   virtual ~GFace();
@@ -205,6 +207,10 @@ class GFace : public GEntity
 
   std::vector<MTriangle*> triangles;
   std::vector<MQuadrangle*> quadrangles;
+
+  // Vertex array to draw the geometry efficiently
+  VertexArray *va_geom_triangles;
+
 };
 
 #endif
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index facabba2153b29f9d87293051196ab6888f85e9c..8342dd2ef1601f1a416e864a6d05a2b62dbbee90 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_OCC.cpp,v 1.24 2008-01-19 22:06:01 geuzaine Exp $
+// $Id: GModelIO_OCC.cpp,v 1.25 2008-02-05 14:40:29 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -496,6 +496,122 @@ void GModel::deleteOCCInternals()
   if(occ_internals) delete occ_internals;
 }
 
+
+/*
+  OCC Creation routines
+*/
+
+// This function has been inspired from SALOME
+// It removes all duplicates from the geometry, starting
+// from vertices, edges, faces, shells and solids
+// This 
+void OCC_Internals::removeAllDuplicates (const double &tolerance){
+}
+
+
+void AddSimpleShapes(TopoDS_Shape theShape, TopTools_ListOfShape& theList)
+{
+  if (theShape.ShapeType() != TopAbs_COMPOUND &&
+      theShape.ShapeType() != TopAbs_COMPSOLID) {
+    theList.Append(theShape);
+    return;
+  }
+
+  TopTools_MapOfShape mapShape;
+  TopoDS_Iterator It (theShape, Standard_True, Standard_True);
+
+  for (; It.More(); It.Next()) {
+    TopoDS_Shape aShape_i = It.Value();
+    if (mapShape.Add(aShape_i)) {
+      if (aShape_i.ShapeType() == TopAbs_COMPOUND ||
+          aShape_i.ShapeType() == TopAbs_COMPSOLID) {
+        AddSimpleShapes(aShape_i, theList);
+      } else {
+        theList.Append(aShape_i);
+      }
+    }
+  }
+}
+
+void OCC_Internals::applyBooleanOperator ( TopoDS_Shape tool ,  const BooleanOperator & op ){
+  if (tool.IsNull())return;
+  if (shape.IsNull())shape = tool;
+  else{
+    switch(op){
+    case OCC_Internals::Add :
+      {
+	TopoDS_Shape theNewShape;	
+	BRep_Builder B;
+	TopoDS_Compound C;
+	B.MakeCompound(C);
+	TopTools_ListOfShape listShape1, listShape2;
+	AddSimpleShapes(shape, listShape1);
+	AddSimpleShapes(tool, listShape2);
+	Standard_Boolean isCompound =
+	  (listShape1.Extent() > 1 || listShape2.Extent() > 1);
+	
+	TopTools_ListIteratorOfListOfShape itSub1 (listShape1);
+	for (; itSub1.More(); itSub1.Next()) {
+	  TopoDS_Shape aValue1 = itSub1.Value();
+	  TopTools_ListIteratorOfListOfShape itSub2 (listShape2);
+	  for (; itSub2.More(); itSub2.Next()) {
+	    TopoDS_Shape aValue2 = itSub2.Value();
+	    BRepAlgoAPI_Common BO (aValue1, aValue2);
+	    if (!BO.IsDone()) {
+	      Msg(GERROR,"Boolean Add Operator can not be performed");
+	    }
+	    if (isCompound) {
+	      TopoDS_Shape aStepResult = BO.Shape();
+	      if (aStepResult.ShapeType() == TopAbs_COMPOUND) {
+		TopoDS_Iterator aCompIter (aStepResult);
+		for (; aCompIter.More(); aCompIter.Next()) {
+		  B.Add(C, aCompIter.Value());
+		}
+	      }
+	      else {
+		B.Add(C, aStepResult);
+	      }
+	    }
+	    else
+	      theNewShape = BO.Shape();
+	  }
+	}
+	if (isCompound) {
+	  TopTools_ListOfShape listShapeC;
+	  AddSimpleShapes(C, listShapeC);
+	  TopTools_ListIteratorOfListOfShape itSubC (listShapeC);
+	  bool isOnlySolids = true;
+	  for (; itSubC.More(); itSubC.Next()) {
+	    TopoDS_Shape aValueC = itSubC.Value();
+	    if (aValueC.ShapeType() != TopAbs_SOLID) isOnlySolids = false;
+	  }
+	  if (isOnlySolids)
+	    throw;
+	    //	    theNewShape = GEOMImpl_GlueDriver::GlueFaces(C, Precision::Confusion());
+	  else
+	    theNewShape = C;
+	}	
+      }
+      break;
+    case OCC_Internals::Cut :
+      {
+      }
+    default :
+      throw;
+    }
+  }
+}
+  
+void OCC_Internals::Sphere  ( const SPoint3 & center, const double & radius, const BooleanOperator & op ){
+  // build a sphere
+  gp_Pnt aP (center.x(), center.y(), center.z());  
+  TopoDS_Shape aShape = BRepPrimAPI_MakeSphere(aP, radius).Shape(); 
+  // either add it to the current shape, or use it as a tool and remove the
+  // sphere from the current shape
+  applyBooleanOperator ( aShape , op );
+}
+
+
 #else
 
 int GModel::readOCCSTEP(const std::string &fn)
diff --git a/Geo/GModelIO_OCC.h b/Geo/GModelIO_OCC.h
index 3c51437dca8756a56a0dfb87e43616c3b64f1607..bdc4569d8c593a70abc950aa14a4448a6635f946 100644
--- a/Geo/GModelIO_OCC.h
+++ b/Geo/GModelIO_OCC.h
@@ -1,7 +1,7 @@
 #ifndef _GMODELIO_OCC_H_
 #define _GMODELIO_OCC_H_
 
-// $Id: GModelIO_OCC.h,v 1.2 2007-04-23 08:04:16 geuzaine Exp $
+// $Id: GModelIO_OCC.h,v 1.3 2008-02-05 14:40:29 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -32,6 +32,9 @@ class OCC_Internals {
   TopoDS_Shape shape;
   TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;
  public:
+  
+  enum BooleanOperator { Add , Cut }; 
+
   OCC_Internals()
   {
     somap.Clear();
@@ -49,6 +52,11 @@ class OCC_Internals {
   void loadBREP(const char *);  
   void buildGModel(GModel *gm);
   void buildLists();
+  void removeAllDuplicates (const double &tolerance);
+
+  void Sphere  ( const SPoint3 & center, const double & radius, const BooleanOperator & op );
+  void Cylinder( const SPoint3 & bottom_center, const SVector3 & dir, const BooleanOperator & op );
+  void applyBooleanOperator ( TopoDS_Shape tool, const BooleanOperator & op);
 };
 
 #endif
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 347abe9fdb629852da68a44480072fb67aa2c2b0..9903385554cf614fb8eeec3c6ca4e0457289008e 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.28 2008-02-03 08:54:28 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.29 2008-02-05 14:40:29 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -26,6 +26,7 @@
 #include "OCCFace.h"
 #include "Message.h"
 #include "Numeric.h"
+#include "VertexArray.h"
 
 #if defined(HAVE_OCC)
 #include "Geom_CylindricalSurface.hxx"
@@ -37,6 +38,7 @@
 #include "Geom_BezierSurface.hxx"
 #include "Geom_Plane.hxx"
 #include "gp_Pln.hxx"
+#include "BRepMesh_FastDiscret.hxx"
 
 OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape &emap)
   : GFace(m, num), s(_s)
@@ -88,6 +90,7 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
   umax += fabs(du)/100.0;
   vmax += fabs(dv)/100.0;
   occface = BRep_Tool::Surface(s);
+  buildSTLTriangulation();
 }
 
 Range<double> OCCFace::parBounds(int i) const
@@ -272,14 +275,79 @@ surface_params OCCFace::getSurfaceParams() const
 }
 
 
-// void OCCFace::buildVisTriangulation ();
-// {
-//   TopLoc_Location loc;
-//   Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (occface, loc);
-//   int ntriangles = triangulation -> NbTriangles();
-//   for (int j = 1; j <= ntriangles; j++){
-//     Poly_Triangle triangle = (triangulation -> Triangles())(j);
-//   }  
-// }
+bool OCCFace::buildSTLTriangulation ()
+{
+
+  if (va_geom_triangles){
+    delete va_geom_triangles;
+    va_geom_triangles = 0;
+  }
+
+   TopLoc_Location loc;
+   int p1,p2,p3;
+   Bnd_Box aBox;
+   Standard_Boolean bWithShare;
+   Standard_Real aDiscret, aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
+   Standard_Real dX, dY, dZ, dMax, aCoeff, aAngle;     
+   BRepBndLib::Add(s, aBox);
+   aBox.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
+   
+   //   if (tag() > 100)return;
+   dX=aXmax-aXmin;
+   dY=aYmax-aYmin;
+   dZ=aZmax-aZmin;
+   dMax=dX;
+   if (dY>dMax) {
+     dMax=dY;
+   }
+   if (dZ>dMax) {
+     dMax=dZ;
+   }
+   //
+   aCoeff=0.01;
+   aDiscret=aCoeff*dMax;     
+   BRepMesh_FastDiscret aMesher(aDiscret,
+				0.5,
+				aBox,
+				bWithShare,
+				Standard_True,
+				Standard_False,
+				Standard_True);
+   aMesher.Add(s);
+   Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (s, loc);
+   if (triangulation.IsNull()){
+     Msg(WARNING,"OCC STL triangulation failed!\n");    //     return;
+     return false;
+   }
+
+   int ntriangles = triangulation -> NbTriangles();
+   if (!triangulation->HasUVNodes())
+     return false;
+   va_geom_triangles  = new VertexArray(3,ntriangles);
+
+
+   unsigned int c = CTX.color.geom.surface;
+   unsigned int col[4] = {c, c, c, c};
+   for (int j = 1; j <= ntriangles; j++){
+     Poly_Triangle triangle = (triangulation -> Triangles())(j);
+     triangle.Get (p1,p2,p3);
+     gp_Pnt2d x1 = (triangulation->UVNodes())(p1);
+     gp_Pnt2d x2 = (triangulation->UVNodes())(p2);
+     gp_Pnt2d x3 = (triangulation->UVNodes())(p3);
+     GPoint gp1 = point (x1.X(),x1.Y());
+     GPoint gp2 = point (x2.X(),x2.Y());
+     GPoint gp3 = point (x3.X(),x3.Y());
+     SVector3 n[3];
+     n[0] = normal (SPoint2(x1.X(),x1.Y()));
+     n[1] = normal (SPoint2(x2.X(),x2.Y()));
+     n[2] = normal (SPoint2(x3.X(),x3.Y()));
+     double x[3] = {gp1.x(),gp2.x(),gp3.x()};
+     double y[3] = {gp1.y(),gp2.y(),gp3.y()};
+     double z[3] = {gp1.z(),gp2.z(),gp3.z()};
+     va_geom_triangles->add (x,y,z,n,col);
+   }  
+   va_geom_triangles->finalize();
+   return true;
+}
 
 #endif
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index b686dc66a892bbf367016370e0acf985047dc479..f91db45cd3b22720d578b518bf5e65fdd9c4a4eb 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -34,6 +34,7 @@ class OCCFace : public GFace {
   Handle(Geom_Surface) occface;
   double umin, umax, vmin, vmax;
   bool _periodic[2];
+  bool buildSTLTriangulation ();
  public:
   OCCFace(GModel *m, TopoDS_Face s, int num, TopTools_IndexedMapOfShape &emap);
   virtual ~OCCFace(){}
diff --git a/Geo/OCCIncludes.h b/Geo/OCCIncludes.h
index 7557292f0403565f880aaaeb8efbefeb919fe3a3..a8b91c7f017d42959ec2ebaf5d305a14a56b1a7a 100644
--- a/Geo/OCCIncludes.h
+++ b/Geo/OCCIncludes.h
@@ -109,7 +109,13 @@ using std::iostream;
 #include "ShapeBuild_ReShape.hxx"
 #include "ShapeFix.hxx"
 #include "ShapeFix_FixSmallFace.hxx"
-
+#include "TopoDS_Compound.hxx"
+#include "TopoDS_Iterator.hxx"
+#include "BRepPrimAPI_MakeSphere.hxx"
+#include "TopTools_ListIteratorOfListOfShape.hxx"
+#include "Precision.hxx"
+#include "BRepAlgoAPI_Common.hxx"
+#include "BRepAlgoAPI_Cut.hxx"
 #endif
 
 #endif
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 321a5a99522ce680be72d499cba00060c9267785..4be7623dfa96766b3a723c725e107ba6bcd95ba1 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.145 2008-02-03 08:25:36 geuzaine Exp $
+// $Id: Geom.cpp,v 1.146 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -200,6 +200,8 @@ class drawGEdge {
   }
 };
 
+void drawArrays(GEntity *e, VertexArray *va, GLint type, bool useNormalArray, 
+		int forceColor=0, unsigned int color=0);
 class drawGFace {
  private:
   void _drawNonPlaneGFace(GFace *f)
@@ -444,7 +446,9 @@ public :
       glColor4ubv((GLubyte *) & CTX.color.geom.surface);
     }
 
-    if(f->geomType() == GEntity::Plane)
+    if(f->va_geom_triangles)
+      drawArrays(f, f->va_geom_triangles, GL_TRIANGLES, CTX.geom.light, CTX.geom.surfaces,CTX.color.geom.surface);
+    else if(f->geomType() == GEntity::Plane)
       _drawPlaneGFace(f);
     else if(f->geomType() == GEntity::ProjectionFace ||
 	    f->geomType() == GEntity::ParametricSurface)
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 1f11c35c37034f9404cc06c247b1ac2741dc3731..80ff495a6ab8472b0eb3f4394caab73cf00c2376 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.211 2008-01-19 22:06:02 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.212 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -440,7 +440,7 @@ static void addElementsInArrays(GEntity *e, std::vector<T*> &elements,
   }
 }
 
-static void drawArrays(GEntity *e, VertexArray *va, GLint type, bool useNormalArray, 
+void drawArrays(GEntity *e, VertexArray *va, GLint type, bool useNormalArray, 
 		       int forceColor=0, unsigned int color=0)
 {
   if(!va) return;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 603ecf9e1c0f65bc5a519bd7f0d9d2eec08316e5..85744518732ae067168057fad05ada7e3e7633e0 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.96 2008-01-30 15:27:41 remacle Exp $
+// $Id: BDS.cpp,v 1.97 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -30,7 +30,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, int param)
+void outputScalarField(std::list < BDS_Face * >t, const char *iii, int param, GFace *gf)
 {
   FILE *f = fopen(iii, "w");
   if(!f) return;
@@ -45,11 +45,21 @@ void outputScalarField(std::list < BDS_Face * >t, const char *iii, int param)
 	      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);
+    else{
+      if (!gf)
+	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);
+      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,
+		gf->curvature (SPoint2(pts[0]->u, pts[0]->v)),
+		gf->curvature (SPoint2(pts[1]->u, pts[1]->v)),
+		gf->curvature (SPoint2(pts[2]->u, pts[2]->v)));
+    }
     ++tit;
   }
   fprintf(f, "};\n");
@@ -182,6 +192,26 @@ BDS_Edge *BDS_Mesh::find_edge(int num1, int num2)
 int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
 		       double x3, double y3, double x4, double y4)
 {
+
+//   double p1[2] = {x1,y1};
+//   double p2[2] = {x2,y2};
+//   double q1[2] = {x3,y3};
+//   double q2[2] = {x4,y4};
+  
+//   double rp1 = gmsh::orient2d(p1, p2, q1);
+//   double rp2 = gmsh::orient2d(p1, p2, q2);
+//   double rq1 = gmsh::orient2d(q1, q2, p1);
+//   double rq2 = gmsh::orient2d(q1, q2, p2);
+
+//   if (rp1*rp2<=0 && rq1*rq2<=0){
+// //      printf("p1 %22.15E %22.15E %22.15E %22.15E \n",x1,y1,x2,y2);
+// //      printf("p2 %22.15E %22.15E %22.15E %22.15E \n",x3,y3,x4,y4);
+// //      printf("or %22.15E %22.15E %22.15E %22.15E\n",rp1,rp2,rq1,rq2);
+//     return 1;
+//   }
+//   return 0;
+
+
   double mat[2][2];
   double rhs[2], x[2];
   mat[0][0] = (x2 - x1);
@@ -235,7 +265,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
 		  {
 		    std::set<EdgeToRecover>::iterator itr1 = e2r->find(EdgeToRecover(e->p1->iD,e->p2->iD,0));		    
 		    std::set<EdgeToRecover>::iterator itr2 = e2r->find(EdgeToRecover(num1,num2,0));		    
-		    Msg(WARNING," edge %d %d on model edge %d cannot be recovered because it intersects %d %d on model edge %d",
+		    Msg(DEBUG2," edge %d %d on model edge %d cannot be recovered because it intersects %d %d on model edge %d",
 			num1,num2,itr2->ge->tag(),
 			e->p1->iD,e->p2->iD,itr1->ge->tag());
 
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 72847258181073ca3a0f8492e4e77b37b5df3b5b..c8882ae66cf363d09dba0110d942e42369994dda 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -22,6 +22,9 @@
 // points may know the normals to the surface they are classified on
 // default values are 0,0,1
 
+#ifndef _BDS_GMSH_H_
+#define _BDS_GMSH_H_
+
 #include <string>
 #include <set>
 #include <map>
@@ -474,5 +477,7 @@ public:
   void recombineIntoQuads (const double angle, GFace *gf);
 };
 
-void outputScalarField(std::list < BDS_Face * >t, const char *fn, int param);
+void outputScalarField(std::list < BDS_Face * >t, const char *fn, int param, GFace *gf = 0);
 void recur_tag(BDS_Face * t, BDS_GeomEntity * g);
+
+#endif
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index c9f60b13126dada1e70ed11ce90bdcf97c1eba05..89717c853005a28f7adfdce6451bc40267ff30a5 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.33 2008-01-30 15:27:41 remacle Exp $
+// $Id: BackgroundMesh.cpp,v 1.34 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -79,6 +79,7 @@ static double max_surf_curvature(const GEdge *ge, double u)
     SPoint2 par = ge->reparamOnFace((*it),u,1);
     double cc = (*it)->curvature(par);
     max_curvature = std::max(cc, max_curvature);
+    //    if (ge->tag() == 66)printf("%12.5E %12.5E %d %12.5E %12.5E\n",par.x(),par.y(),(*it)->tag(),cc,max_curvature);
     ++it;
   }  
   return max_curvature;
@@ -195,10 +196,15 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
 
   if(CTX.mesh.lc_from_curvature && ge->dim() <= 2 )
     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 (ge->tag() == 200 || ge->tag() == 202|| ge->tag() == 162|| ge->tag() == 163|| ge->tag() == 161|| ge->tag() == 141){
+//     printf("%d %d %12.5E %12.5E %12.5E\n",ge->dim(),ge->tag(),LC_MVertex_CURV(ge, U, V),LC_MVertex_PNTS(ge, U, V),lc);
+//   }
+
   if(lc <= 0.){
     Msg(GERROR, "Incorrect char. length lc = %g: using default instead", lc);
     return l3 * CTX.mesh.lc_factor;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 39496c40544e846d8835606a86c92403456ef369..84f7f0f892238546f1a2bf6a4186059f3e876505 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.133 2008-01-28 16:03:19 geuzaine Exp $
+// $Id: Generator.cpp,v 1.134 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -28,6 +28,7 @@
 #include "MElement.h"
 #include "meshGEdge.h"
 #include "meshGFace.h"
+#include "meshGFaceBDS.h"
 #include "meshGRegion.h"
 #include "BackgroundMesh.h"
 #include "BoundaryLayer.h"
@@ -261,6 +262,8 @@ void Mesh2D(GModel *m)
     }
   }
 
+  gmshCollapseSmallEdges (*m);
+
   double t2 = Cpu();
   CTX.mesh_timer[1] = t2 - t1;
   Msg(INFO, "Mesh 2D complete (%g s)", CTX.mesh_timer[1]);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 5d618b1be9dbdef7e61e66e796aa8443998611da..98922988219d907e45aa2cc1b23f44be9b0bc094 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.52 2008-01-30 15:27:41 remacle Exp $
+// $Id: meshGEdge.cpp,v 1.53 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -322,7 +322,7 @@ void meshGEdge::operator() (GEdge *ge)
 		  CTX.mesh.lc_integration_precision);
       buildInterpLc(lcPoints);
       //      printInterpLc("toto1.dat");
-      smoothInterpLc(20);
+      //      smoothInterpLc(20);
       //      printInterpLc("toto2.dat");
       a = Integration(ge, t_begin, t_end, F_Lc_usingInterpLc, Points, 1.e-8);
     }
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 13e8c9fa4a6202a6fddf78e9bd5bdf82267565e9..88b004f56a8396543692b7d88aaa438bc2922afa 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.114 2008-01-30 15:27:41 remacle Exp $
+// $Id: meshGFace.cpp,v 1.115 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -616,6 +616,16 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
     Msg(WARNING,"8-| Gmsh splits those edges and tries again");
     std::list<GFace *> facesToRemesh;
     remeshUnrecoveredEdges ( edgesNotRecovered, facesToRemesh);
+    std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
+    for ( ; itr != edgesNotRecovered.end() ; ++itr)
+    {
+      int p1 = itr->p1;
+      int p2 = itr->p2;
+      int tag = itr->ge->tag();
+      Msg(WARNING,"MEdge %d %d in GEdge %d",p1,p2,tag);
+    }
+
+
     delete m;
     delete [] U_;
     delete [] V_;
@@ -793,7 +803,7 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
   else if (debug){
     char name[256];
     sprintf(name,"real%d.pos",gf->tag());
-    outputScalarField(m->triangles, name,0);
+    outputScalarField(m->triangles, name,0,gf);
     sprintf(name,"param%d.pos",gf->tag());
     outputScalarField(m->triangles, name,1);
   }
@@ -1423,7 +1433,7 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
     {
       char name[245];
       sprintf(name,"surface%d-final-real.pos",gf->tag());
-      outputScalarField(m->triangles, name,0);
+      outputScalarField(m->triangles, name,0,gf);
       sprintf(name,"surface%d-final-param.pos",gf->tag());
       outputScalarField(m->triangles, name,1);
     }
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 13de94d874965ac24da35a7c97a2ef02cfe664df..6d8b475ef3baf95c924217ac3816bf60c335c4df 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceBDS.cpp,v 1.3 2008-01-30 15:27:41 remacle Exp $
+// $Id: meshGFaceBDS.cpp,v 1.4 2008-02-05 14:40:30 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -29,6 +29,7 @@
 #include "MElement.h"
 #include "Context.h"
 #include "GPoint.h"
+#include "GModel.h"
 #include "Message.h"
 #include "Numeric.h"
 #include "BDS.h"
@@ -711,8 +712,6 @@ void gmshOptimizeMeshBDS(GFace *gf,
 
 // DELAUNAY BDS
 
-
-
 void delaunayPointInsertionBDS ( GFace *gf, BDS_Mesh &m, BDS_Point *v, BDS_Face *f){
   const double p[2] = {v->u,v->v};
   
@@ -724,3 +723,51 @@ void delaunayPointInsertionBDS ( GFace *gf, BDS_Mesh &m, BDS_Point *v, BDS_Face
   int nb_swap = 0;
   gmshDelaunayizeBDS ( gf, m, nb_swap );
 }
+
+// build the BDS from a list of GFace
+// This is a TRUE copy
+BDS_Mesh * gmsh2BDS ( std::list<GFace*> & l){
+  BDS_Mesh *m = new BDS_Mesh;
+  for (std::list<GFace*>::iterator it = l.begin();it!=l.end();++it){
+    GFace *gf = *it;
+    m->add_geom (gf->tag(), 2);
+    BDS_GeomEntity *g2 = m->get_geom(gf->tag(), 2);
+    for (int i=0;i<gf->triangles.size();i++){
+      MTriangle *e = gf->triangles[i];
+      BDS_Point *p[3];
+      for (int j=0;j<3;j++){
+	p[j] = m->find_point(e->getVertex(j)->getNum());
+	if (!p[j]) {
+	  p[j] = m->add_point(e->getVertex(j)->getNum(),e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
+	  double u0,v0;
+	  parametricCoordinates ( e->getVertex(j), gf, u0, v0);
+	  p[j]->u = u0;
+	  p[j]->v = v0;
+	  m->add_geom (e->getVertex(j)->onWhat()->tag(), e->getVertex(j)->onWhat()->dim());
+	  BDS_GeomEntity *g = m->get_geom(e->getVertex(j)->onWhat()->tag(), e->getVertex(j)->onWhat()->dim());
+	  p[j]->g = g;
+	}
+      }
+      BDS_Face *f = m->add_triangle ( p[0]->iD,p[1]->iD,p[2]->iD);
+      f->g = g2;
+    }
+  }
+  return m;
+}
+
+void gmshCollapseSmallEdges (GModel &gm){
+  return;
+  gm.renumberMeshVertices(true);
+  std::list<GFace*> faces;
+  for (GModel::fiter fit = gm.firstFace(); fit != gm.lastFace(); fit++){
+    faces.push_back(*fit);
+  }
+  BDS_Mesh *pm = gmsh2BDS (faces);
+  outputScalarField(pm->triangles,"all.pos",0);
+
+  
+  for (GModel::eiter eit = gm.firstEdge(); eit != gm.lastEdge(); eit++){
+  }
+  
+
+}
diff --git a/Mesh/meshGFaceBDS.h b/Mesh/meshGFaceBDS.h
index ef11e51f9b01b9d9161e1297ae74dfd64876d77f..23aeae62e1c2f2c321f6759f9d71288fee2e8d9f 100644
--- a/Mesh/meshGFaceBDS.h
+++ b/Mesh/meshGFaceBDS.h
@@ -22,6 +22,7 @@
 
 #include <map>
 class GFace;
+class GModel;
 class BDS_Mesh;
 class BDS_Point;
 class MVertex;
@@ -38,4 +39,7 @@ void gmshOptimizeMeshBDS(GFace *gf,
 			 std::map<BDS_Point*,MVertex*> *recover_map = 0);
 
 void gmshDelaunayizeBDS ( GFace *gf, BDS_Mesh &m, int &nb_swap );
+
+void gmshCollapseSmallEdges (GModel &gm);
+
 #endif