diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index e9cbe03a0a0eba93d0cdcc73a777d8eb2734ebfb..b77c086457bf29b3c23ff81cb51bd597527cf921 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -838,12 +838,7 @@ StringXNumber GeometryOptions_Number[] = {
 StringXNumber MeshOptions_Number[] = {
   { F|O, "Algorithm" , opt_mesh_algo2d , ALGO_2D_MESHADAPT ,
     "2D mesh algorithm (1=meshadapt, 2=delaunay)" }, 
-  { F|O, "Algorithm3D" , opt_mesh_algo3d ,
-#if defined(HAVE_TETGEN)
-    ALGO_3D_DELAUNAY,
-#else
-    ALGO_3D_NETGEN,
-#endif
+  { F|O, "Algorithm3D" , opt_mesh_algo3d , ALGO_3D_NETGEN ,
     "3D mesh algorithm (1=delaunay, 4=netgen)" }, 
   { F|O, "AngleSmoothNormals" , opt_mesh_angle_smooth_normals , 30.0 ,
     "Threshold angle below which normals are not smoothed" }, 
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 4b26e78911d9107caf7a82e0e10edb3da8818d75..ee7d04166d6a23284ec07c55c8fa33f720a7e748 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.502 2007-01-12 19:47:52 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.503 2007-01-16 11:31:40 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -3758,8 +3758,6 @@ void mesh_remesh_cb(CALLBACK_ARGS)
 
 void mesh_inspect_cb(CALLBACK_ARGS)
 {
-  char *str = (char*)data;
-
   std::vector<GVertex*> vertices;
   std::vector<GEdge*> edges;
   std::vector<GFace*> faces;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 6368d14030362ee508f4dd6d775f3e9ab88b3119..1e9042ee9b6d49d9811bcf56ec7550d7caf952d0 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -51,8 +51,8 @@ class GModel
   std::map<int, std::string> physicalNames;
 
  public:
-  GModel() : modelName("Untitled"), normals(0), meshSize(-1.) {}
-  GModel(const std::string &name) : modelName(name), normals(0), meshSize(-1.) {}
+  GModel() : meshSize(-1.), modelName("Untitled"), normals(0) {}
+  GModel(const std::string &name) : meshSize(-1.), modelName(name), normals(0) {}
   virtual ~GModel(){ deleteOCCInternals(); destroy(); }
   
   typedef std::set<GRegion*, GEntityLessThan>::iterator riter;
@@ -149,7 +149,7 @@ class GModel
   // Get or set the global mesh size for the model (default value is a
   // tenth of the size of the domain)
   double getMeshSize();
-  void setMeshSize(const double s) {meshSize = s;}
+  void setMeshSize(const double s) { meshSize = s; }
 
   // A container for smooth normals
   smooth_normals *normals;
diff --git a/Geo/GModelIO_Fourier.cpp b/Geo/GModelIO_Fourier.cpp
index 6641b1eb584d306110c44024c2d47fbfaf7992ff..5aa3bccd4771fa340becf4610b8517ea1e9390fb 100644
--- a/Geo/GModelIO_Fourier.cpp
+++ b/Geo/GModelIO_Fourier.cpp
@@ -170,7 +170,7 @@ public:
       else{
 	double allPou = 0.;
 	for(unsigned int i = 0; i < param.size(); i++){
-	  double u2 = param[i].second[0], v2 = param[i].second[1];
+	  //double u2 = param[i].second[0], v2 = param[i].second[1];
 	  double pou2=1;
 	  //FM->GetPou(param[i].first->tag(), u2, v2, pou2);
 	  allPou += pou2;
@@ -637,7 +637,7 @@ void fourierFace::meshBoundary()
   std::vector<double> u, v;
   FM->GetBoundary_Points(tag(), u, v);
 
-  if(2*nu+2*nv != u.size()){
+  if(2*nu+2*nv != (int)u.size()){
     Msg(INFO, "Special planar patch from YoungAe: %d", tag());
 #if 1 // transfinite, by hand -- WARNING
     _plane = 1; // to enable smoothing of transfinite meshes
@@ -673,13 +673,13 @@ void fourierFace::meshBoundary()
     _e[2]->addFace(this);
     _e[3] = new fourierEdge(model(), 4, _v[3], _v[0]);
     _e[3]->addFace(this);
-    for(unsigned int i = corners[0] + 1; i < corners[1]; i++)
+    for(int i = corners[0] + 1; i < corners[1]; i++)
       _e[0]->mesh_vertices.push_back(verts[i]);
-    for(unsigned int i = corners[1] + 1; i < corners[2]; i++)
+    for(int i = corners[1] + 1; i < corners[2]; i++)
       _e[1]->mesh_vertices.push_back(verts[i]);
-    for(unsigned int i = corners[2] + 1; i < corners[3]; i++)
+    for(int i = corners[2] + 1; i < corners[3]; i++)
       _e[2]->mesh_vertices.push_back(verts[i]);
-    for(unsigned int i = corners[3] + 1; i < verts.size(); i++)
+    for(int i = corners[3] + 1; i < (int)verts.size(); i++)
       _e[3]->mesh_vertices.push_back(verts[i]);
     l_edges.push_back(_e[0]); l_dirs.push_back(1);
     l_edges.push_back(_e[1]); l_dirs.push_back(1);
@@ -733,19 +733,19 @@ void fourierFace::meshBoundary()
   _e[3] = new fourierEdge(model(), 4, _v[3], _v[0]);
   _e[3]->addFace(this);
 
-  for(unsigned int i = corners[0] + 1; i < corners[1] - 1; i++){
+  for(int i = corners[0] + 1; i < corners[1] - 1; i++){
     GPoint p = point(u[i], v[i]);
     _e[0]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
   }
-  for(unsigned int i = corners[1] + 1; i < corners[2] - 1; i++){
+  for(int i = corners[1] + 1; i < corners[2] - 1; i++){
     GPoint p = point(u[i], v[i]);
     _e[1]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
   }
-  for(unsigned int i = corners[2] + 1; i < corners[3] - 1; i++){
+  for(int i = corners[2] + 1; i < corners[3] - 1; i++){
     GPoint p = point(u[i], v[i]);
     _e[2]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
   }
-  for(unsigned int i = corners[3] + 1; i < u.size() - 1; i++){
+  for(int i = corners[3] + 1; i < (int)u.size() - 1; i++){
     GPoint p = point(u[i], v[i]);
     _e[3]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
   }
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index bf2ac3c8fa943355a3ffb1eb687b16594357265d..e4cbd77feca9b390cc87aeb623a083cd5d6d4cfd 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_OCC.cpp,v 1.14 2006-11-29 16:57:00 remacle Exp $
+// $Id: GModelIO_OCC.cpp,v 1.15 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -40,13 +40,12 @@
 #include "ShapeFix.hxx"
 #include "ShapeFix_FixSmallFace.hxx"
 
-class OCC_Internals
-{
+class OCC_Internals {
   TopoDS_Shape shape;
   TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;
   double tolerance;
 public:
-  OCC_Internals ()
+  OCC_Internals()
   {
     somap.Clear();
     shmap.Clear();
@@ -56,15 +55,16 @@ public:
     vmap.Clear();
     tolerance = 1.e-3;
   }
-  void HealGeometry (bool fixsmalledges = true, bool fixspotstripfaces = true, bool sewfaces = false, bool makesolids = false);
+  void HealGeometry(bool fixsmalledges = true, bool fixspotstripfaces = true,
+		    bool sewfaces = false, bool makesolids = false);
   void loadSTEP(const char *);
   void loadIGES(const char *);
   void loadBREP(const char *);  
-  void buildGModel (GModel *gm);
-  void buildLists ();
+  void buildGModel(GModel *gm);
+  void buildLists();
 };
 
-void OCC_Internals :: buildLists ()
+void OCC_Internals::buildLists()
 {
   TopExp_Explorer exp0, exp1, exp2, exp3, exp4, exp5;
   somap.Clear();
@@ -74,517 +74,419 @@ void OCC_Internals :: buildLists ()
   emap.Clear();
   vmap.Clear();
   
-  for (exp0.Init(shape, TopAbs_SOLID);
-       exp0.More(); exp0.Next())
-    {
-      TopoDS_Solid solid = TopoDS::Solid (exp0.Current());
-      
-      if (somap.FindIndex(TopoDS::Solid (exp0.Current())) < 1)
-	{
-	  somap.Add (TopoDS::Solid (exp0.Current()));
-	  
-	  for (exp1.Init(exp0.Current(), TopAbs_SHELL);
-	       exp1.More(); exp1.Next())
-	    {
-	      TopoDS_Shell shell = TopoDS::Shell (exp1.Current().Composed (exp0.Current().Orientation()));
-	      if (shmap.FindIndex(shell) < 1)
-		{
-		  shmap.Add (shell);
-		  
-		  for (exp2.Init(shell, TopAbs_FACE);
-		       exp2.More(); exp2.Next())
-		    {
-		      TopoDS_Face face = TopoDS::Face(exp2.Current().Composed(shell.Orientation()));
-		      if (fmap.FindIndex(face) < 1)
-			{
-			  fmap.Add (face);
-			  
-			  for (exp3.Init(exp2.Current(), TopAbs_WIRE);
-			       exp3.More(); exp3.Next())
-			    {
-			      TopoDS_Wire wire = TopoDS::Wire (exp3.Current().Composed(face.Orientation()));
-			      if (wmap.FindIndex(wire) < 1)
-				{
-				  wmap.Add (wire);
-				  
-				  for (exp4.Init(exp3.Current(), TopAbs_EDGE);
-				       exp4.More(); exp4.Next())
-				    {
-				      TopoDS_Edge edge = TopoDS::Edge(exp4.Current().Composed(wire.Orientation()));
-				      if (emap.FindIndex(edge) < 1)
-					{
-					  emap.Add (edge);
-					  for (exp5.Init(exp4.Current(), TopAbs_VERTEX);
-					       exp5.More(); exp5.Next())
-					    {
-					      TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
-					      if (vmap.FindIndex(vertex) < 1)
-						vmap.Add (vertex);
-					    }
-					}
-				    }
-				}
-			    }
-			}
+  for(exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()){
+    TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
+    if(somap.FindIndex(TopoDS::Solid(exp0.Current())) < 1){
+      somap.Add(TopoDS::Solid(exp0.Current()));
+
+      for(exp1.Init(exp0.Current(), TopAbs_SHELL); exp1.More(); exp1.Next()){
+	TopoDS_Shell shell = TopoDS::Shell(exp1.Current().Composed(exp0.Current().Orientation()));
+	if(shmap.FindIndex(shell) < 1){
+	  shmap.Add(shell);
+
+	  for(exp2.Init(shell, TopAbs_FACE); exp2.More(); exp2.Next()){
+	    TopoDS_Face face = TopoDS::Face(exp2.Current().Composed(shell.Orientation()));
+	    if(fmap.FindIndex(face) < 1){
+	      fmap.Add(face);
+
+	      for(exp3.Init(exp2.Current(), TopAbs_WIRE); exp3.More(); exp3.Next()){
+		TopoDS_Wire wire = TopoDS::Wire(exp3.Current().Composed(face.Orientation()));
+		if(wmap.FindIndex(wire) < 1){
+		  wmap.Add(wire);
+
+		  for(exp4.Init(exp3.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()){
+		    TopoDS_Edge edge = TopoDS::Edge(exp4.Current().Composed(wire.Orientation()));
+		    if(emap.FindIndex(edge) < 1){
+		      emap.Add(edge);
+
+		      for(exp5.Init(exp4.Current(), TopAbs_VERTEX); exp5.More(); exp5.Next()){
+			TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
+			if(vmap.FindIndex(vertex) < 1)
+			  vmap.Add(vertex);
+		      }
 		    }
+		  }
 		}
+	      }
 	    }
+	  }
 	}
+      }
     }
+  }
   
   // Free Shells
-  for (exp1.Init(exp0.Current(), TopAbs_SHELL, TopAbs_SOLID);
-       exp1.More(); exp1.Next())
-    {
-      TopoDS_Shape shell = exp1.Current().Composed (exp0.Current().Orientation());
-      if (shmap.FindIndex(shell) < 1)
-	{
-	  shmap.Add (shell);
-	  
-	  for (exp2.Init(shell, TopAbs_FACE);
-	       exp2.More(); exp2.Next())
-	    {
-	      TopoDS_Face face = TopoDS::Face(exp2.Current().Composed(shell.Orientation()));
-	      if (fmap.FindIndex(face) < 1)
-		{
-		  fmap.Add (face);
+  for(exp1.Init(exp0.Current(), TopAbs_SHELL, TopAbs_SOLID); exp1.More(); exp1.Next()){
+    TopoDS_Shape shell = exp1.Current().Composed(exp0.Current().Orientation());
+    if(shmap.FindIndex(shell) < 1){
+      shmap.Add(shell);
+
+      for(exp2.Init(shell, TopAbs_FACE); exp2.More(); exp2.Next()){
+	TopoDS_Face face = TopoDS::Face(exp2.Current().Composed(shell.Orientation()));
+	if(fmap.FindIndex(face) < 1){
+	  fmap.Add(face);
 		  
-		  for (exp3.Init(exp2.Current(), TopAbs_WIRE);
-		       exp3.More(); exp3.Next())
-		    {
-		      TopoDS_Wire wire = TopoDS::Wire (exp3.Current());
-		      if (wmap.FindIndex(wire) < 1)
-			{
-			  wmap.Add (wire);
-			  
-			  for (exp4.Init(exp3.Current(), TopAbs_EDGE);
-			       exp4.More(); exp4.Next())
-			    {
-			      TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
-			      if (emap.FindIndex(edge) < 1)
-				{
-				  emap.Add (edge);
-				  for (exp5.Init(exp4.Current(), TopAbs_VERTEX);
-				       exp5.More(); exp5.Next())
-				    {
-				      TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
-				      if (vmap.FindIndex(vertex) < 1)
-					vmap.Add (vertex);
-				    }
-				}
-			    }
-			}
-		    }
+	  for(exp3.Init(exp2.Current(), TopAbs_WIRE); exp3.More(); exp3.Next()){
+	    TopoDS_Wire wire = TopoDS::Wire(exp3.Current());
+	    if(wmap.FindIndex(wire) < 1){
+	      wmap.Add(wire);
+
+	      for(exp4.Init(exp3.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()){
+		TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
+		if(emap.FindIndex(edge) < 1){
+		  emap.Add(edge);
+
+		  for(exp5.Init(exp4.Current(), TopAbs_VERTEX); exp5.More(); exp5.Next()){
+		    TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
+		    if(vmap.FindIndex(vertex) < 1)
+		      vmap.Add(vertex);
+		  }
 		}
+	      }
 	    }
+	  }
 	}
+      }
     }
-  
-  
+  }
+    
   // Free Faces
-  
-  for (exp2.Init(shape, TopAbs_FACE, TopAbs_SHELL);
-       exp2.More(); exp2.Next())
-    {
-      TopoDS_Face face = TopoDS::Face(exp2.Current());
-      if (fmap.FindIndex(face) < 1)
-	{
-	  fmap.Add (face);
+  for(exp2.Init(shape, TopAbs_FACE, TopAbs_SHELL); exp2.More(); exp2.Next()){
+    TopoDS_Face face = TopoDS::Face(exp2.Current());
+    if(fmap.FindIndex(face) < 1){
+      fmap.Add(face);
 	  
-	  for (exp3.Init(exp2.Current(), TopAbs_WIRE);
-	       exp3.More(); exp3.Next())
-	    {
-	      TopoDS_Wire wire = TopoDS::Wire (exp3.Current());
-	      if (wmap.FindIndex(wire) < 1)
-		{
-		  wmap.Add (wire);
-		  
-		  for (exp4.Init(exp3.Current(), TopAbs_EDGE);
-		       exp4.More(); exp4.Next())
-		    {
-		      TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
-		      if (emap.FindIndex(edge) < 1)
-			{
-			  emap.Add (edge);
-			  for (exp5.Init(exp4.Current(), TopAbs_VERTEX);
-			       exp5.More(); exp5.Next())
-			    {
-			      TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
-			      if (vmap.FindIndex(vertex) < 1)
-				vmap.Add (vertex);
-			    }
-			}
-		    }
-		}
+      for(exp3.Init(exp2.Current(), TopAbs_WIRE); exp3.More(); exp3.Next()){
+	TopoDS_Wire wire = TopoDS::Wire(exp3.Current());
+	if(wmap.FindIndex(wire) < 1){
+	  wmap.Add(wire);
+	  
+	  for(exp4.Init(exp3.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()){
+	    TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
+	    if(emap.FindIndex(edge) < 1){
+	      emap.Add(edge);
+
+	      for(exp5.Init(exp4.Current(), TopAbs_VERTEX); exp5.More(); exp5.Next()){
+		TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
+		if(vmap.FindIndex(vertex) < 1)
+		  vmap.Add(vertex);
+	      }
 	    }
+	  }
 	}
+      }
     }
-
+  }
 
   // Free Wires
-  
-  for (exp3.Init(shape, TopAbs_WIRE, TopAbs_FACE);
-       exp3.More(); exp3.Next())
-    {
-      TopoDS_Wire wire = TopoDS::Wire (exp3.Current());
-      if (wmap.FindIndex(wire) < 1)
-	{
-	  wmap.Add (wire);
-	  
-	  for (exp4.Init(exp3.Current(), TopAbs_EDGE);
-	       exp4.More(); exp4.Next())
-	    {
-	      TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
-	      if (emap.FindIndex(edge) < 1)
-		{
-		  emap.Add (edge);
-		  for (exp5.Init(exp4.Current(), TopAbs_VERTEX);
-		       exp5.More(); exp5.Next())
-		    {
-		      TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
-		      if (vmap.FindIndex(vertex) < 1)
-			vmap.Add (vertex);
-		    }
-		}
-	    }
+  for(exp3.Init(shape, TopAbs_WIRE, TopAbs_FACE); exp3.More(); exp3.Next()){
+    TopoDS_Wire wire = TopoDS::Wire(exp3.Current());
+    if(wmap.FindIndex(wire) < 1){
+      wmap.Add(wire);
+      
+      for(exp4.Init(exp3.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()){
+	TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
+	if(emap.FindIndex(edge) < 1){
+	  emap.Add(edge);
+
+	  for(exp5.Init(exp4.Current(), TopAbs_VERTEX); exp5.More(); exp5.Next()){
+	    TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
+	    if(vmap.FindIndex(vertex) < 1)
+	      vmap.Add(vertex);
+	  }
 	}
+      }
     }
-
+  }
 
   // Free Edges
-  
-  for (exp4.Init(shape, TopAbs_EDGE, TopAbs_WIRE);
-       exp4.More(); exp4.Next())
-    {
-      TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
-      if (emap.FindIndex(edge) < 1)
-	{
-	  emap.Add (edge);
-	  for (exp5.Init(exp4.Current(), TopAbs_VERTEX);
-	       exp5.More(); exp5.Next())
-	    {
-	      TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
-	      if (vmap.FindIndex(vertex) < 1)
-		vmap.Add (vertex);
-	    }
-	}
+  for(exp4.Init(shape, TopAbs_EDGE, TopAbs_WIRE); exp4.More(); exp4.Next()){
+    TopoDS_Edge edge = TopoDS::Edge(exp4.Current());
+    if(emap.FindIndex(edge) < 1){
+      emap.Add(edge);
+
+      for(exp5.Init(exp4.Current(), TopAbs_VERTEX); exp5.More(); exp5.Next()){
+	TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
+	if(vmap.FindIndex(vertex) < 1)
+	  vmap.Add(vertex);
+      }
     }
-
+  }
 
   // Free Vertices
+  for(exp5.Init(shape, TopAbs_VERTEX, TopAbs_EDGE); exp5.More(); exp5.Next()){
+    TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
+    if(vmap.FindIndex(vertex) < 1)
+      vmap.Add(vertex);
+  }    
   
-  for (exp5.Init(shape, TopAbs_VERTEX, TopAbs_EDGE);
-       exp5.More(); exp5.Next())
-    {
-      TopoDS_Vertex vertex = TopoDS::Vertex(exp5.Current());
-      if (vmap.FindIndex(vertex) < 1)
-	vmap.Add (vertex);
-    }    
-
 }
 
-
-void OCC_Internals ::   HealGeometry (bool fixsmalledges , bool fixspotstripfaces, bool sewfaces, bool makesolids)
+void OCC_Internals::HealGeometry(bool fixsmalledges, bool fixspotstripfaces, 
+				 bool sewfaces, bool makesolids)
 {
   int nrc = 0, nrcs = 0;
   TopExp_Explorer e;
-  for (e.Init(shape, TopAbs_COMPOUND); e.More(); e.Next()) nrc++;
-  for (e.Init(shape, TopAbs_COMPSOLID); e.More(); e.Next()) nrcs++;
+  for(e.Init(shape, TopAbs_COMPOUND); e.More(); e.Next()) nrc++;
+  for(e.Init(shape, TopAbs_COMPSOLID); e.More(); e.Next()) nrcs++;
 
   double surfacecont = 0;
 
-  for (int i = 1; i <= fmap.Extent(); i++)
-    {
-      GProp_GProps system;
-      BRepGProp::LinearProperties(fmap(i), system);
-      surfacecont += system.Mass();
-    }
-
-  cout << "Starting geometry healing procedure (tolerance: " << tolerance << ")" << endl
-       << "-----------------------------------" << endl;
-
-  if (fixsmalledges)
-    {
-      cout << endl << "- fixing small edges" << endl;
+  for(int i = 1; i <= fmap.Extent(); i++){
+    GProp_GProps system;
+    BRepGProp::LinearProperties(fmap(i), system);
+    surfacecont += system.Mass();
+  }
 
-      Handle(ShapeFix_Wire) sfw;
-      Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape;
-      rebuild->Apply(shape);
+  Msg(INFO, "Healing geometry (tolerance=%g)", tolerance);
 
-      for (int i = 1; i <= fmap.Extent(); i++)
-	{
-	  TopExp_Explorer exp1;
-	  for (exp1.Init (fmap(i), TopAbs_WIRE); exp1.More(); exp1.Next())
-	    {
-	      TopoDS_Wire oldwire = TopoDS::Wire(exp1.Current());
-	      sfw = new ShapeFix_Wire (oldwire, TopoDS::Face(fmap(i)),tolerance);
-	      sfw->ModifyTopologyMode() = Standard_True;
-
-	      if (sfw->FixSmall (false, tolerance))
-		{
-		  cout << "Fixed small edge in wire " << wmap.FindIndex (oldwire) << endl;
-		  TopoDS_Wire newwire = sfw->Wire();
-		  rebuild->Replace(oldwire, newwire, Standard_False);
-		}
-	      if ((sfw->StatusSmall(ShapeExtend_FAIL1)) ||
-		  (sfw->StatusSmall(ShapeExtend_FAIL2)) ||
-		  (sfw->StatusSmall(ShapeExtend_FAIL3)))
-		cout << "Failed to fix small edge in wire " << wmap.FindIndex (oldwire) << endl;
+  if(fixsmalledges){
+    Msg(INFO, "- fixing small edges");
 
-	      
-	    }
+    Handle(ShapeFix_Wire) sfw;
+    Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape;
+    rebuild->Apply(shape);
+    
+    for(int i = 1; i <= fmap.Extent(); i++){
+      TopExp_Explorer exp1;
+      for(exp1.Init(fmap(i), TopAbs_WIRE); exp1.More(); exp1.Next()){
+	TopoDS_Wire oldwire = TopoDS::Wire(exp1.Current());
+	sfw = new ShapeFix_Wire(oldwire, TopoDS::Face(fmap(i)),tolerance);
+	sfw->ModifyTopologyMode() = Standard_True;
+	
+	if(sfw->FixSmall(false, tolerance)){
+	  Msg(INFO, "Fixed small edge in wire %d", wmap.FindIndex(oldwire));
+	  TopoDS_Wire newwire = sfw->Wire();
+	  rebuild->Replace(oldwire, newwire, Standard_False);
 	}
-
-      shape = rebuild->Apply(shape);
-
-
-
-      {
+	if((sfw->StatusSmall(ShapeExtend_FAIL1)) ||
+	   (sfw->StatusSmall(ShapeExtend_FAIL2)) ||
+	   (sfw->StatusSmall(ShapeExtend_FAIL3)))
+	  Msg(INFO, "Failed to fix small edge in wire %d",  wmap.FindIndex(oldwire));
+      }
+    }
+    shape = rebuild->Apply(shape);
+    
+    {
       Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape;
       rebuild->Apply(shape);
       TopExp_Explorer exp1;
-      for (exp1.Init (shape, TopAbs_EDGE); exp1.More(); exp1.Next())
-	{
-	  TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
-	  if (vmap.FindIndex(TopExp::FirstVertex (edge)) == 
-	      vmap.FindIndex(TopExp::LastVertex (edge)))
-	    {
-	      GProp_GProps system;
-	      BRepGProp::LinearProperties(edge, system);
-	      if (system.Mass() < tolerance)
-		{
-		  cout << "removing degenerated edge " << emap.FindIndex(edge) << endl;
-		  rebuild->Remove(edge, false);
-		}
-	    }
+      for(exp1.Init(shape, TopAbs_EDGE); exp1.More(); exp1.Next()){
+	TopoDS_Edge edge = TopoDS::Edge(exp1.Current());
+	if(vmap.FindIndex(TopExp::FirstVertex(edge)) == 
+	   vmap.FindIndex(TopExp::LastVertex(edge))){
+	  GProp_GProps system;
+	  BRepGProp::LinearProperties(edge, system);
+	  if(system.Mass() < tolerance){
+	    Msg(INFO, "removing degenerated edge %d", emap.FindIndex(edge));
+	    rebuild->Remove(edge, false);
+	  }
 	}
-      shape = rebuild->Apply(shape);
       }
-
-
-      Handle(ShapeFix_Wireframe) sfwf = new ShapeFix_Wireframe;
-      sfwf->SetPrecision(tolerance);
-      sfwf->Load (shape);
-
-      if (sfwf->FixSmallEdges())
-	{
-	  cout << endl << "- fixing wire frames" << endl;  
-	  if (sfwf->StatusSmallEdges(ShapeExtend_OK)) cout << "no small edges found" << endl;
-	  if (sfwf->StatusSmallEdges(ShapeExtend_DONE1)) cout << "some small edges fixed" << endl;
-	  if (sfwf->StatusSmallEdges(ShapeExtend_FAIL1)) cout << "failed to fix some small edges" << endl;
-	}
+      shape = rebuild->Apply(shape);
+    }
+    
+    Handle(ShapeFix_Wireframe) sfwf = new ShapeFix_Wireframe;
+    sfwf->SetPrecision(tolerance);
+    sfwf->Load(shape);
+    
+    if(sfwf->FixSmallEdges()){
+      Msg(INFO, "- fixing wire frames");
+      if(sfwf->StatusSmallEdges(ShapeExtend_OK)) Msg(INFO, "no small edges found");
+      if(sfwf->StatusSmallEdges(ShapeExtend_DONE1)) Msg(INFO, "some small edges fixed");
+      if(sfwf->StatusSmallEdges(ShapeExtend_FAIL1)) Msg(INFO, "failed to fix some small edges");
+    }
   
-
-      if (sfwf->FixWireGaps())
-	{
-	  cout << endl << "- fixing wire gaps" << endl;
-	  if (sfwf->StatusWireGaps(ShapeExtend_OK)) cout << "no gaps found" << endl;
-	  if (sfwf->StatusWireGaps(ShapeExtend_DONE1)) cout << "some 2D gaps fixed" << endl;
-	  if (sfwf->StatusWireGaps(ShapeExtend_DONE2)) cout << "some 3D gaps fixed" << endl;
-	  if (sfwf->StatusWireGaps(ShapeExtend_FAIL1)) cout << "failed to fix some 2D gaps" << endl;
-	  if (sfwf->StatusWireGaps(ShapeExtend_FAIL2)) cout << "failed to fix some 3D gaps" << endl;
-	}
-      
-
-      shape = sfwf->Shape();
+    if(sfwf->FixWireGaps()){
+      Msg(INFO, "- fixing wire gaps");
+      if(sfwf->StatusWireGaps(ShapeExtend_OK)) Msg(INFO, "no gaps found");
+      if(sfwf->StatusWireGaps(ShapeExtend_DONE1)) Msg(INFO, "some 2D gaps fixed");
+      if(sfwf->StatusWireGaps(ShapeExtend_DONE2)) Msg(INFO, "some 3D gaps fixed");
+      if(sfwf->StatusWireGaps(ShapeExtend_FAIL1)) Msg(INFO, "failed to fix some 2D gaps");
+      if(sfwf->StatusWireGaps(ShapeExtend_FAIL2)) Msg(INFO, "failed to fix some 3D gaps");
     }
-
-
-
-
-
-  if (fixspotstripfaces)
-    {
+    
+    shape = sfwf->Shape();
+  }
   
-      cout << endl << "- fixing spot and strip faces" << endl;
-      Handle(ShapeFix_FixSmallFace) sffsm = new ShapeFix_FixSmallFace();
-      sffsm -> Init (shape);
-      sffsm -> SetPrecision (tolerance);
-      sffsm -> Perform();
-      
-      shape = sffsm -> FixShape();
+  if(fixspotstripfaces){
+    Msg(INFO, "- fixing spot and strip faces");
+    Handle(ShapeFix_FixSmallFace) sffsm = new ShapeFix_FixSmallFace();
+    sffsm->Init(shape);
+    sffsm->SetPrecision(tolerance);
+    sffsm->Perform();
+    
+    shape = sffsm->FixShape();
+  }
+  
+  if(sewfaces){
+    Msg(INFO, "- sewing faces");
+
+    TopExp_Explorer exp0;
+    
+    BRepOffsetAPI_Sewing sewedObj(tolerance);
+    
+    for(exp0.Init(shape, TopAbs_FACE); exp0.More(); exp0.Next()){
+      TopoDS_Face face = TopoDS::Face(exp0.Current());
+      sewedObj.Add(face);
     }
-
-  if (sewfaces)
-    {
-      cout << endl << "- sewing faces" << endl;
-
-      TopExp_Explorer exp0;
-
-      BRepOffsetAPI_Sewing sewedObj(tolerance);
-
-      for (exp0.Init (shape, TopAbs_FACE); exp0.More(); exp0.Next())
-	{
-	  TopoDS_Face face = TopoDS::Face (exp0.Current());
-	  sewedObj.Add (face);
-	}
-      
-      sewedObj.Perform();
+    
+    sewedObj.Perform();
+    
+    if(!sewedObj.SewedShape().IsNull())
+      shape = sewedObj.SewedShape();
+    else
+      Msg(INFO, " not possible");
+  }
   
-      if (!sewedObj.SewedShape().IsNull())
-	shape = sewedObj.SewedShape();
-      else
-	cout << " not possible";
+  if(makesolids){  
+    Msg(INFO, "- making solids");
+    
+    TopExp_Explorer exp0;
+    
+    BRepBuilderAPI_MakeSolid ms;
+    int count = 0;
+    for(exp0.Init(shape, TopAbs_SHELL); exp0.More(); exp0.Next()){
+      count++;
+      ms.Add(TopoDS::Shell(exp0.Current()));
     }
-
-  if (makesolids)
-    {  
-      cout << endl << "- making solids" << endl;
-      
-      TopExp_Explorer exp0;
-
-      BRepBuilderAPI_MakeSolid ms;
-      int count = 0;
-      for (exp0.Init(shape, TopAbs_SHELL); exp0.More(); exp0.Next())
-	{
-	  count++;
-	  ms.Add (TopoDS::Shell(exp0.Current()));
-	}
-      
-      if (!count)
-	{
-	  cout << " not possible (no shells)" << endl;
+    
+    if(!count){
+      Msg(INFO, " not possible (no shells)");
+    }
+    else{
+      BRepCheck_Analyzer ba(ms);
+      if(ba.IsValid()){
+	Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
+	sfs->Init(ms);
+	sfs->SetPrecision(tolerance);
+	sfs->SetMaxTolerance(tolerance);
+	sfs->Perform();
+	shape = sfs->Shape();
+	
+	for(exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()){
+	  TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
+	  TopoDS_Solid newsolid = solid;
+	  BRepLib::OrientClosedSolid(newsolid);
+	  Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape;
+	  // rebuild->Apply(shape);
+	  rebuild->Replace(solid, newsolid, Standard_False);
+	  TopoDS_Shape newshape = rebuild->Apply(shape, TopAbs_COMPSOLID, 1);
+	  // TopoDS_Shape newshape = rebuild->Apply(shape);
+	  shape = newshape;
 	}
+      }
       else
-	{
-	  BRepCheck_Analyzer ba(ms);
-	  if (ba.IsValid ())
-	    {
-	      Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
-	      sfs->Init (ms);
-	      sfs->SetPrecision(tolerance);
-	      sfs->SetMaxTolerance(tolerance);
-	      sfs->Perform();
-	      shape = sfs->Shape();
-	      
-	      for (exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next())
-		{
-		  TopoDS_Solid solid = TopoDS::Solid(exp0.Current());
-		  TopoDS_Solid newsolid = solid;
-		  BRepLib::OrientClosedSolid (newsolid);
-		  Handle_ShapeBuild_ReShape rebuild = new ShapeBuild_ReShape;
-		  //		  rebuild->Apply(shape);
-		  rebuild->Replace(solid, newsolid, Standard_False);
-		  TopoDS_Shape newshape = rebuild->Apply(shape, TopAbs_COMPSOLID, 1);
-		  //		  TopoDS_Shape newshape = rebuild->Apply(shape);
-		  shape = newshape;
-		}
-	    }
-	  else
-	    cout << " not possible" << endl;
-	}
+	Msg(INFO, " not possible");
     }
+  }
   buildLists();
 }
 
-
-void OCC_Internals :: loadBREP (const char *fn)
+void OCC_Internals::loadBREP(const char *fn)
 {
   BRep_Builder aBuilder;
-  Standard_Boolean result = BRepTools::Read( shape, (char*)fn, aBuilder );
-  BRepTools::Clean (shape);
+  BRepTools::Read(shape, (char*)fn, aBuilder);
+  BRepTools::Clean(shape);
   buildLists();
   HealGeometry();
-  BRepTools::Clean (shape);
+  BRepTools::Clean(shape);
 }
 
-void OCC_Internals :: loadSTEP (const char *fn)
+void OCC_Internals::loadSTEP(const char *fn)
 {
   STEPControl_Reader reader;
-  Standard_Integer stat = reader.ReadFile((char*)fn);
-  Standard_Integer nb = reader.NbRootsForTransfer();
-  reader.TransferRoots (); 
+  reader.ReadFile((char*)fn);
+  reader.NbRootsForTransfer();
+  reader.TransferRoots(); 
   shape = reader.OneShape();  
-  BRepTools::Clean (shape);
+  BRepTools::Clean(shape);
   buildLists();
   HealGeometry();
-  BRepTools::Clean (shape);
+  BRepTools::Clean(shape);
 }
 
-
-void OCC_Internals :: loadIGES (const char *fn)
+void OCC_Internals::loadIGES(const char *fn)
 {
   IGESControl_Reader reader;
-  Standard_Integer stat = reader.ReadFile((char*)fn);
-  Standard_Integer nb = reader.NbRootsForTransfer();
-  reader.TransferRoots (); 
+  reader.ReadFile((char*)fn);
+  reader.NbRootsForTransfer();
+  reader.TransferRoots(); 
   shape = reader.OneShape();  
-   BRepTools::Clean (shape);
-   buildLists();
-   HealGeometry();
-   BRepTools::Clean (shape);
-   buildLists();
-   HealGeometry();
-  BRepTools::Clean (shape);
+  BRepTools::Clean(shape);
+  buildLists();
+  HealGeometry();
+  BRepTools::Clean(shape);
+  buildLists();
+  HealGeometry();
+  BRepTools::Clean(shape);
 }
 
-void OCC_Internals :: buildGModel (GModel *model)
+void OCC_Internals::buildGModel(GModel *model)
 {
   // building geom vertices
   int nvertices = vmap.Extent();
-  for (int i = 1; i <= nvertices; i++)
-    {
-      OCCVertex *v = new OCCVertex (model, i, TopoDS::Vertex(vmap(i)));
-      model->add(v);
-    }
+  for(int i = 1; i <= nvertices; i++){
+    OCCVertex *v = new OCCVertex(model, i, TopoDS::Vertex(vmap(i)));
+    model->add(v);
+  }
   // building geom edges
   int nedges = emap.Extent();
-  for (int i = 1; i <= nedges; i++)
-    {
-      TopoDS_Edge edge = TopoDS::Edge(emap(i));
-      int i1 = vmap.FindIndex(TopExp::FirstVertex (edge)); 
-      int i2 = vmap.FindIndex(TopExp::LastVertex (edge));
-      GVertex *v1 = model->vertexByTag ( i1);
-      GVertex *v2 = model->vertexByTag ( i2);
-      OCCEdge *e = new OCCEdge (model, edge, i, v1, v2);
-      model->add(e);
-    }
+  for(int i = 1; i <= nedges; i++){
+    TopoDS_Edge edge = TopoDS::Edge(emap(i));
+    int i1 = vmap.FindIndex(TopExp::FirstVertex(edge)); 
+    int i2 = vmap.FindIndex(TopExp::LastVertex(edge));
+    GVertex *v1 = model->vertexByTag(i1);
+    GVertex *v2 = model->vertexByTag(i2);
+    OCCEdge *e = new OCCEdge(model, edge, i, v1, v2);
+    model->add(e);
+  }
   // building geom faces
   int nfaces = fmap.Extent();
-  for (int i = 1; i <= nfaces; i++)
-    {
-      TopoDS_Face face = TopoDS::Face(fmap(i));
-      OCCFace *f = new OCCFace (model, face, i, emap);
-      model->add(f);
-    }
+  for(int i = 1; i <= nfaces; i++){
+    TopoDS_Face face = TopoDS::Face(fmap(i));
+    OCCFace *f = new OCCFace(model, face, i, emap);
+    model->add(f);
+  }
   // building geom regions
   int nvolumes = somap.Extent();
-  for (int i = 1; i <= nvolumes; i++)
-    {
-      TopoDS_Solid solid = TopoDS::Solid(somap(i));
-      OCCRegion *r = new OCCRegion (model, solid, i, fmap);
-      model->add(r);
-    }
-
+  for(int i = 1; i <= nvolumes; i++){
+    TopoDS_Solid solid = TopoDS::Solid(somap(i));
+    OCCRegion *r = new OCCRegion(model, solid, i, fmap);
+    model->add(r);
+  }
 }
 
 int GModel::readOCCSTEP(const std::string &fn)
 {
   occ_internals = new OCC_Internals;
-  occ_internals->loadSTEP (fn.c_str());
-  occ_internals->buildLists ();
-  occ_internals->buildGModel (this);
+  occ_internals->loadSTEP(fn.c_str());
+  occ_internals->buildLists();
+  occ_internals->buildGModel(this);
   return 1;
 }
+
 int GModel::readOCCIGES(const std::string &fn)
 {
   occ_internals = new OCC_Internals;
-  occ_internals->loadIGES (fn.c_str());
-  occ_internals->buildLists ();
-  occ_internals->buildGModel (this);
+  occ_internals->loadIGES(fn.c_str());
+  occ_internals->buildLists();
+  occ_internals->buildGModel(this);
   return 1;
 }
+
 int GModel::readOCCBREP(const std::string &fn)
 {
   occ_internals = new OCC_Internals;
-  occ_internals->loadBREP (fn.c_str());
-  occ_internals->buildLists ();
-  occ_internals->buildGModel (this);
+  occ_internals->loadBREP(fn.c_str());
+  occ_internals->buildLists();
+  occ_internals->buildGModel(this);
   return 1;
 }
+
 void GModel::deleteOCCInternals()
 {
-  if(occ_internals)delete occ_internals;
+  if(occ_internals) delete occ_internals;
 }
 
 #else
@@ -595,18 +497,21 @@ int GModel::readOCCSTEP(const std::string &fn)
       fn.c_str());
   return 0;
 }
+
 int GModel::readOCCIGES(const std::string &fn)
 {
   Msg(GERROR,"Gmsh has to be compiled with OpenCascade support to load %s",
       fn.c_str());
   return 0;
 }
+
 int GModel::readOCCBREP(const std::string &fn)
 {
   Msg(GERROR,"Gmsh has to be compiled with OpenCascade support to load %s",
       fn.c_str());
   return 0;
 }
+
 void GModel::deleteOCCInternals()
 {
 }
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 55229f7cf47c268c4cefe8a53b656eed4a931d3d..5ef10eb41723e3f93cd72c53f49b867fe035fe5b 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.14 2007-01-12 08:10:32 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.15 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -217,10 +217,9 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
     return V;
   }
 
-  int N, i, j;
+  int N, i;
   Vertex *v[5];
   double theta, t1, t2, t;
-  double vec[4];
   Vertex temp1, temp2;
 
   switch (c->Typ) {
@@ -392,79 +391,51 @@ void TransfiniteSph(Vertex S, Vertex center, Vertex * T)
 
 Vertex InterpolateRuledSurface(Surface * s, double u, double v)
 {
-  Vertex *c1, *c2, T, D[4], V[4], *S[4];
   Curve *C[4];
-  int i, issphere;
-  double eps = 1.e-6;
-
-  switch (s->Typ) {
-
-  case MSH_SURF_REGL:
-    issphere = 1;
-    for(i = 0; i < 4; i++) {
-      List_Read(s->Generatrices, i, &C[i]);
-      if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV) {
-        issphere = 0;
+  Vertex *c1, *c2;
+  int issphere = 1;
+  for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
+    List_Read(s->Generatrices, i, &C[i]);
+    if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV) {
+      issphere = 0;
+    }
+    else if(issphere) {
+      if(!i) {
+	List_Read(C[i]->Control_Points, 1, &c1);
       }
-      else if(issphere) {
-        if(!i) {
-          List_Read(C[i]->Control_Points, 1, &c1);
-        }
-        else {
-          List_Read(C[i]->Control_Points, 1, &c2);
-          if(compareVertex(&c1, &c2))
-            issphere = 0;
-        }
+      else {
+	List_Read(C[i]->Control_Points, 1, &c2);
+	if(compareVertex(&c1, &c2))
+	  issphere = 0;
       }
     }
+  }
 
+  Vertex *S[4], V[4], T;
+
+  if(s->Typ == MSH_SURF_REGL){
     S[0] = C[0]->beg;
     S[1] = C[1]->beg;
     S[2] = C[2]->beg;
     S[3] = C[3]->beg;
-
     V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
     V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
     V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
     V[3] = InterpolateCurve(C[3], C[3]->ubeg + (C[3]->uend - C[3]->ubeg) * (1. - v), 0);
-    
     T = TransfiniteQua(V[0], V[1], V[2], V[3], *S[0], *S[1], *S[2], *S[3], u, v);
-    if(issphere)
-      TransfiniteSph(*S[0], *c1, &T);
-    return T;
-    
-  case MSH_SURF_TRIC:
-    issphere = 1;
-    for(i = 0; i < 3; i++) {
-      List_Read(s->Generatrices, i, &C[i]);
-      if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV) {
-	issphere = 0;
-      }
-      else if(issphere) {
-	if(!i) {
-          List_Read(C[i]->Control_Points, 1, &c1);
-        }
-        else {
-          List_Read(C[i]->Control_Points, 1, &c2);
-          if(compareVertex(&c1, &c2))
-            issphere = 0;
-        }
-      }
-    }
-    
+  }
+  else{
     S[0] = C[0]->beg;
     S[1] = C[1]->beg;
     S[2] = C[2]->beg;
-
     V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
     V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
     V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1.-u), 0);
-    
     T = TransfiniteTri(V[0], V[1], V[2], *S[0], *S[1], *S[2], u, v);
-    if(issphere)
-      TransfiniteSph(*S[0], *c1, &T);    
-    return T;
   }
+
+  if(issphere) TransfiniteSph(*S[0], *c1, &T);
+  return T;
 }
 
 Vertex InterpolateExtrudedSurface(Surface * s, double u, double v)
diff --git a/Geo/GeoInterpolation.h b/Geo/GeoInterpolation.h
index 32ad551c139a0bf96526127f3e9b5eff90b025bd..45317748eac7eac9f3243bc5452d5ebb47933d2a 100644
--- a/Geo/GeoInterpolation.h
+++ b/Geo/GeoInterpolation.h
@@ -22,9 +22,7 @@
 
 #include "Geo.h"
 
-Vertex InterpolateCurve (Curve * Curve, double u, int derivee);
-
-Vertex InterpolateSurface (Surface * s, double u, double v, 
-                           int derivee, int u_v);
+Vertex InterpolateCurve(Curve *Curve, double u, int derivee);
+Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v);
 
 #endif
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 2c8ce8ed8915957dfe885cbd85b2f0ccf1669040..2be7850984e17dcf37367e939e1072a41b2c5526 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.26 2007-01-12 13:16:59 remacle Exp $
+// $Id: MElement.cpp,v 1.27 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -407,7 +407,7 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
 bool MTriangle::invertmappingXY(double *p, double *uv, double tol)
 {
   double mat[2][2];
-  double b[2], dum;
+  double b[2];
   getMat(mat);
   b[0] = p[0] - getVertex(0)->x();
   b[1] = p[1] - getVertex(0)->y();
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index c5136f5e37366f0a8e154e778a430bb467e55f28..219152728322fbcf01921edf71b1f3cc1c36c773 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.15 2006-12-11 20:12:33 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.16 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -42,20 +42,17 @@ OCCEdge::OCCEdge(GModel *model, TopoDS_Edge edge, int num, GVertex *v1, GVertex
 
 Range<double> OCCEdge::parBounds(int i) const
 { 
-  //  double a,b;
-  //  BRep_Tool::Range (c,a,b); 
   return(Range<double>(s0,s1));
 }
 
 void OCCEdge::setTrimmed (OCCFace *f)
 {
-  if (!trimmed)
-    {
-      trimmed = f;
-      const TopoDS_Face *s = (TopoDS_Face*) trimmed->getNativePtr();
-      curve2d = BRep_Tool::CurveOnSurface(c, *s, s0, s1);
-      if (curve2d.IsNull())  trimmed = 0;
-    }
+  if (!trimmed){
+    trimmed = f;
+    const TopoDS_Face *s = (TopoDS_Face*) trimmed->getNativePtr();
+    curve2d = BRep_Tool::CurveOnSurface(c, *s, s0, s1);
+    if (curve2d.IsNull())  trimmed = 0;
+  }
 }
 
 SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
@@ -63,46 +60,41 @@ SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
   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())
-    {
-      Msg(GERROR,"Reparam on face failed : curve %d is not on surface %d\n",tag(),face->tag());
-      return GEdge::reparamOnFace(face, epar,dir);
-    }
+  if(c2d.IsNull()){
+    Msg(GERROR,"Reparam on face failed : curve %d is not on surface %d\n",tag(),face->tag());
+    return GEdge::reparamOnFace(face, epar,dir);
+  }
   double u,v;
   c2d->Value(epar).Coord(u,v);
 
-  if (! isSeam ( face ) )
-    {
-      // sometimes OCC miserably fails ...
-      GPoint p1 = point(epar);
-      GPoint p2 = face->point(u,v);
-      const double dx = p1.x()-p2.x();
-      const double dy = p1.y()-p2.y();
-      const double dz = p1.z()-p2.z();
-      if (sqrt(dx*dx+dy*dy+dz*dz) > 1.e-7 * CTX.lc)
-	{
-	  GPoint ppp = face->closestPoint(SPoint3(p1.x(),p1.y(),p1.z()));
-	  return SPoint2(ppp.u(),ppp.v());
-	}
-      return SPoint2 (u,v);
+  if (!isSeam(face)){
+    // sometimes OCC miserably fails ...
+    GPoint p1 = point(epar);
+    GPoint p2 = face->point(u,v);
+    const double dx = p1.x()-p2.x();
+    const double dy = p1.y()-p2.y();
+    const double dz = p1.z()-p2.z();
+    if(sqrt(dx*dx+dy*dy+dz*dz) > 1.e-7 * CTX.lc){
+      GPoint ppp = face->closestPoint(SPoint3(p1.x(),p1.y(),p1.z()));
+      return SPoint2(ppp.u(),ppp.v());
     }
-  else
-    {
-      BRepAdaptor_Surface surface( *s );
-      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());
-      }
+    return SPoint2(u,v);
+  }
+  else{
+    BRepAdaptor_Surface surface(*s);
+    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. */
@@ -110,32 +102,27 @@ 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 );
-    }
+  if(surface.IsUPeriodic() || surface.IsVPeriodic()){
+    return BRep_Tool::IsClosed( c, *s );
+  }
   return 0;
 }
 
-
 GPoint OCCEdge::point(double par) const
 {
-  if (!curve.IsNull())
-    {
-      gp_Pnt pnt = curve->Value (par);
-      return GPoint(pnt.X(),pnt.Y(),pnt.Z());
-    }
-  else if (trimmed)
-    {
-      double u,v;
-      curve2d->Value(par).Coord(u,v);
-      return trimmed->point(u,v);
-    }
-  else
-    {
-      Msg(WARNING,"OCC Curve %d is neither a 3D curve not a trimmed curve",tag());
-      return GPoint (0,0,0);
-    }
+  if(!curve.IsNull()){
+    gp_Pnt pnt = curve->Value (par);
+    return GPoint(pnt.X(),pnt.Y(),pnt.Z());
+  }
+  else if(trimmed){
+    double u,v;
+    curve2d->Value(par).Coord(u,v);
+    return trimmed->point(u,v);
+  }
+  else{
+    Msg(WARNING,"OCC Curve %d is neither a 3D curve not a trimmed curve",tag());
+    return GPoint (0,0,0);
+  }
 }
 
 GPoint OCCEdge::closestPoint(const SPoint3 & qp)
@@ -165,47 +152,41 @@ double OCCEdge::parFromPoint(const SPoint3 &pt) const
 
 GEntity::GeomType OCCEdge::geomType() const
 {
-  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 (curve2d->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve))
-	return BSpline;
-      else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_BezierCurve))
-	return Bezier;
-      //   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 (curve->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve))
-	return BSpline;
-      else if (curve->DynamicType() == STANDARD_TYPE(Geom_BezierCurve))
-	return Bezier;
-      //   else if (occface->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
-      //     return Cone;
-      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 (curve2d->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve))
+      return BSpline;
+    else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_BezierCurve))
+      return Bezier;
+    //   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 (curve->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve))
+      return BSpline;
+    else if (curve->DynamicType() == STANDARD_TYPE(Geom_BezierCurve))
+      return Bezier;
+    //   else if (occface->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
+    //     return Cone;
+    return Unknown;
+  }
 }
 
 int OCCEdge::minimumMeshSegments () const
 {
-//   if(geomType() == Circle || geomType() == Ellipse)
-//     return (int)(fabs(s1 - s0) *
-// 		 (double)CTX.mesh.min_circ_points / Pi) - 1;
-//   else
-    return GEdge::minimumMeshSegments () ;
+  return GEdge::minimumMeshSegments () ;
 }
 
 int OCCEdge::minimumDrawSegments () const
@@ -220,43 +201,41 @@ int OCCEdge::minimumDrawSegments () const
     return 20 * n;
 }
 
-
 double OCCEdge::curvature(double par) const 
 {
   const double eps = 1.e-15;
   Standard_Real Crv;
-  if (curve.IsNull())
-    {
-      Geom2dLProp_CLProps2d aCLProps(curve2d, 2, eps);
-      aCLProps.SetParameter (par);
-      if(!aCLProps.IsTangentDefined())
-	Crv =eps;
-      else
-	Crv = aCLProps.Curvature();
-    }
-  else
-    {
-      BRepAdaptor_Curve brepc(c);
-      BRepLProp_CLProps prop(brepc, 2, eps);
-      prop.SetParameter (par); 
-      if (!prop.IsTangentDefined())
-	Crv = eps;
-      else
-	Crv = prop.Curvature();
-    }
+  if (curve.IsNull()){
+    Geom2dLProp_CLProps2d aCLProps(curve2d, 2, eps);
+    aCLProps.SetParameter (par);
+    if(!aCLProps.IsTangentDefined())
+      Crv =eps;
+    else
+      Crv = aCLProps.Curvature();
+  }
+  else{
+    BRepAdaptor_Curve brepc(c);
+    BRepLProp_CLProps prop(brepc, 2, eps);
+    prop.SetParameter (par); 
+    if (!prop.IsTangentDefined())
+      Crv = eps;
+    else
+      Crv = prop.Curvature();
+  }
   if (Crv <= eps)Crv = eps;
-
-//   std::list<GFace*> ff = faces();
-//   std::list<GFace *>::iterator it =  ff.begin();
-//   while (it != ff.end())
-//     {
-//       SPoint2 par2 = reparamOnFace((*it),par,1);
-//       const double cc = (*it)->curvature ( par2 );
-//       if (cc > 0)
-//     Crv = std::max( Crv, cc);  
-//     ++it;
-// }  
-// printf("curvature = %12.5E\n",Crv); 
-return Crv;
+  
+  // std::list<GFace*> ff = faces();
+  // std::list<GFace *>::iterator it =  ff.begin();
+  // while (it != ff.end()){
+  //   SPoint2 par2 = reparamOnFace((*it),par,1);
+  //   const double cc = (*it)->curvature ( par2 );
+  //   if (cc > 0)
+  //     Crv = std::max( Crv, cc);  
+  //   ++it;
+  // }  
+  // printf("curvature = %12.5E\n",Crv); 
+
+  return Crv;
 }
+
 #endif
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index cc7257fb15baddf6d45a61bdc537af30a9516ab4..c407ba3723e7d799d2bfe6d13b14b7629a90d84e 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.17 2006-12-20 15:50:57 remacle Exp $
+// $Id: OCCFace.cpp,v 1.18 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -39,43 +39,38 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
   : GFace(m, num), s(_s)
 {
   TopExp_Explorer exp0, exp01, exp1, exp2, exp3;
-  for (exp2.Init (s, TopAbs_WIRE); exp2.More(); exp2.Next())
-    {
-      TopoDS_Shape wire = exp2.Current();
-      Msg(INFO,"OCC Face %d - New Wire",num);
-      std::list<GEdge*> l_wire;
-      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);
-	  e->addFace(this);
-	  if (!e->is3D())
-	    {
-	      OCCEdge *occe = (OCCEdge*)e;
-	      occe->setTrimmed(this);
-	    }
-	}      
-
-      GEdgeLoop el (l_wire);
-
-      for (GEdgeLoop::citer it = el.begin() ; it != el.end() ; ++it)
-	{
-	  l_edges.push_back(it->ge);
-	  l_dirs.push_back(it->_sign);
-	}
-      
-      edgeLoops.push_back(el);
+  for(exp2.Init(s, TopAbs_WIRE); exp2.More(); exp2.Next()){
+    TopoDS_Shape wire = exp2.Current();
+    Msg(INFO,"OCC Face %d - New Wire",num);
+    std::list<GEdge*> l_wire;
+    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);
+      e->addFace(this);
+      if(!e->is3D()){
+	OCCEdge *occe = (OCCEdge*)e;
+	occe->setTrimmed(this);
+      }
+    }      
+    
+    GEdgeLoop el(l_wire);
+    
+    for(GEdgeLoop::citer it = el.begin() ; it != el.end() ; ++it){
+      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);
+  Msg(INFO, "OCC Face %d with %d edges", num, l_edges.size());
+  ShapeAnalysis::GetFaceUVBounds(s, umin, umax, vmin, vmax);
   // we do that for the projections to converge on the 
   // borders of the surface
   umin -= fabs(umax-umin)/100.0;
@@ -89,8 +84,8 @@ Range<double> OCCFace::parBounds(int i) const
 {  
   double umin2, umax2, vmin2, vmax2;
 
-  ShapeAnalysis::GetFaceUVBounds (s, umin2, umax2, vmin2, vmax2);
-  if (i==0)
+  ShapeAnalysis::GetFaceUVBounds(s, umin2, umax2, vmin2, vmax2);
+  if(i == 0)
     return Range<double>(umin2, umax2);
   return Range<double>(vmin2, vmax2);
 }
@@ -100,13 +95,13 @@ SVector3 OCCFace::normal(const SPoint2 &param) const
   gp_Pnt pnt;
   gp_Vec du, dv;
 
-  occface->D1(param.x(), param.y(),pnt,du,dv);
+  occface->D1(param.x(), param.y(), pnt, du, dv);
 
-  SVector3 t1 (du.X(),du.Y(),du.Z());
-  SVector3 t2 (dv.X(),dv.Y(),dv.Z());
-  SVector3 n ( crossprod (t1,t2));
+  SVector3 t1(du.X(), du.Y(), du.Z());
+  SVector3 t2(dv.X(), dv.Y(), dv.Z());
+  SVector3 n(crossprod(t1, t2));
   n.normalize();
-  if (s.Orientation() == TopAbs_REVERSED) return n*(-1.);  
+  if (s.Orientation() == TopAbs_REVERSED) return n * (-1.);  
   return n;  
 }
 
@@ -127,24 +122,23 @@ GPoint OCCFace::point(const SPoint2 &pt) const
 
 GPoint OCCFace::point(double par1, double par2) const
 {
-  double pp[2] = {par1,par2};
-  gp_Pnt val = occface->Value (par1, par2);
-  return GPoint(val.X(),val.Y(),val.Z(), this, pp);
+  double pp[2] = {par1, par2};
+  gp_Pnt val = occface->Value(par1, par2);
+  return GPoint(val.X(), val.Y(), val.Z(), this, pp);
 }
 
 GPoint OCCFace::closestPoint(const SPoint3 & qp) const
 {
-  gp_Pnt pnt(qp.x(),qp.y(),qp.z());
+  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");
-      return GPoint(0,0);
-    }
+  if(!proj.NbPoints()){
+    Msg(GERROR,"OCC Project Point on Surface FAIL");
+    return GPoint(0, 0);
+  }
   pnt = proj.NearestPoint();
   double pp[2];
   proj.LowerDistanceParameters (pp[0], pp[1]);
-  return GPoint(pnt.X(),pnt.Y(),pnt.Z(), this, pp);
+  return GPoint(pnt.X(), pnt.Y(), pnt.Z(), this, pp);
 }
 
 int OCCFace::containsParam(const SPoint2 &pt) const
@@ -160,16 +154,15 @@ int OCCFace::containsParam(const SPoint2 &pt) const
 
 SPoint2 OCCFace::parFromPoint(const SPoint3 &qp) const
 {
-  gp_Pnt pnt(qp.x(),qp.y(),qp.z());
+  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");
-      return GFace::parFromPoint(qp);
-    }   
+  if(!proj.NbPoints()){
+    Msg(GERROR,"OCC Project Point on Surface FAIL");
+    return GFace::parFromPoint(qp);
+  }   
   double U,V;
   proj.LowerDistanceParameters (U, V);
-  return SPoint2(U,V);
+  return SPoint2(U, V);
 }
 
 GEntity::GeomType OCCFace::geomType() const
@@ -180,14 +173,13 @@ GEntity::GeomType OCCFace::geomType() const
     return Cylinder;
   else if (occface->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
     return Cone;
-   else if (occface->DynamicType() == STANDARD_TYPE(Geom_SphericalSurface))
-     return Sphere;
-   else if (occface->DynamicType() == STANDARD_TYPE(Geom_BSplineSurface))
-     return BSplineSurface;
+  else if (occface->DynamicType() == STANDARD_TYPE(Geom_SphericalSurface))
+    return Sphere;
+  else if (occface->DynamicType() == STANDARD_TYPE(Geom_BSplineSurface))
+    return BSplineSurface;
   return Unknown;
 }
 
-
 double OCCFace::curvature (const SPoint2 &param) const
 {
   const double eps = 1.e-12;
@@ -195,59 +187,56 @@ double OCCFace::curvature (const SPoint2 &param) const
   BRepLProp_SLProps prop(sf, 2, eps);
   prop.SetParameters (param.x(),param.y());
 
-  if (!prop.IsCurvatureDefined())
-    {
-      return eps;
-    }
+  if (!prop.IsCurvatureDefined()){
+    return eps;
+  }
   return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
 }
 
-
 int OCCFace::containsPoint(const SPoint3 &pt) const
 { 
-  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.;
-      double v[3] = {pt.x(), pt.y(), 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() ));
-	  double v1[3] = {pp1.x(), pp1.y(), pp1.z()};
-	  double v2[3] = {pp2.x(), pp2.y(), pp2.z()};
-	  angle += angle_plan(v, v1, v2, n);
-	}
+  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.;
+    double v[3] = {pt.x(), pt.y(), 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() ));
+	double v1[3] = {pp1.x(), pp1.y(), pp1.z()};
+	double v2[3] = {pp2.x(), pp2.y(), 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;
     }
+    // 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 ...");
   return false;
 }
+
 // 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);
-//     }  
+//   for (int j = 1; j <= ntriangles; j++){
+//     Poly_Triangle triangle = (triangulation -> Triangles())(j);
+//   }  
 // }
 
 #endif
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index 99b73dffcf476ee1afd873436eeb8f4a2f227cfe..d0b82fc46b7fb1344634c2376920bc6521151c2c 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCRegion.cpp,v 1.4 2006-11-27 22:22:14 geuzaine Exp $
+// $Id: OCCRegion.cpp,v 1.5 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -32,21 +32,19 @@ OCCRegion::OCCRegion(GModel *m, TopoDS_Solid _s, int num, TopTools_IndexedMapOfS
   : GRegion(m, num), s(_s)
 {
   TopExp_Explorer exp0, exp01, exp1, exp2, exp3;
-  for (exp2.Init (s, TopAbs_SHELL); exp2.More(); exp2.Next())
-    {
-      TopoDS_Shape shell = exp2.Current();
-      Msg(INFO,"OCC Region %d - New Shell",num);
-      for (exp3.Init (shell, TopAbs_FACE); exp3.More(); exp3.Next())
-	{	  
-	  TopoDS_Face face = TopoDS::Face (exp3.Current());
-	  int index = fmap.FindIndex(face);
-	  GFace *f = m->faceByTag(index);
-	  if(!f) throw;
-	  l_faces.push_back(f);
-	  f->addRegion(this);
-	}      
-    }
-  Msg(INFO,"OCC Region %d with %d edges",num,l_faces.size());
+  for(exp2.Init(s, TopAbs_SHELL); exp2.More(); exp2.Next()){
+    TopoDS_Shape shell = exp2.Current();
+    Msg(INFO,"OCC Region %d - New Shell",num);
+    for(exp3.Init(shell, TopAbs_FACE); exp3.More(); exp3.Next()){	  
+      TopoDS_Face face = TopoDS::Face(exp3.Current());
+      int index = fmap.FindIndex(face);
+      GFace *f = m->faceByTag(index);
+      if(!f) throw;
+      l_faces.push_back(f);
+      f->addRegion(this);
+    }      
+  }
+  Msg(INFO, "OCC Region %d with %d edges", num, l_faces.size());
 }
 
 GEntity::GeomType OCCRegion::geomType() const
diff --git a/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index 84294282ff9acf37183c6b97bdab4beae6c87fbf..fe92a1eb88005e696f0aac4e6fe53b61b499d2ba 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCVertex.cpp,v 1.9 2006-11-29 16:57:01 remacle Exp $
+// $Id: OCCVertex.cpp,v 1.10 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -26,76 +26,65 @@
 
 #if defined(HAVE_OCC)
 
-double max_surf_curvature ( const GVertex *gv, double x, double y, double z , const GEdge *_myGEdge)
+double max_surf_curvature(const GVertex *gv, double x, double y, double z, const GEdge *_myGEdge)
 {
   std::list<GFace *> faces = _myGEdge->faces();
   std::list<GFace *>::iterator it =  faces.begin();
   double curv = 1.e-22;
-  while (it != faces.end())
-    {
-      SPoint2 par = gv->reparamOnFace((*it),1);
-      double cc = (*it)->curvature ( par );
-      if (cc > 0)curv = std::max(curv, cc );					      
-      ++it;
-    }  
+  while(it != faces.end()){
+    SPoint2 par = gv->reparamOnFace((*it),1);
+    double cc = (*it)->curvature(par);
+    if(cc > 0) curv = std::max(curv, cc);					      
+    ++it;
+  }  
   return curv;
 }
 
 double OCCVertex::max_curvature_of_surfaces() const
 {  
-  if (max_curvature <0)
-    {
-      for (std::list<GEdge*> :: const_iterator it = l_edges.begin() ; it != l_edges.end() ; ++it )
-	{
-	  max_curvature = std::max ( max_surf_curvature (this, x(), y(), z(), *it) , max_curvature);  
-	}
-      //      printf("max curvature (%d) = %12.5E lc = %12.5E\n",tag(),max_curvature,prescribedMeshSizeAtVertex());
-
+  if(max_curvature < 0){
+    for(std::list<GEdge*>::const_iterator it = l_edges.begin(); it != l_edges.end(); ++it ){
+      max_curvature = std::max(max_surf_curvature(this, x(), y(), z(), *it), max_curvature);
     }
+    // printf("max curvature (%d) = %12.5E lc = %12.5E\n",tag(),max_curvature,prescribedMeshSizeAtVertex());
+  }
   return max_curvature;
 }
 
-SPoint2 OCCVertex::reparamOnFace ( GFace *gf , int dir) 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;
-    }  
+  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();
-	  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;
+  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();
+      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;
+  }
 }
 
-
-
 #endif
diff --git a/Geo/projectionFace.cpp b/Geo/projectionFace.cpp
index d59a5124b7e391fb7818690d7f08563aac41fa44..5cef9ce07d0eabfbb4b9b34423cb30d248468e96 100644
--- a/Geo/projectionFace.cpp
+++ b/Geo/projectionFace.cpp
@@ -136,11 +136,6 @@ SPoint2 parabolicCylinder::parFromPoint(const SPoint3 &p) const
   pos = rotatePoint(pos,rot);
   pos = scalePoint(pos, scalar);
   
-  // since the y co-ordinate is completely dependent on u, I actually
-  // don't need it to compute the u-v point?
-  double u = pos[0] + .5;
-  double v = pos[1] + .5;
-  
   SPoint2 q(pos[0],pos[1]);	
   return q;
 }
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 588bdaf7b2fb33a25833cf52227034b07731bcac..fe4c0ce920a2eab7ca56be6f3626c7f87a744ea3 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.71 2006-12-11 20:12:33 remacle Exp $
+// $Id: BDS.cpp,v 1.72 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -802,8 +802,8 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
 int BDS_Edge :: numTriangles() const 
 {
   int NT = 0;
-  for (int i=0;i<_faces.size();i++) 
-    if (faces (i)->numEdges() == 3) NT++ ;
+  for (unsigned int i=0;i<_faces.size();i++) 
+    if (faces(i)->numEdges() == 3) NT++ ;
   return NT;
 }
 
@@ -895,9 +895,6 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
   if(o->g != p->g)
     return false;
 
-  double u = p->u;
-  double v = p->v;
-
   // printf("collapsing an edge :");
   // print_edge(e);
 
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index e62e5b48f1e44e8e15d0961fcfaff6603cbb1ce9..b5c5decff454c87e3eec8abb8291f3a77e7029b3 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -361,6 +361,7 @@ class BDS_SwapEdgeTest
  public:
   virtual bool operator() (BDS_Point *p1,BDS_Point *p2,
 			   BDS_Point *q1,BDS_Point *q2) const = 0; 
+  virtual ~BDS_SwapEdgeTest(){}
 };
 
 class BDS_SwapEdgeTestParametric : public BDS_SwapEdgeTest
@@ -368,6 +369,7 @@ class BDS_SwapEdgeTestParametric : public BDS_SwapEdgeTest
  public:
   virtual bool operator() (BDS_Point *p1,BDS_Point *p2,
 			   BDS_Point *q1,BDS_Point *q2) const ; 
+  virtual ~BDS_SwapEdgeTestParametric(){}
 };
 
 class BDS_Mesh 
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index aeb38945894c44cc80d35626d5e996af3be85493..8d91a2707b70610414808696bd01c5b93633fb53 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.11 2007-01-12 13:16:59 remacle Exp $
+// $Id: BackgroundMesh.cpp,v 1.12 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -99,7 +99,7 @@ double LC_MVertex_CURV ( GEntity *ge, double U, double V )
       {
 	GEdge *ged = (GEdge *)ge;
 	//Crv = ged->curvature(  U );
-	Crv = max_surf_curvature ( (const GEdge *)ge, U);
+	Crv = max_surf_curvature ( ged, U);
       }
       break;
     case 2:
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index afe33c4cc54da2b710a5ba1262ba552cb020c599..6d4c6fe6449cdd6ae4aeabdb5c79c5e476dd3126 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -1,4 +1,4 @@
-// $Id: DivideAndConquer.cpp,v 1.8 2007-01-08 16:42:42 geuzaine Exp $
+// $Id: DivideAndConquer.cpp,v 1.9 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -531,8 +531,6 @@ int Conversion(DocPeek doc)
   int n, i, j;
   int count = 0, count2 = 0;
   PointNumero aa, bb, cc;
-  PointRecord *ptemp;
-
   PointRecord *gPointArray = doc->points;
 
   n = doc->numPoints;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index b4f9ecdb555ed7c54adc3ca65599841c6f646b34..916b28f3e7eabe1668272bd442bd9cbcc12b82f6 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.111 2006-12-17 12:44:27 geuzaine Exp $
+// $Id: Generator.cpp,v 1.112 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -202,7 +202,9 @@ void Mesh1D()
   if(TooManyElements(1)) return;
   Msg(STATUS1, "Meshing 1D...");
   double t1 = Cpu();
+
   std::for_each(GMODEL->firstEdge(), GMODEL->lastEdge(), meshGEdge());
+
   double t2 = Cpu();
   CTX.mesh_timer[0] = t2 - t1;
   Msg(INFO, "Mesh 1D complete (%g s)", CTX.mesh_timer[0]);
@@ -214,7 +216,9 @@ void Mesh2D()
   if(TooManyElements(2)) return;
   Msg(STATUS1, "Meshing 2D...");
   double t1 = Cpu();
+
   std::for_each(GMODEL->firstFace(), GMODEL->lastFace(), meshGFace());
+
   double t2 = Cpu();
   CTX.mesh_timer[1] = t2 - t1;
   Msg(INFO, "Mesh 2D complete (%g s)", CTX.mesh_timer[1]);
@@ -226,7 +230,15 @@ void Mesh3D()
   if(TooManyElements(3)) return;
   Msg(STATUS1, "Meshing 3D...");
   double t1 = Cpu();
+
+  // mesh the extruded volumes first
+  std::for_each(GMODEL->firstRegion(), GMODEL->lastRegion(), meshGRegionExtruded());
+  // then subdivide if necessary (unfortunately the subdivision is a
+  // global operation, which can require changing the surface mesh!)
+  SubdivideExtrudedMesh(GMODEL);
+  // then mesh the rest
   std::for_each(GMODEL->firstRegion(), GMODEL->lastRegion(), meshGRegion());
+
   double t2 = Cpu();
   CTX.mesh_timer[2] = t2 - t1;
   Msg(INFO, "Mesh 3D complete (%g s)", CTX.mesh_timer[2]);
@@ -237,7 +249,9 @@ void OptimizeMesh()
 {
   Msg(STATUS1, "Optimizing 3D...");
   double t1 = Cpu();
+
   std::for_each(GMODEL->firstRegion(), GMODEL->lastRegion(), optimizeMeshGRegion());
+
   double t2 = Cpu();
   Msg(INFO, "Mesh 3D optimization complete (%g s)", t2 - t1);
   Msg(STATUS1, "Mesh");
diff --git a/Mesh/meshGEdgeExtruded.cpp b/Mesh/meshGEdgeExtruded.cpp
index 9fb967856982e67cc673da48f7a9d984cd42ee2e..58f6560f1bc7040366e7562719af875ca666c11b 100644
--- a/Mesh/meshGEdgeExtruded.cpp
+++ b/Mesh/meshGEdgeExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdgeExtruded.cpp,v 1.5 2006-11-27 22:22:17 geuzaine Exp $
+// $Id: meshGEdgeExtruded.cpp,v 1.6 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -76,7 +76,10 @@ int MeshExtrudedCurve(GEdge *ge)
   else {
     // curve is a copy of another curve (the "top" of the extrusion)
     GEdge *from = ge->model()->edgeByTag(std::abs(ep->geo.Source));
-    if(!from) return 0;
+    if(!from){
+      Msg(GERROR, "Unknown source curve %d for extrusion", ep->geo.Source);
+      return 0;
+    }
     copyMesh(from, ge);
   }
 
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 3ad90a14c93d9b54d25f2e292cfe5c077cc6da99..2ae5a72370c33828d1092a584c9565d5f4fbeefa 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.47 2007-01-12 19:47:52 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.48 2007-01-16 11:31:41 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -795,7 +795,7 @@ bool gmsh2DMeshGenerator ( GFace *gf )
       }
   }
 
-   char name[245];
+  //char name[245];
    //sprintf(name,"param%d.pos",gf->tag());
    //outputScalarField(m->triangles, name,1);
    //   sprintf(name,"real%d.pos",gf->tag());
@@ -1088,11 +1088,11 @@ bool buildConsecutiveListOfVertices (  GFace *gf,
      result.insert(result.end(),edgeLoop_BDS.begin(),edgeLoop_BDS.end());	         
    }
 
-  if (gf->tag() == 280)
-    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);
-      }
+//   if (gf->tag() == 280)
+//     for (unsigned 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);
+//       }
 
   return true;
 }
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index 0892a9fdf89c75a7512cea99cf43c43bb367d125..3cb51edea00160107370a98a93c10dd9cf366170 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -21,6 +21,7 @@
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include <vector>
+#include <set>
 
 class MVertex;
 class GFace;
@@ -55,6 +56,7 @@ void computeEdgeLoops(const GFace *gf,
 		      std::vector<int> &indices);
 
 int MeshTransfiniteSurface(GFace *gf);
-int MeshExtrudedSurface(GFace *gf);
+int MeshExtrudedSurface(GFace *gf, 
+			std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges=0);
 
 #endif
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index d85388f3c3f98c8f9dc7085ad27d7912b81775fa..91d0fdaeda091a3386c465d37badb35a7d2b4ab1 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceExtruded.cpp,v 1.14 2006-12-03 02:05:47 geuzaine Exp $
+// $Id: meshGFaceExtruded.cpp,v 1.15 2007-01-16 11:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -27,34 +27,58 @@
 
 extern Context_T CTX;
 
-void createQuaTri(std::vector<MVertex*> &v, GFace *to)
+void createQuaTri(std::vector<MVertex*> &v, GFace *to,
+		  std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges)
 {
+  ExtrudeParams *ep = to->meshAttributes.extrude;
+
   if(v[0] == v[1] || v[1] == v[3])
     to->triangles.push_back(new MTriangle(v[0], v[3], v[2]));
   else if(v[0] == v[2] || v[2] == v[3])
     to->triangles.push_back(new MTriangle(v[0], v[1], v[3]));
   else if(v[0] == v[3] || v[1] == v[2])
     Msg(GERROR, "Uncoherent extruded quadrangle in surface %d", to->tag());
-  else
-    to->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[3], v[2]));
+  else{
+    if(ep->mesh.Recombine){
+      to->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[3], v[2]));
+    }
+    else if(!constrainedEdges){
+      to->triangles.push_back(new MTriangle(v[0], v[1], v[3]));
+      to->triangles.push_back(new MTriangle(v[0], v[3], v[2]));
+    }
+    else{
+      std::pair<MVertex*, MVertex*> p(std::min(v[1], v[2]), std::max(v[1], v[2]));
+      if(constrainedEdges->count(p)){
+	to->triangles.push_back(new MTriangle(v[2], v[1], v[0]));
+	to->triangles.push_back(new MTriangle(v[2], v[3], v[1]));
+      }
+      else {
+	to->triangles.push_back(new MTriangle(v[2], v[3], v[0]));
+	to->triangles.push_back(new MTriangle(v[0], v[3], v[1]));
+      }
+    }
+  }
 }
 
 void extrudeMesh(GEdge *from, GFace *to,
-		 std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+		 std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+		 std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges)
 {
   ExtrudeParams *ep = to->meshAttributes.extrude;
 
-  // create vertices
-  for(unsigned int i = 0; i < from->mesh_vertices.size(); i++){
-    MVertex *v = from->mesh_vertices[i];
-    for(int j = 0; j < ep->mesh.NbLayer; j++) {
-      for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
-	double x = v->x(), y = v->y(), z = v->z();
-	ep->Extrude(j, k + 1, x, y, z);
-	if(j != ep->mesh.NbLayer - 1 || k != ep->mesh.NbElmLayer[j] - 1){
-	  MVertex *newv = new MVertex(x, y, z, to);
-	  to->mesh_vertices.push_back(newv);
-	  pos.insert(newv);
+  // create vertices (if the edges are constrained, they already exist)
+  if(!constrainedEdges){
+    for(unsigned int i = 0; i < from->mesh_vertices.size(); i++){
+      MVertex *v = from->mesh_vertices[i];
+      for(int j = 0; j < ep->mesh.NbLayer; j++) {
+	for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
+	  double x = v->x(), y = v->y(), z = v->z();
+	  ep->Extrude(j, k + 1, x, y, z);
+	  if(j != ep->mesh.NbLayer - 1 || k != ep->mesh.NbElmLayer[j] - 1){
+	    MVertex *newv = new MVertex(x, y, z, to);
+	    to->mesh_vertices.push_back(newv);
+	    pos.insert(newv);
+	  }
 	}
       }
     }
@@ -86,7 +110,7 @@ void extrudeMesh(GEdge *from, GFace *to,
 	  }
 	  verts.push_back(*itp);
 	}
-	createQuaTri(verts, to);
+	createQuaTri(verts, to, constrainedEdges);
       }
     }
   }
@@ -144,20 +168,14 @@ void copyMesh(GFace *from, GFace *to,
   }
 }
 
-int MeshExtrudedSurface(GFace *gf)
+int MeshExtrudedSurface(GFace *gf, 
+			std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges)
 {
   ExtrudeParams *ep = gf->meshAttributes.extrude;
 
   if(!ep || !ep->mesh.ExtrudeMesh)
     return 0;
 
-  // FIXME
-  static bool warn_recombine = 1;
-  if(!ep->mesh.Recombine && warn_recombine){
-    Msg(WARNING, "Non-recombined surface extrusion has not been reimplemented yet");
-    warn_recombine = false;
-  }
-
   // build a set with all the vertices on the boundary of gf
   double old_tol = MVertexLessThanLexicographic::tolerance; 
   MVertexLessThanLexicographic::tolerance = 1.e-12 * CTX.lc;
@@ -173,16 +191,27 @@ int MeshExtrudedSurface(GFace *gf)
     ++it;
   }
 
+  // if the edges are constrained, the vertices already exist on the
+  // face--so we add them to the set
+  if(constrainedEdges)
+    pos.insert(gf->mesh_vertices.begin(), gf->mesh_vertices.end());
+  
   if(ep->geo.Mode == EXTRUDED_ENTITY) {
     // surface is extruded from a curve
     GEdge *from = gf->model()->edgeByTag(std::abs(ep->geo.Source));
-    if(!from) return 0;
-    extrudeMesh(from, gf, pos);
+    if(!from){
+      Msg(GERROR, "Unknown source curve %d for extrusion", ep->geo.Source);
+      return 0;
+    }
+    extrudeMesh(from, gf, pos, constrainedEdges);
   }
   else {
     // surface is a copy of another surface (the "top" of the extrusion)
     GFace *from = gf->model()->faceByTag(std::abs(ep->geo.Source));
-    if(!from) return 0;
+    if(!from){ 
+      Msg(GERROR, "Unknown source surface %d for extrusion", ep->geo.Source);
+      return 0;
+    }
     copyMesh(from, gf, pos);
   }
 
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 52040844139eca6fcd4ab7cfd886117548b20935..2beeccc3e4c80645c33138b44de8fb6f52f73967 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.22 2007-01-12 19:47:52 geuzaine Exp $
+// $Id: meshGRegion.cpp,v 1.23 2007-01-16 11:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -431,6 +431,10 @@ void meshGRegion::operator() (GRegion *gr)
 {  
   if(gr->geomType() == GEntity::DiscreteVolume) return;
 
+  ExtrudeParams *ep = gr->meshAttributes.extrude;
+
+  if(ep && ep->mesh.ExtrudeMesh) return;
+
   // Send a messsage to the GMSH environment
   Msg(STATUS2, "Meshing volume %d", gr->tag());
 
@@ -439,7 +443,6 @@ void meshGRegion::operator() (GRegion *gr)
   dem(gr);
 
   if(MeshTransfiniteVolume(gr)) return;
-  if(MeshExtrudedVolume(gr)) return;
 
   std::list<GFace*> faces = gr->faces();
 
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index 2f93cf5140cde44128605033e06dae7b98f3d767..9c15ea346660e4f281a20e3f31587749283f7fe6 100644
--- a/Mesh/meshGRegion.h
+++ b/Mesh/meshGRegion.h
@@ -20,6 +20,7 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
+class GModel;
 class GRegion;
 
 // Create the mesh of the region
@@ -28,6 +29,11 @@ class meshGRegion {
   void operator () (GRegion *);
 };
 
+class meshGRegionExtruded {
+ public :
+  void operator () (GRegion *);
+};
+
 // Optimize the mesh of the region
 class optimizeMeshGRegion {
  public :
@@ -36,11 +42,11 @@ class optimizeMeshGRegion {
 
 // destroy the mesh of the region
 class deMeshGRegion {
-public :
+ public :
   void operator () (GRegion *);
 };
 
 int MeshTransfiniteVolume(GRegion *gr);
-int MeshExtrudedVolume(GRegion *gr);
+int SubdivideExtrudedMesh(GModel *m);
 
 #endif
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index 8bf3b05be23decf13cb3be21973e33b80821a21f..c7cae7449cc4987a7399e662da8335664a35a53d 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionExtruded.cpp,v 1.7 2006-12-03 02:05:47 geuzaine Exp $
+// $Id: meshGRegionExtruded.cpp,v 1.8 2007-01-16 11:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -22,6 +22,8 @@
 #include <set>
 #include "ExtrudeParams.h"
 #include "GModel.h"
+#include "meshGFace.h"
+#include "meshGRegion.h"
 #include "Context.h"
 #include "Message.h"
 
@@ -83,6 +85,40 @@ void createHexPri(std::vector<MVertex*> &v, GRegion *to)
   }
 }
 
+void createTet(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, GRegion *to)
+{
+  if(v1 != v2 && v1 != v3 && v1 != v4 && v2 != v3 && v2 != v4 && v3 != v4)
+    to->tetrahedra.push_back(new MTetrahedron(v1, v2, v3, v4));
+}
+
+int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k, 
+			std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+			std::vector<MVertex*> &verts)
+{
+  std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
+  double x[8], y[8], z[8];
+  int n = ele->getNumVertices();
+  for(int p = 0; p < n; p++){
+    MVertex *v = ele->getVertex(p);
+    x[p] = x[p + n] = v->x();
+    y[p] = y[p + n] = v->y();
+    z[p] = z[p + n] = v->z();
+  }
+  for(int p = 0; p < n; p++){
+    ep->Extrude(j, k, x[p], y[p], z[p]);
+    ep->Extrude(j, k + 1, x[p + n], y[p + n], z[p + n]);
+  }
+  for(int p = 0; p < 2 * n; p++){
+    MVertex tmp(x[p], y[p], z[p], 0, -1);
+    itp = pos.find(&tmp);
+    if(itp == pos.end())
+      Msg(GERROR, "Could not find extruded vertex");
+    else
+      verts.push_back(*itp);
+  }
+  return verts.size();
+}
+
 void extrudeMesh(GFace *from, GRegion *to,
 		 std::set<MVertex*, MVertexLessThanLexicographic> &pos)
 {
@@ -107,91 +143,35 @@ void extrudeMesh(GFace *from, GRegion *to,
   // create elements (note that it would be faster to access the
   // *interior* nodes by direct indexing, but it's just simpler to
   // query everything by position)
-  std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
   for(unsigned int i = 0; i < from->triangles.size(); i++){
-    MVertex *v0 = from->triangles[i]->getVertex(0);
-    MVertex *v1 = from->triangles[i]->getVertex(1);
-    MVertex *v2 = from->triangles[i]->getVertex(2);
     for(int j = 0; j < ep->mesh.NbLayer; j++) {
       for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
 	std::vector<MVertex*> verts;
-	double x[6] = {v0->x(), v1->x(), v2->x(), v0->x(), v1->x(), v2->x()};
-	double y[6] = {v0->y(), v1->y(), v2->y(), v0->y(), v1->y(), v2->y()};
-	double z[6] = {v0->z(), v1->z(), v2->z(), v0->z(), v1->z(), v2->z()};
-	for(int p = 0; p < 3; p++){
-	  ep->Extrude(j, k, x[p], y[p], z[p]);
-	  ep->Extrude(j, k + 1, x[p + 3], y[p + 3], z[p + 3]);
-	}
-	for(int p = 0; p < 6; p++){
-	  MVertex tmp(x[p], y[p], z[p], 0, -1);
-	  itp = pos.find(&tmp);
-	  if(itp == pos.end()) {
-	    Msg(GERROR, "Could not find extruded vertex in volume %d", to->tag());
-	    return;
-	  }
-	  verts.push_back(*itp);
-	}
-	createPriPyrTet(verts, to);
+	if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, verts) == 6)
+	  createPriPyrTet(verts, to);
       }
     }
   }
   for(unsigned int i = 0; i < from->quadrangles.size(); i++){
-    MVertex *v0 = from->quadrangles[i]->getVertex(0);
-    MVertex *v1 = from->quadrangles[i]->getVertex(1);
-    MVertex *v2 = from->quadrangles[i]->getVertex(2);
-    MVertex *v3 = from->quadrangles[i]->getVertex(3);
     for(int j = 0; j < ep->mesh.NbLayer; j++) {
       for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
 	std::vector<MVertex*> verts;
-	double x[8] = {v0->x(), v1->x(), v2->x(), v3->x(), 
-		       v0->x(), v1->x(), v2->x(), v3->x()};
-	double y[8] = {v0->y(), v1->y(), v2->y(), v3->y(), 
-		       v0->y(), v1->y(), v2->y(), v3->y()};
-	double z[8] = {v0->z(), v1->z(), v2->z(), v3->z(), 
-		       v0->z(), v1->z(), v2->z(), v3->z()};
-	for(int p = 0; p < 4; p++){
-	  ep->Extrude(j, k, x[p], y[p], z[p]);
-	  ep->Extrude(j, k + 1, x[p + 4], y[p + 4], z[p + 4]);
-	}
-	for(int p = 0; p < 8; p++){
-	  MVertex tmp(x[p], y[p], z[p], 0, -1);
-	  itp = pos.find(&tmp);
-	  if(itp == pos.end()) {
-	    Msg(GERROR, "Could not find extruded vertex in volume %d", to->tag());
-	    return;
-	  }
-	  verts.push_back(*itp);
-	}
-	createHexPri(verts, to);
+	if(getExtrudedVertices(from->quadrangles[i], ep, j, k, pos, verts) == 8)
+	  createHexPri(verts, to);
       }
     }
   }
 }
 
-int MeshExtrudedVolume(GRegion *gr)
+void insertAllVertices(GRegion *gr, 
+		       std::set<MVertex*, MVertexLessThanLexicographic> &pos)
 {
-  ExtrudeParams *ep = gr->meshAttributes.extrude;
-
-  if(!ep || !ep->mesh.ExtrudeMesh || ep->geo.Mode != EXTRUDED_ENTITY)
-    return 0;
-
-  // FIXME
-  static bool warn_recombine = 1;
-  if(!ep->mesh.Recombine && warn_recombine){
-    Msg(WARNING, "Non-recombined volume extrusion has not been reimplemented yet");
-    warn_recombine = false;
-  }
-
-
-  // build a set with all the vertices on the boundary of gr
-  double old_tol = MVertexLessThanLexicographic::tolerance; 
-  MVertexLessThanLexicographic::tolerance = 1.e-12 * CTX.lc;
-  std::set<MVertex*, MVertexLessThanLexicographic> pos;
+  pos.insert(gr->mesh_vertices.begin(), gr->mesh_vertices.end());
   std::list<GFace*> faces = gr->faces();
-  std::list<GFace*>::iterator it = faces.begin();
-  while(it != faces.end()){
-    pos.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
-    std::list<GEdge*> edges = (*it)->edges();
+  std::list<GFace*>::iterator itf = faces.begin();
+  while(itf != faces.end()){
+    pos.insert((*itf)->mesh_vertices.begin(), (*itf)->mesh_vertices.end());
+    std::list<GEdge*> edges = (*itf)->edges();
     std::list<GEdge*>::iterator ite = edges.begin();
     while(ite != edges.end()){
       pos.insert((*ite)->mesh_vertices.begin(), (*ite)->mesh_vertices.end());
@@ -201,16 +181,297 @@ int MeshExtrudedVolume(GRegion *gr)
 		 (*ite)->getEndVertex()->mesh_vertices.end());
       ++ite;
     }
-    ++it;
+    ++itf;
   }
+}
+
+void meshGRegionExtruded::operator() (GRegion *gr) 
+{  
+  if(gr->geomType() == GEntity::DiscreteVolume) return;
+
+  ExtrudeParams *ep = gr->meshAttributes.extrude;
+
+  if(!ep || !ep->mesh.ExtrudeMesh || ep->geo.Mode != EXTRUDED_ENTITY) return;
+
+  // Send a messsage to the GMSH environment
+  Msg(STATUS2, "Meshing volume %d (extruded)", gr->tag());
+
+  // destroy the mesh if it exists
+  deMeshGRegion dem;
+  dem(gr);
+
+  // build a set with all the vertices on the boundary of gr
+  double old_tol = MVertexLessThanLexicographic::tolerance; 
+  MVertexLessThanLexicographic::tolerance = 1.e-12 * CTX.lc;
+  std::set<MVertex*, MVertexLessThanLexicographic> pos;
+  insertAllVertices(gr, pos);
 
   // volume is extruded from a surface
   GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
-  if(!from) return 0;
+  if(!from){
+    Msg(GERROR, "Unknown source surface %d for extrusion", ep->geo.Source);
+    return;
+  }
 
   extrudeMesh(from, gr, pos);
 
   MVertexLessThanLexicographic::tolerance = old_tol; 
+}
+
+int edgeExists(MVertex *v1, MVertex *v2, 
+	      std::set<std::pair<MVertex*, MVertex*> > &edges)
+{
+  std::pair<MVertex*, MVertex*> p(std::min(v1, v2), std::max(v1, v2));
+  return edges.count(p);
+}
+
+void createEdge(MVertex *v1, MVertex *v2, 
+	      std::set<std::pair<MVertex*, MVertex*> > &edges)
+{
+  std::pair<MVertex*, MVertex*> p(std::min(v1, v2), std::max(v1, v2));
+  edges.insert(p);
+}
+
+void deleteEdge(MVertex *v1, MVertex *v2, 
+	      std::set<std::pair<MVertex*, MVertex*> > &edges)
+{
+  std::pair<MVertex*, MVertex*> p(std::min(v1, v2), std::max(v1, v2));
+  edges.erase(p);
+}
+
+// subdivide the 3 lateral faces of each prism
+void phase1(GRegion *gr,
+	    std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+	    std::set<std::pair<MVertex*, MVertex*> > &edges)
+{
+  ExtrudeParams *ep = gr->meshAttributes.extrude;
+  GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
+  if(!from) return;
+
+  for(unsigned int i = 0; i < from->triangles.size(); i++){
+    for(int j = 0; j < ep->mesh.NbLayer; j++) {
+      for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
+	std::vector<MVertex*> v;
+	if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, v) == 6){
+	  if(!edgeExists(v[0], v[4], edges))
+	    createEdge(v[1], v[3], edges);
+	  if(!edgeExists(v[4], v[2], edges))
+	    createEdge(v[1], v[5], edges);
+	  if(!edgeExists(v[3], v[2], edges))
+	    createEdge(v[0], v[5], edges);
+	}
+      }
+    }
+  }
+}
+
+// modify lateral edges to make them "tet-compatible"
+void phase2(GRegion *gr,
+	    std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+	    std::set<std::pair<MVertex*, MVertex*> > &edges,
+	    std::set<std::pair<MVertex*, MVertex*> > &edges_swap,
+	    int &swap)
+{
+  ExtrudeParams *ep = gr->meshAttributes.extrude;
+  GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
+  if(!from) return;
+
+  for(unsigned int i = 0; i < from->triangles.size(); i++){
+    for(int j = 0; j < ep->mesh.NbLayer; j++) {
+      for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
+	std::vector<MVertex*> v;
+	if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, v) == 6){
+	  if(edgeExists(v[3], v[1], edges) &&
+	     edgeExists(v[4], v[2], edges) && 
+	     edgeExists(v[0], v[5], edges)) {
+	    swap++;
+	    if(!edgeExists(v[3], v[1], edges_swap)) {
+	      deleteEdge(v[3], v[1], edges);
+	      createEdge(v[0], v[4], edges);
+	      createEdge(v[0], v[4], edges_swap);
+	      createEdge(v[3], v[1], edges_swap);
+	    }
+	    else if(!edgeExists(v[4], v[2], edges_swap)) {
+	      deleteEdge(v[4], v[2], edges);
+	      createEdge(v[1], v[5], edges);
+	      createEdge(v[4], v[2], edges_swap);
+	      createEdge(v[1], v[5], edges_swap);
+	    }
+	    else if(!edgeExists(v[0], v[5], edges_swap)) {
+	      deleteEdge(v[0], v[5], edges);
+	      createEdge(v[3], v[2], edges);
+	      createEdge(v[0], v[5], edges_swap);
+	      createEdge(v[3], v[2], edges_swap);
+	    }
+	  }
+	  else if(edgeExists(v[0], v[4], edges) &&
+		  edgeExists(v[1], v[5], edges) &&
+		  edgeExists(v[3], v[2], edges)) {
+	    swap++;
+	    if(!edgeExists(v[0], v[4], edges_swap)) {
+	      deleteEdge(v[0], v[4], edges);
+	      createEdge(v[3], v[1], edges);
+	      createEdge(v[0], v[4], edges_swap);
+	      createEdge(v[3], v[1], edges_swap);
+	    }
+	    else if(!edgeExists(v[1], v[5], edges_swap)) {
+	      deleteEdge(v[1], v[5], edges);
+	      createEdge(v[4], v[2], edges);
+	      createEdge(v[4], v[2], edges_swap);
+	      createEdge(v[1], v[5], edges_swap);
+	    }
+	    else if(!edgeExists(v[3], v[2], edges_swap)) {
+	      deleteEdge(v[3], v[2], edges);
+	      createEdge(v[0], v[5], edges);
+	      createEdge(v[0], v[5], edges_swap);
+	      createEdge(v[3], v[2], edges_swap);
+	    }
+	  }
+	}
+      }
+    }
+  }
+}
+ 
+// create tets
+void phase3(GRegion *gr,
+	    std::set<MVertex*, MVertexLessThanLexicographic> &pos,
+	    std::set<std::pair<MVertex*, MVertex*> > &edges)
+{
+  ExtrudeParams *ep = gr->meshAttributes.extrude;
+  GFace *from = gr->model()->faceByTag(std::abs(ep->geo.Source));
+  if(!from) return;
+
+  for(unsigned int i = 0; i < from->triangles.size(); i++){
+    for(int j = 0; j < ep->mesh.NbLayer; j++) {
+      for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
+	std::vector<MVertex*> v;
+	if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, v) == 6){
+	  if(edgeExists(v[3], v[1], edges) &&
+	     edgeExists(v[4], v[2], edges) &&
+	     edgeExists(v[3], v[2], edges)) {
+	    createTet(v[0], v[1], v[2], v[3], gr);
+	    createTet(v[3], v[4], v[5], v[2], gr);
+	    createTet(v[1], v[3], v[4], v[2], gr);
+	  }
+	  else if(edgeExists(v[3], v[1], edges) &&
+		  edgeExists(v[1], v[5], edges) &&
+		  edgeExists(v[3], v[2], edges)) {
+	    createTet(v[0], v[1], v[2], v[3], gr);
+	    createTet(v[3], v[4], v[5], v[1], gr);
+	    createTet(v[3], v[1], v[5], v[2], gr);
+	  }
+	  else if(edgeExists(v[3], v[1], edges) &&
+		  edgeExists(v[1], v[5], edges) &&
+		  edgeExists(v[5], v[0], edges)) {
+	    createTet(v[0], v[1], v[2], v[5], gr);
+	    createTet(v[3], v[4], v[5], v[1], gr);
+	    createTet(v[1], v[3], v[5], v[0], gr);
+	  }
+	  else if(edgeExists(v[4], v[0], edges) &&
+		  edgeExists(v[4], v[2], edges) &&
+		  edgeExists(v[3], v[2], edges)) {
+	    createTet(v[0], v[1], v[2], v[4], gr);
+	    createTet(v[3], v[4], v[5], v[2], gr);
+	    createTet(v[0], v[3], v[4], v[2], gr);
+	  }
+	  else if(edgeExists(v[4], v[0], edges) &&
+		  edgeExists(v[4], v[2], edges) &&
+		  edgeExists(v[5], v[0], edges)) {
+	    createTet(v[0], v[1], v[2], v[4], gr);
+	    createTet(v[3], v[4], v[5], v[0], gr);
+	    createTet(v[0], v[2], v[4], v[5], gr);
+	  }
+	  else if(edgeExists(v[4], v[0], edges) &&
+		  edgeExists(v[1], v[5], edges) &&
+		  edgeExists(v[5], v[0], edges)) {
+	    createTet(v[0], v[1], v[2], v[5], gr);
+	    createTet(v[3], v[4], v[5], v[0], gr);
+	    createTet(v[0], v[1], v[4], v[5], gr);
+	  }
+	}
+      }
+    }
+  }
+}
+
+int SubdivideExtrudedMesh(GModel *m)
+{
+  // get all non-recombined extruded regions and vertices
+  std::vector<GRegion*> regions;
+  double old_tol = MVertexLessThanLexicographic::tolerance; 
+  MVertexLessThanLexicographic::tolerance = 1.e-12 * CTX.lc;
+  std::set<MVertex*, MVertexLessThanLexicographic> pos;
+  for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++){
+    ExtrudeParams *ep = (*it)->meshAttributes.extrude;
+    if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY && 
+       !ep->mesh.Recombine){
+      regions.push_back(*it);
+      insertAllVertices(*it, pos);
+    }
+  }
+
+  if(regions.empty()) return 0;
+
+  Msg(INFO, "Subdividing extruded mesh");
+
+  // create edges on lateral sides of "prisms"
+  std::set<std::pair<MVertex*, MVertex*> > edges;
+  for(unsigned int i = 0; i < regions.size(); i++)
+    phase1(regions[i], pos, edges);
+
+  // swap lateral edges to make them "tet-compatible"
+  int j = 0, swap;
+  std::set<std::pair<MVertex*, MVertex*> > edges_swap;
+  do {
+    swap = 0;
+    for(unsigned int i = 0; i < regions.size(); i++)
+      phase2(regions[i], pos, edges, edges_swap, swap);
+    Msg(INFO, "Swapping %d", swap);
+    if(j && j == swap) {
+      Msg(GERROR, "Unable to subdivide extruded mesh: use 'Recombine' instead");
+      return -1;
+    }
+    j = swap;
+  } while(swap);
+
+  // delete "recombined" volume elements and create tetrahedra instead
+  for(unsigned int i = 0; i < regions.size(); i++){
+    GRegion *gr = regions[i];
+    for(unsigned int i = 0; i < gr->hexahedra.size(); i++) 
+      delete gr->hexahedra[i];
+    gr->hexahedra.clear();
+    for(unsigned int i = 0; i < gr->prisms.size(); i++) 
+      delete gr->prisms[i];
+    gr->prisms.clear();
+    for(unsigned int i = 0; i < gr->pyramids.size(); i++)
+      delete gr->pyramids[i];
+    gr->pyramids.clear();
+    phase3(gr, pos, edges);
+  }
+
+  // re-Extrude bounding surfaces using edges as constraint
+  std::set<GFace*> faces;
+  for(unsigned int i = 0; i < regions.size(); i++){
+    std::list<GFace*> f = regions[i]->faces();
+    faces.insert(f.begin(), f.end());
+  }
+  for(std::set<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
+    ExtrudeParams *ep = (*it)->meshAttributes.extrude;
+    if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY && 
+       !ep->mesh.Recombine){
+      GFace *gf = *it;
+      Msg(INFO, "Remeshing surface %d", gf->tag());
+      for(unsigned int i = 0; i < gf->triangles.size(); i++) 
+	delete gf->triangles[i];
+      gf->triangles.clear();
+      for(unsigned int i = 0; i < gf->quadrangles.size(); i++) 
+	delete gf->quadrangles[i];
+      gf->quadrangles.clear();
+      MeshExtrudedSurface(gf, &edges);
+    }
+  }
 
-  return 1;
+  MVertexLessThanLexicographic::tolerance = old_tol;   
+  return 0;
 }
diff --git a/Mesh/meshGRegionTransfinite.cpp b/Mesh/meshGRegionTransfinite.cpp
index ad613ca8f4d6359257a4b8a1e1fa8a3b5bd0367e..fb4cc6cdc93712350902b7fdf7aefc59b95b0bda 100644
--- a/Mesh/meshGRegionTransfinite.cpp
+++ b/Mesh/meshGRegionTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionTransfinite.cpp,v 1.3 2006-12-02 21:16:49 geuzaine Exp $
+// $Id: meshGRegionTransfinite.cpp,v 1.4 2007-01-16 11:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -162,8 +162,8 @@ MVertex *transfiniteHex(MVertex *f1, MVertex *f2, MVertex *f3, MVertex *f4,
 class GOrientedTransfiniteFace {
 private:
   GFace *_gf;
-  int _permutation, _index;
   int _L, _H;
+  int _permutation, _index;
   std::vector<MVertex*> _list;
 
 public:
@@ -269,8 +269,8 @@ public:
     case 7: index = (m + M * n); break;
     }
     MVertex *v = 0;
-    if(index >= 0 && index < _list.size()) v = _list[index];
-    if(index < 0 || index >= _list.size() || !v){
+    if(index >= 0 && index < (int)_list.size()) v = _list[index];
+    if(index < 0 || index >= (int)_list.size() || !v){
       Msg(GERROR, "Wrong index in transfinite mesh of surface %d: "
 	  "m=%d n=%d M=%d N=%d perm=%d", _gf->tag(), m, n, M, N, _permutation);
       return _list[0];
diff --git a/Plugin/Plugin.cpp b/Plugin/Plugin.cpp
index 168912b007da1d75a15af9ac7ea0d233ed2298e9..f96fe17712ced6127121762f5c45449377b15aeb 100644
--- a/Plugin/Plugin.cpp
+++ b/Plugin/Plugin.cpp
@@ -1,4 +1,4 @@
-// $Id: Plugin.cpp,v 1.85 2006-11-27 22:22:32 geuzaine Exp $
+// $Id: Plugin.cpp,v 1.86 2007-01-16 11:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -185,8 +185,10 @@ void GMSH_PluginManager::registerDefaultPlugins()
 		      ("Extract", GMSH_RegisterExtractPlugin()));
     allPlugins.insert(std::pair < char *, GMSH_Plugin * >
 		      ("ExtractElements", GMSH_RegisterExtractElementsPlugin()));
+#if 0 // waiting for BDS rewrite
     allPlugins.insert(std::pair < char *, GMSH_Plugin * >
 		      ("ExtractEdges", GMSH_RegisterExtractEdgesPlugin()));
+#endif
     allPlugins.insert(std::pair < char *, GMSH_Plugin * >
 		      ("DecomposeInSimplex", GMSH_RegisterDecomposeInSimplexPlugin()));
     allPlugins.insert(std::pair < char *, GMSH_Plugin * >
diff --git a/Plugin/Triangulate.cpp b/Plugin/Triangulate.cpp
index dc9057aee424c8e269c7c789ba011d13eae3cd38..2c49d0f5f8325a497e5ae7c7093ead82a5911ff5 100644
--- a/Plugin/Triangulate.cpp
+++ b/Plugin/Triangulate.cpp
@@ -1,4 +1,4 @@
-// $Id: Triangulate.cpp,v 1.33 2006-11-27 22:22:32 geuzaine Exp $
+// $Id: Triangulate.cpp,v 1.34 2007-01-16 11:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -151,13 +151,13 @@ static void Triangulate(int nbIn, List_T *inList, int *nbOut, List_T *outList,
   in.holelist = NULL;
 
   j = 0;
-  for(int i = 0; i < points.size(); i++) {
+  for(unsigned int i = 0; i < points.size(); i++) {
     in.pointlist[j] = points[i]->x();
     in.pointlist[j + 1] = points[i]->y();
     j += 2;
   }
 
-  for(int i = 0; i < points.size(); i++) 
+  for(unsigned int i = 0; i < points.size(); i++) 
     delete points[i];
 
   struct triangulateio out;