diff --git a/Mesh/2D_Mesh_Aniso.cpp b/Mesh/2D_Mesh_Aniso.cpp
index f6287785022101354b991b572300afbfcaaccd78..e30d4a5b28b9f65d91c37b3dc521ac512627c261 100644
--- a/Mesh/2D_Mesh_Aniso.cpp
+++ b/Mesh/2D_Mesh_Aniso.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Mesh_Aniso.cpp,v 1.20 2001-08-20 07:38:30 geuzaine Exp $
+// $Id: 2D_Mesh_Aniso.cpp,v 1.21 2001-08-28 22:15:15 geuzaine Exp $
 
 /*
    Jean-Francois Remacle
@@ -666,17 +666,16 @@ void Convex_Hull_Mesh_2D (List_T * Points, Surface * s){
       } 
     */
     if (!THES){
-      Msg(GERROR, "Vertex (%g,%g,%g) in no simplex",
-          THEV->Pos.X, THEV->Pos.Y, THEV->Pos.Z);
-      THEV->Pos.X += CTX.mesh.rand_factor * LC2D * rand()/RAND_MAX;
-      THEV->Pos.Y += CTX.mesh.rand_factor * LC2D * rand()/RAND_MAX;
-      THEV->Pos.Z += CTX.mesh.rand_factor * LC2D * rand()/RAND_MAX;
+      Msg(GERROR, "Vertex (%g,%g,%g) in no simplex", THEV->Pos.X, THEV->Pos.Y, THEV->Pos.Z);
+      THEV->Pos.X += CTX.mesh.rand_factor * LC2D * (1.-(double)rand()/(double)RAND_MAX);
+      THEV->Pos.Y += CTX.mesh.rand_factor * LC2D * (1.-(double)rand()/(double)RAND_MAX);
+      THEV->Pos.Z += CTX.mesh.rand_factor * LC2D * (1.-(double)rand()/(double)RAND_MAX);
       Tree_Action (s->Simplexes, Action_First_Simplexes_2D);
     }
-    bool  ca_marche = Bowyer_Watson_2D (s, THEV, THES, 1);
+    bool ca_marche = Bowyer_Watson_2D (s, THEV, THES, 1);
     while(!ca_marche){
-      double dx = CTX.mesh.rand_factor * LC2D * rand()/RAND_MAX;
-      double dy = CTX.mesh.rand_factor * LC2D * rand()/RAND_MAX;
+      double dx = CTX.mesh.rand_factor * LC2D * (1.-(double)rand()/(double)RAND_MAX);
+      double dy = CTX.mesh.rand_factor * LC2D * (1.-(double)rand()/(double)RAND_MAX);
       THEV->Pos.X += dx;
       THEV->Pos.Y += dy;
       ca_marche = Bowyer_Watson_2D (s, THEV, THES, 1);
diff --git a/Mesh/3D_Coherence.cpp b/Mesh/3D_Coherence.cpp
index 008cce077d147992ac9b1fad50681bd1d173714d..b660744757992b1050153ec9400f0020e4bec453 100644
--- a/Mesh/3D_Coherence.cpp
+++ b/Mesh/3D_Coherence.cpp
@@ -1,4 +1,4 @@
-// $Id: 3D_Coherence.cpp,v 1.18 2001-08-28 15:09:47 geuzaine Exp $
+// $Id: 3D_Coherence.cpp,v 1.19 2001-08-28 22:15:15 geuzaine Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -1235,19 +1235,19 @@ int Coherence (Volume * v, Mesh * m){
     /* il reste des intersections */
     
     if (List_Nbr (ListFaces) != Nh - 2){
-      /*
-        printf("Recherche des intersections\n");
-        printf("La face initiale comprend %d faces existantes\n",List_Nbr(ListFaces));
-        printf("La face est divisee en %d points\n",List_Nbr(NewPoints));
-      */
+      
+      Msg(INFO, "Intersections left");
+      Msg(INFO, "   Initial face contains %d existing faces", List_Nbr(ListFaces));
+      Msg(INFO, "   Face divided in %d points",List_Nbr(NewPoints));
+      
       TreexNewv = Tree_Create (sizeof (xNewv), compxNewv);
       TetAdd = Tree_Create (sizeof (Simplex *), compareSimplex);
       TetDel = Tree_Create (sizeof (Simplex *), compareSimplex);
       Tree_Action (v->Simplexes, Recover_Face);
-      /*
-        printf("La face est divisee en %d points %d %d \n",
-               List_Nbr(NewPoints),Tree_Nbr(TetAdd),Tree_Nbr(TetDel));
-      */
+
+      Msg(INFO, "The face is divided in %d points %d %d",
+	  List_Nbr(NewPoints),Tree_Nbr(TetAdd),Tree_Nbr(TetDel));
+
       Actual_Tree = v->Simplexes;
       Tree_Action (TetAdd, _Add);
       Tree_Action (TetDel, _Del);
@@ -1258,6 +1258,7 @@ int Coherence (Volume * v, Mesh * m){
     Np = List_Nbr (NewPoints);
 
     if (1 || List_Nbr (ListFaces) == 2 * (Np - 1) - Nh){
+      // on passe tjs ici pour le moment, meme si ca ne marche pas...
       
       Msg(INFO, "Recoverable face (%d <--> %d=2*(%d-1)-%d)",
           List_Nbr (ListFaces), 2 * (Np - 1) - Nh, Np, Nh);
@@ -1265,14 +1266,6 @@ int Coherence (Volume * v, Mesh * m){
       for (j = 0; j < List_Nbr (v->Surfaces); j++){
         List_Read (v->Surfaces, j, &s);
         if (Tree_Search (s->Simplexes, &simp)){
-	  /*
-	  if(List_Nbr(ListFaces)>2){
-	    printf("suppressed tri %d\n", simp->Num);
-	    printf(" (%g %g %g)\n", simp->V[0]->Pos.X, simp->V[0]->Pos.Y,simp->V[0]->Pos.Z);
-	    printf(" (%g %g %g)\n", simp->V[1]->Pos.X, simp->V[1]->Pos.Y,simp->V[1]->Pos.Z);
-	    printf(" (%g %g %g)\n", simp->V[2]->Pos.X, simp->V[2]->Pos.Y,simp->V[2]->Pos.Z);
-	  }
-	  */
           for (k = 0; k < List_Nbr (ListFaces); k++){
             List_Read (ListFaces, k, &Face);
             simp1 = Create_Simplex_MemeSens (simp, Face.V[0], Face.V[1], Face.V[2]);
@@ -1283,14 +1276,6 @@ int Coherence (Volume * v, Mesh * m){
             Tree_Replace (v->Vertices, &Face.V[0]);
             Tree_Replace (v->Vertices, &Face.V[1]);
             Tree_Replace (v->Vertices, &Face.V[2]);
-	    /*
-	    if(List_Nbr(ListFaces)>2){
-	      printf("replaced by tri %d\n", simp1->Num);
-	      printf(" (%g %g %g)\n", simp1->V[0]->Pos.X, simp1->V[0]->Pos.Y,simp1->V[0]->Pos.Z);
-	      printf(" (%g %g %g)\n", simp1->V[1]->Pos.X, simp1->V[1]->Pos.Y,simp1->V[1]->Pos.Z);
-	      printf(" (%g %g %g)\n", simp1->V[2]->Pos.X, simp1->V[2]->Pos.Y,simp1->V[2]->Pos.Z);
-	    }
-	    */
           }
           Tree_Suppress (s->Simplexes, &simp);
         }
diff --git a/Mesh/3D_Mesh.cpp b/Mesh/3D_Mesh.cpp
index 98d88808b5a68fa4f30a7d2ec178a045c7e1919d..6d024fbed511332e08d16040bf785ec6dedca605 100644
--- a/Mesh/3D_Mesh.cpp
+++ b/Mesh/3D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 3D_Mesh.cpp,v 1.28 2001-08-24 06:58:19 geuzaine Exp $
+// $Id: 3D_Mesh.cpp,v 1.29 2001-08-28 22:15:15 geuzaine Exp $
 
 /*
  
@@ -470,7 +470,6 @@ int recur_bowyer (Simplex * s){
 bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){
   int i;
   Simplex *s;
-  //  static int init = 1;
   double volumeold, volumenew;
 
   THEV = v;
@@ -486,13 +485,9 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){
 
   Tsd = Tree_Create (sizeof (Simplex *), compareSimplex);
   Sim_Sur_Le_Bord = Tree_Create (sizeof (Simplex *), compareSimplex);
-  //  if (init){
-    //    init = 0;
-    //  }
   List_Reset (Simplexes_Destroyed);
   List_Reset (Simplexes_New);
 
-
   if (Methode){
     Tree_Action (m->Simplexes, Fill_Sim_Des);
   }
@@ -550,7 +545,6 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){
       List_Read (Simplexes_Destroyed, i, &s);
       if (!Tree_Suppress (m->Simplexes, &s))
         Msg(GERROR, "Impossible to delete simplex");
-      // CORRECTION FROM Free(s) to that
       Free_Simplex (&s,0);
     }
     
@@ -565,11 +559,6 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){
   return true;
 }
 
-double rand_sign(){
-  double d = ((double)rand()/(double)RAND_MAX) ;
-  return (d < 0.5)?-1.0:1.0;
-}
-
 void Convex_Hull_Mesh (List_T * Points, Mesh * m){
   int i, j, N, n;
   int Nbr_OK = 0, Nbr_NOTOK = 0;
@@ -603,25 +592,19 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){
       Msg(STATUS1, "Vol=%g",volume); 
     }
     if (!THES){
-      Msg(WARNING, "Vertex (%g,%g,%g) in no simplex",
-          THEV->Pos.X,THEV->Pos.Y,THEV->Pos.Z); 
-      THEV->Pos.X += 10 * CTX.mesh.rand_factor * LC3D * (double)rand()/(double)RAND_MAX;
-      THEV->Pos.Y += 10 * CTX.mesh.rand_factor * LC3D * (double)rand()/(double)RAND_MAX;
-      THEV->Pos.Z += 10 * CTX.mesh.rand_factor * LC3D * (double)rand()/(double)RAND_MAX;
+      Msg(WARNING, "Vertex (%g,%g,%g) in no simplex", THEV->Pos.X,THEV->Pos.Y,THEV->Pos.Z); 
+      THEV->Pos.X += CTX.mesh.rand_factor * LC3D * (1.-(double)rand()/(double)RAND_MAX);
+      THEV->Pos.Y += CTX.mesh.rand_factor * LC3D * (1.-(double)rand()/(double)RAND_MAX);
+      THEV->Pos.Z += CTX.mesh.rand_factor * LC3D * (1.-(double)rand()/(double)RAND_MAX);
       Tree_Action (m->Simplexes, Action_First_Simplexes);
     }
-    bool  ca_marche = Bowyer_Watson (m, THEV, THES, 1);
+    bool ca_marche = Bowyer_Watson (m, THEV, THES, 1);
     int count = 0;
     while(!ca_marche){
-      //Msg(INFO, "Unable to add point %d (%g,%g,%g)",
-      //	  THEV->Num, THEV->Pos.X,THEV->Pos.Y,THEV->Pos.Z );
       count ++;
-      double dx = rand_sign() * 1000 * CTX.mesh.rand_factor * LC3D *
-	(double)rand()/(double)RAND_MAX;
-      double dy = rand_sign() * 1000 * CTX.mesh.rand_factor * LC3D *
-	(double)rand()/(double)RAND_MAX;
-      double dz = rand_sign() * 1000 * CTX.mesh.rand_factor * LC3D *
-	(double)rand()/(double)RAND_MAX;
+      double dx = CTX.mesh.rand_factor * LC3D * (1.-(double)rand()/(double)RAND_MAX);
+      double dy = CTX.mesh.rand_factor * LC3D * (1.-(double)rand()/(double)RAND_MAX);
+      double dz = CTX.mesh.rand_factor * LC3D * (1.-(double)rand()/(double)RAND_MAX);
       THEV->Pos.X += dx;
       THEV->Pos.Y += dy;
       THEV->Pos.Z += dz;
@@ -634,8 +617,7 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){
       if(count > 5){
         N++;
         List_Add(POINTS,&THEV);
-        Msg(WARNING, "Unable to add point %d (will do it later)",
-            THEV->Num);
+        Msg(WARNING, "Unable to add point %d (will try it later)", THEV->Num);
         break;
       }
     }
@@ -654,16 +636,7 @@ void suppress_simplex (void *data, void *dum){
   Simplex **pv;
 
   pv = (Simplex **) data;
-  if ((*pv)->iEnt == 0)
-    List_Add (Suppress, pv);
-
-  /*
-  else{
-    for(i=0;i<List_Nbr(TrsfVolNum);i++)
-      if((*pv)->iEnt == (*(int*)List_Pointer(TrsfVolNum,i))->Num)
-        List_Add(Suppress,pv);
-  }
-  */
+  if ((*pv)->iEnt == 0) List_Add (Suppress, pv);
 }
 
 void add_in_bgm (void *a, void *b){
diff --git a/benchmarks/1d/nurbs.geo b/benchmarks/1d/nurbs.geo
index 3a6771d4412f36da189b5f91ab827bc168256674..28010ddd043707b29eabb87313bb71a102cb852c 100644
--- a/benchmarks/1d/nurbs.geo
+++ b/benchmarks/1d/nurbs.geo
@@ -1,10 +1,10 @@
 
-lc = 1;
+lc = 0.01;
 
 Point(1) = {0,0,0,lc};
-Point(2) = {1,0,0,lc};
-Point(3) = {1,1,0,lc};
-Point(4) = {0,1,0,lc};
+Point(2) = {1,0,0,2*lc};
+Point(3) = {1,1,0,0.1*lc};
+Point(4) = {0,1,0,4*lc};
 
 // A knot represents a parameter interval (and can be 0).
 // There are (N + D + 1) knots
diff --git a/benchmarks/2d/Square-Attr4.geo b/benchmarks/2d/Square-Attr4.geo
index e15277cc7318b8a048b6929fcc5ee66e98b9ab63..8c5cf71d4c8eeb9bbd42b3248572668c5eb4546c 100644
--- a/benchmarks/2d/Square-Attr4.geo
+++ b/benchmarks/2d/Square-Attr4.geo
@@ -13,6 +13,6 @@ Line(4) = {4,3};
 Point(55) = {0.2,.5,0,lc};                 
 Line Loop(5) = {1,2,3,4};                
 Plane Surface(6) = {5};                
-Attractor Point {55} = {.01,.1,3.0};       
+Attractor Point {55} = {.01,.1,3.0};
 Mesh.Algorithm = 2;
 
diff --git a/benchmarks/bugs/p19-bug.geo b/benchmarks/bugs/p19-bug.geo
new file mode 100644
index 0000000000000000000000000000000000000000..928cbd540c7789f6eee997d54c819ce2d2626d2f
--- /dev/null
+++ b/benchmarks/bugs/p19-bug.geo
@@ -0,0 +1,122 @@
+
+Mesh.CharacteristicLengthFactor = 2;
+
+fact     = 0.4 ;
+rondelle = fact * 0.01;
+iris     = fact * 0.004;
+size     = fact * 0.01;
+
+larg = 86.36e-3 / 2.0 ;
+long = 45.0e-3 ;
+l    = 16.e-3 / 2.0 ;
+re   = 90.e-3 / 2.0 ;
+re2  = re + 3.e-3 ;
+ri   = 8.e-3 / 2.0 ;
+ll   = 60.0e-3 / 2.0;
+a    = larg - ll;
+c45  = 0.5^0.5  ;
+
+hg   = 21.59e-3 ;
+hcav = 42.5e-3 ;
+
+x2   = long;
+y1   = ri * c45;
+y2   = ri ;
+y3   = l ;
+ss1  = l / re ;
+xtemp = re * ( ( 1.0 - ss1 * ss1 ) ^ 0.5 ) ;
+ss2  = ( 1.0 - ( xtemp / re2 ) * ( xtemp / re2 ) ) ^ 0.5 ;
+y4   = ss2 * xtemp ;
+ss3  = ll / re2;
+xtemp2 = re2 *  ( ( 1.0 - ss3 * ss3 ) ^ 0.5 ) ;
+x1   = long - ( xtemp - xtemp2 ) ;
+x6   = x1 + xtemp ;
+x3   = x6 - re * c45 ;
+x4   = x6 - ri ;
+x5   = x6 - ri * c45 ;
+x7   = x6 + ri * c45 ;
+x8   = x6 + ri ;
+x9   = x6 + re * c45 ;
+x10  = x6 + re ;
+y5   = ll ;
+y6   = re * c45 ;
+y7   = re ;
+y8   = larg ;
+
+
+Point(1) = {0.0,0.0,0.0,size};
+Point(2) = {0.0,y8,0.0,size};
+Point(3) = {x2,y8,0.0,size};
+Point(4) = {x2,y5,0.0,size};
+Point(5) = {x1,y4,0.0,size};
+Point(6) = {x1,y3,0.0,iris};
+Point(7) = {x1,0.0,0.0,iris};
+Point(8) = {x3,y6,0.0,size};
+Point(9) = {x6,y7,0.0,size};
+Point(10) = {x9,y6,0.0,size};
+Point(11) = {x10,0.0,0.0,size};
+Point(12) = {x8,0.0,0.0,rondelle};
+Point(13) = {x6,0.0,0.0,rondelle};
+Point(14) = {x4,0.0,0.0,rondelle};
+Point(15) = {x7,y1,0.0,rondelle};
+Point(16) = {x6,y2,0.0,rondelle};
+Point(17) = {x5,y1,0.0,rondelle};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {5,6};
+Line(5) = {1,7};
+Line(6) = {7,14};
+Line(10) = {6,7};
+Line(11) = {14,13};
+Line(12) = {13,12};
+Circle(16) = {11,13,10};
+Circle(17) = {10,13,9};
+Circle(18) = {9,13,8};
+Circle(19) = {8,13,6};
+Circle(20) = {4,13,5};
+Line(21) = {11,12};
+Circle(22) = {17,13,14};
+Circle(23) = {16,13,17};
+Circle(24) = {15,13,16};
+Circle(25) = {12,13,15};
+
+Line Loop(26) = {-5,1,2,3,20,4,10};
+Plane Surface(27) = {26};
+
+Line Loop(28) = {-6,-10,-19,-18,-17,-16,21,25,24,23,22};
+Plane Surface(29) = {28};
+
+Line Loop(30) = {11,12,25,24,23,22};
+Plane Surface(31) = {30};
+
+Extrude Surface{27, {0,0,hg} };
+Coherence;
+
+Extrude Surface{29, {0,0,hg} };
+Coherence;
+
+Extrude Surface{31, {0,0,hg} };
+Coherence;
+
+Extrude Surface {105, {0,0,hcav-hg} };
+Coherence;
+
+Extrude Surface {126, {0,0,hcav-hg} };
+Coherence;
+
+//Characteristic Length {58,71} = 0.01; 
+
+Point(85) = {0.0,0.0,0.0,1.0};
+Surface Loop(158) = {67,27,43,47,51,55,59,63,68};
+Volume(159) = {158};
+Surface Loop(159) = {112,29,84,67,125,92,96,100,104,108,116,120,124};
+Volume(160) = {159};
+
+GO      = 1 ;
+CAV     = 2 ;
+DIS     = 3 ; 
+CLDSRC  = 4 ;
+CLD     = 5 ;
+
diff --git a/benchmarks/3d/old_extrude.geo b/benchmarks/extrude/old_extrude.geo
similarity index 100%
rename from benchmarks/3d/old_extrude.geo
rename to benchmarks/extrude/old_extrude.geo
diff --git a/benchmarks/3d/old_extrude.par b/benchmarks/extrude/old_extrude.par
similarity index 100%
rename from benchmarks/3d/old_extrude.par
rename to benchmarks/extrude/old_extrude.par