diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index e975a888211f0dba549e0e6b5a2b6f1d03f49063..ebc66decf3c11c7c1d036d49b6bc95c2315d6e3c 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -35,6 +35,7 @@ class GRegion : public GEntity {
 
   virtual int dim() const {return 3;}
   virtual void setVisibility(char val, bool recursive=false);
+  virtual std::list<GFace*> faces() const{return l_faces;}
 
   // The bounding box
   virtual SBoundingBox3d bounds() const; 
diff --git a/Geo/MElement.h b/Geo/MElement.h
index dd146319d2020dc93012f11cf76138fb7f47dec7..0645fd93a9d4a1ae2bfcb90f1b196b35e9de0b0c 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -226,6 +226,12 @@ class MTriangle : public MElement {
   virtual int getDim(){ return 2; }
   inline int getNumVertices(){ return 3; }
   inline MVertex *getVertex(int num){ return _v[num]; }
+  inline void revert () 
+    {
+      MVertex *vv = _v[0];
+      _v[0] = _v[1];
+      _v[1] = vv;
+    }
   virtual int getNumEdges(){ return 3; }
   virtual MEdge getEdge(int num)
   {
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 6896150dcab317ee479c72ae2bf2c62040051cf1..69f8c0266100d9c3d16a802f301c9f46d3c00798 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.60 2006-08-29 10:39:48 remacle Exp $
+// $Id: BDS.cpp,v 1.61 2006-09-01 10:10:05 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -29,19 +29,6 @@
 
 bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS_Face *t);
 
-double BDS_Point::min_edge_length()
-{
-  std::list<BDS_Edge*>::iterator it  = edges.begin();
-  std::list<BDS_Edge*>::iterator ite = edges.end();
-  double L = 1.e245;
-  while(it!=ite){
-    double l = (*it)->length();
-    if (l<L)L=l;
-    ++it;
-  }
-  return L;
-}
-
 void outputScalarField(std::list < BDS_Face * >t, const char *iii)
 {
   FILE *f = fopen(iii, "w");
@@ -506,11 +493,12 @@ BDS_Face *BDS_Mesh::add_quadrangle(BDS_Edge * e1, BDS_Edge * e2,
 //   return t;
 // }
 
-void BDS_Mesh::del_triangle(BDS_Face * t)
+void BDS_Mesh::del_face(BDS_Face * t)
 {
   t->e1->del(t);
   t->e2->del(t);
   t->e3->del(t);
+  if(t->e4)t->e4->del(t);
   t->deleted = true;
 }
 
@@ -740,12 +728,12 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
 
   if(e->faces(0)) {
     g1 = e->faces(0)->g;
-    del_triangle(e->faces(0));
+    del_face(e->faces(0));
   }
   // not a bug !!!
   if(e->faces(0)) {
     g2 = e->faces(0)->g;
-    del_triangle(e->faces(0));
+    del_face(e->faces(0));
   }
 
   //  double t_l = e->target_length;
@@ -899,12 +887,12 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
 
   if(e->faces(0)) {
     g1 = e->faces(0)->g;
-    del_triangle(e->faces(0));
+    del_face(e->faces(0));
   }
   // not a bug !!!
   if(e->faces(0)) {
     g2 = e->faces(0)->g;
-    del_triangle(e->faces(0));
+    del_face(e->faces(0));
   }
   del_edge(e);
 
@@ -990,11 +978,11 @@ bool BDS_Mesh::recombine_edge(BDS_Edge * e)
   BDS_GeomEntity *g=0;
   if(e->faces(0)) {
     g = e->faces(0)->g;
-    del_triangle(e->faces(0));
+    del_face(e->faces(0));
   }
   // not a bug !!!
   if(e->faces(0)) {
-    del_triangle(e->faces(0));
+    del_face(e->faces(0));
   }
   del_edge(e);
 
@@ -1070,7 +1058,7 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
 
     while(it != ite) {
       BDS_Face *t = *it;
-      del_triangle(t);
+      del_face(t);
       ++it;
     }
   }
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 3ba289e297616828ea0b6f94efd5c70f082b8ca1..3048a9416e05f62dc7ab0d5471b8a3255fcb7810 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -38,6 +38,8 @@ class BDS_Mesh;
 class BDS_Point;
 class BDS_Vector;
 class GFace;
+class GEdge;
+class GVertex;
 
 void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]); 
 void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]); 
@@ -186,9 +188,7 @@ public:
       ++it;
     }
   }
-  double min_edge_length();
   void getTriangles(std::list<BDS_Face *> &t) const; 	
-  void compute_curvature();
   BDS_Point(int id, double x=0, double y=0, double z=0)
     : _lc(1.e22),X(x),Y(y),Z(z),u(0),v(0),config_modified(true),iD(id),g(0)
   {	    
@@ -386,6 +386,9 @@ public:
   int MAXPOINTNUMBER,SNAP_SUCCESS,SNAP_FAILURE;
   double Min[3],Max[3],LC;
   BDS_Mesh(int _MAXX = 0) :  MAXPOINTNUMBER(_MAXX){}
+  void load(GVertex   *gv); // load in BDS all the meshes of the vertex 
+  void load(GEdge     *ge); // load in BDS all the meshes of the edge 
+  void load(GFace     *gf); // load in BDS all the meshes of the surface 
   virtual ~BDS_Mesh();
   BDS_Mesh(const BDS_Mesh &other);
   std::set<BDS_GeomEntity*,GeomLessThan> geom; 
@@ -407,7 +410,7 @@ public:
   BDS_Face *add_quadrangle(int p1, int p2, int p3,int p4); 
   BDS_Face *add_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
   BDS_Face *add_quadrangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3,BDS_Edge *e4);
-  void del_triangle(BDS_Face *t);
+  void del_face(BDS_Face *t);
   BDS_Face  *find_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
   BDS_Face  *find_quadrangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3,BDS_Edge *e4);
   // Geom entities
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 550ca47179e2d63d72d84d410be1391a46fe9d68..c6ae1946dc30ad56ec4fe8152b624e03f209b939 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.94 2006-08-22 01:58:35 geuzaine Exp $
+// $Id: Generator.cpp,v 1.95 2006-09-01 10:10:05 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -31,6 +31,7 @@
 #include "OS.h"
 #include "meshGEdge.h"
 #include "meshGFace.h"
+#include "meshGRegion.h"
 #include "GModel.h"
 
 extern Mesh *THEM;
@@ -266,8 +267,6 @@ void Maillage_Dimension_1()
 
   double t1 = Cpu();
 
-  //  Tree_Action(THEM->Curves, Maillage_Curve);
-
   std::for_each(GMODEL->firstEdge(), GMODEL->lastEdge(), meshGEdge());
 
   double t2 = Cpu();
@@ -278,99 +277,23 @@ void Maillage_Dimension_2()
 {
   if(TooManyElements(2)) return;
 
-  double shortest = 1.e300;
-
   double t1 = Cpu();
-
-  // create reverse 1D meshes
-
-  List_T *Curves = Tree2List(THEM->Curves);
-  for(int i = 0; i < List_Nbr(Curves); i++) {
-    Curve *c;
-    List_Read(Curves, i, &c);
-    if(c->Num > 0) {
-      if(c->l < shortest)
-        shortest = c->l;
-      Curve C;
-      Curve *neew = &C;
-      neew->Num = -c->Num;
-      Tree_Query(THEM->Curves, &neew);
-      neew->Vertices =
-        List_Create(List_Nbr(c->Vertices), 1, sizeof(Vertex *));
-      List_Invert(c->Vertices, neew->Vertices);
-    }
-  }
-  List_Delete(Curves);
-
-  Msg(DEBUG, "Shortest curve has length %g", shortest);
-
-  // mesh 2D
-
-  //  Tree_Action(THEM->Surfaces, Maillage_Surface);
-
   std::for_each(GMODEL->firstFace(), GMODEL->lastFace(), meshGFace());
 
-  // global "all-quad" recombine
-
-  if(CTX.mesh.algo_recombine == 2)
-    Recombine_All(THEM);
+  // 2 BE DONE
+  //  if(CTX.mesh.algo_recombine == 2)
+  //    Recombine_All(THEM);
 
   double t2 = Cpu();
   CTX.mesh_timer[1] = t2 - t1;
 }
 
-static Volume *IVOL;
-
-void TransferData(void *a, void *b)
-{
-  Simplex *s = *(Simplex**)a;
-  if(s->iEnt == IVOL->Num){
-    Tree_Add(IVOL->Simplexes, &s);
-    for(int i = 0; i < 4; i++)
-      Tree_Insert(IVOL->Vertices, &s->V[i]);
-  }
-}
-
 void Maillage_Dimension_3()
 {
   if(TooManyElements(3)) return;
 
   double t1 = Cpu();
-
-  // merge all the delaunay parts in a single special volume
-  Volume *v = Create_Volume(99999, 99999);
-  List_T *list = Tree2List(THEM->Volumes);
-  for(int i = 0; i < List_Nbr(list); i++) {
-    Volume *vol;
-    List_Read(list, i, &vol);
-    if((!vol->Extrude || !vol->Extrude->mesh.ExtrudeMesh) &&
-       (vol->Method != TRANSFINI)) {
-      for(int j = 0; j < List_Nbr(vol->Surfaces); j++) {
-        List_Replace(v->Surfaces, List_Pointer(vol->Surfaces, j),
-                     compareSurface);
-      }
-    }
-  }
-  Tree_Insert(THEM->Volumes, &v);
-
-  if(CTX.mesh.oldxtrude) {
-    Extrude_Mesh_Old(); // old extrusion
-  }
-  else {
-    Extrude_Mesh(THEM->Volumes); // new extrusion
-    Tree_Action(THEM->Volumes, Maillage_Volume); // delaunay of remaining parts
-  }
-
-  // transfer data back to individual volumes and remove special volume
-  for(int i = 0; i < List_Nbr(list); i++){
-    List_Read(list, i, &IVOL);
-    Tree_Action(v->Simplexes, TransferData);
-  }
-  Tree_Suppress(THEM->Volumes, &v);
-  Free_Volume_But_Not_Elements(&v, NULL);
-
-  List_Delete(list);
-
+  std::for_each(GMODEL->firstRegion(), GMODEL->lastRegion(), meshGRegion());
   double t2 = Cpu();
   CTX.mesh_timer[2] = t2 - t1;
 }
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 363d2a88237eab9292b524181a78cbb91db5cab6..f8a30e427ab4a179f9a1948bd5790b7c024e3e65 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.126 2006-08-26 17:16:55 geuzaine Exp $
+# $Id: Makefile,v 1.127 2006-09-01 10:10:05 remacle Exp $
 #
 # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 #
@@ -65,6 +65,7 @@ SRC = 1D_Mesh.cpp \
       Metric.cpp \
       meshGEdge.cpp \
       meshGFace.cpp \
+      meshGRegion.cpp \
       Nurbs.cpp \
       Interpolation.cpp \
       SecondOrder.cpp \
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 2812775edba75bd66e68b9ac1c98ebfe653e4b2c..d7481e218ce1bb5c1a893d2ef430904a1262875c 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.9 2006-08-29 10:39:49 remacle Exp $
+// $Id: meshGFace.cpp,v 1.10 2006-09-01 10:10:05 remacle Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -870,7 +870,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
       {
 	BDS_Face *t = *itt;
 	if (!t->g)
-	  m->del_triangle (t);
+	  m->del_face (t);
 	++itt;
       }
   }
@@ -991,6 +991,8 @@ void deMeshGFace :: operator() (GFace *gf)
   gf->mesh_vertices.clear();
   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();
 }
 
 
@@ -1018,15 +1020,8 @@ void meshGFace :: operator() (GFace *gf)
   // compute the mean plane, this is sometimes useful 
   gf->computeMeanPlane(points);
 
-  Msg(DEBUG1, "points were put in parametric coords ...");
-
-  try
-    {
-      gmsh2DMeshGenerator ( gf ) ;
-    }
-  catch (...)
-    {
-    }
+  // mesh the face
+  gmsh2DMeshGenerator ( gf ) ;
 
   Msg(DEBUG1, "type %d %d triangles generated, %d internal vertices",
       gf->geomType(),gf->triangles.size(),gf->mesh_vertices.size());
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index a776a61577d6d5237396d3ac47362598103eb890..a3c4191f045c2d2da3e6b1d6b312f287807e709e 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -20,10 +20,6 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
-#include <map>
-#include "SPoint2.h"
-
-class MVertex;
 class GFace;
 // Create the mesh of the face
 class meshGFace 
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..42ea7dc5f595b07a0160215f385c15065261a2f0
--- /dev/null
+++ b/Mesh/meshGRegion.cpp
@@ -0,0 +1,275 @@
+#include "meshGRegion.h"
+#include "GRegion.h"
+#include "GFace.h"
+#include "BDS.h"
+#include "Message.h"
+#include <vector>
+
+#if defined(HAVE_NETGEN)
+namespace nglib {
+#include "nglib.h"
+#include "nglib_addon.h"
+}
+using namespace nglib;
+
+Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<MVertex*> & numberedV)
+{
+  NgAddOn_Init();
+  Ng_Mesh *_ngmesh = Ng_NewMesh();
+
+  std::list<GFace*> faces = gr->faces();
+  std::list<GFace*>::iterator it = faces.begin();
+  std::set<MVertex*> allBoundingVertices;
+  
+  if (importVolumeMesh)
+    allBoundingVertices.insert ( gr->mesh_vertices.begin() , gr->mesh_vertices.end() );
+
+  while (it != faces.end())
+    {
+      GFace *gf = (*it);      
+      for (int i = 0; i< gf->triangles.size(); i++)
+	{
+	  MTriangle *t = gf->triangles[i];
+	  allBoundingVertices.insert (t->getVertex(0));
+	  allBoundingVertices.insert (t->getVertex(1));
+	  allBoundingVertices.insert (t->getVertex(2));
+	}
+      ++it;
+    }
+  std::set<MVertex*>::iterator itv =  allBoundingVertices.begin();
+  int I=1;
+  while (itv != allBoundingVertices.end())
+    {
+      double tmp[3];
+      tmp[0] = (*itv)->x();
+      tmp[1] = (*itv)->y();
+      tmp[2] = (*itv)->z();
+      (*itv)->setNum(I++);
+      numberedV.push_back(*itv);
+      Ng_AddPoint(_ngmesh, tmp);
+      ++itv;
+    }
+  
+  it = faces.begin();
+  while (it != faces.end())
+    {
+      GFace *gf = (*it);      
+      for (int i = 0; i< gf->triangles.size(); i++)
+	{
+	  MTriangle *t = gf->triangles[i];
+	  int tmp[3];
+	  tmp[0] = t->getVertex(0)->getNum();
+	  tmp[1] = t->getVertex(1)->getNum();
+	  tmp[2] = t->getVertex(2)->getNum();	  
+	  Ng_AddSurfaceElement(_ngmesh, NG_TRIG, tmp);
+	}
+      ++it;
+    }
+
+  if (importVolumeMesh)
+    {
+      for (int i = 0; i< gr->tetrahedra.size(); i++)
+	{
+	  MTetrahedron *t = gr->tetrahedra[i];
+	  int tmp[4];
+	  tmp[0] = t->getVertex(0)->getNum();
+	  tmp[1] = t->getVertex(1)->getNum();
+	  tmp[2] = t->getVertex(2)->getNum();	  
+	  tmp[3] = t->getVertex(2)->getNum();	  
+	  Ng_AddVolumeElement(_ngmesh, NG_TET, tmp);
+	}
+    }
+
+  return _ngmesh;
+}
+
+void TransferVolumeMesh(GRegion *gr, Ng_Mesh * _ngmesh, std::vector<MVertex*> &numberedV)
+{
+  // Gets total number of vertices of Netgen's mesh
+  int nbv = Ng_GetNP(_ngmesh);  
+  if(!nbv) return;
+
+  int nbpts = numberedV.size();
+
+  // Create new volume vertices
+  for(int i = nbpts; i < nbv; i++) 
+    {
+      double tmp[3];
+      Ng_GetPoint(_ngmesh, i+1, tmp);
+      MVertex *v = new MVertex (tmp[0],tmp[1],tmp[2]);
+      numberedV.push_back(v);
+      gr->mesh_vertices.push_back(v);
+    }
+  
+  // Get total number of simplices of Netgen's mesh
+  int nbe = Ng_GetNE(_ngmesh);
+  
+  // Create new volume simplices
+  for(int i = 0; i < nbe; i++) {
+    int tmp[4];
+    Ng_GetVolumeElement(_ngmesh, i+1, tmp);
+    MTetrahedron *t = new MTetrahedron  ( numberedV [tmp[0]-1],
+					  numberedV [tmp[1]-1],
+					  numberedV [tmp[2]-1],
+					  numberedV [tmp[3]-1]);
+    gr->tetrahedra.push_back(t);
+  }  
+}
+#endif
+
+void deMeshGRegion :: operator() (GRegion *gr) 
+{
+  for (unsigned int i=0;i<gr->mesh_vertices.size();i++) delete gr->mesh_vertices[i];
+  gr->mesh_vertices.clear();
+  for (unsigned int i=0;i<gr->tetrahedra.size();i++) delete gr->tetrahedra[i];
+  gr->tetrahedra.clear();
+  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();
+}
+
+int intersect_line_triangle ( double X[3],
+			      double Y[3],
+			      double Z[3] , 
+			      double P[3], 
+			      double N[3] )
+{
+  double mat[3][3], det;
+  double b[3], res[3];
+  const double eps_prec = 1.e-5;
+
+  mat[0][0] = X[1] - X[0];
+  mat[0][1] = X[2] - X[0];
+  mat[0][2] = N[0];
+
+  mat[1][0] = Y[1] - Y[0];
+  mat[1][1] = Y[2] - Y[0];
+  mat[1][2] = N[1];
+
+  mat[2][0] = Z[1] - Z[0];
+  mat[2][1] = Z[2] - Z[0];
+  mat[2][2] = N[2];
+
+  b[0] = P[0] - X[0];
+  b[1] = P[1] - Y[0];
+  b[2] = P[2] - Z[0];
+
+  if(!sys3x3_with_tol(mat, b, res, &det))       
+    return 0;
+
+  //  Msg(INFO,"going there %g %g %g",res[0],res[1],res[2]);
+
+  if(res[0] >=  eps_prec && res[0] <=  1.0 - eps_prec && 
+     res[1] >=  eps_prec && res[1] <=  1.0 - eps_prec && 
+     1-res[0]-res[1] >=  eps_prec && 1-res[0]-res[1] <=  1.0 - eps_prec )
+    {
+      // the line clearly intersects the triangle
+      return (res[2] > 0) ? 1 : 0;
+    }
+  else if (res[0] <  -eps_prec || res[0] >   1.0 + eps_prec || 
+	   res[1] <  -eps_prec || res[1] >   1.0 + eps_prec || 
+	   1-res[0]-res[1] <  -eps_prec || 1-res[0]-res[1] >  1.0 + eps_prec )
+    {
+      // the line clearly does NOT intersect the triangle
+      return 0;
+    }
+  else
+    {
+      // the intersection is not robust, try another triangle
+      return -10000;
+    }
+
+}
+
+
+void meshNormalsPointOutOfTheRegion (GRegion *gr) 
+{
+  std::list<GFace*> faces = gr->faces();
+  std::list<GFace*>::iterator it = faces.begin();
+  while (it != faces.end())
+    {
+      GFace *gf = (*it);      
+      int nb_intersect = 0;
+      for (int i = 0; i< gf->triangles.size(); i++)
+	{
+	  MTriangle *t = gf->triangles[i];
+	  double X[3] = {t->getVertex(0)->x(),t->getVertex(1)->x(),t->getVertex(2)->x()};
+	  double Y[3] = {t->getVertex(0)->y(),t->getVertex(1)->y(),t->getVertex(2)->y()};
+	  double Z[3] = {t->getVertex(0)->z(),t->getVertex(1)->z(),t->getVertex(2)->z()};
+	  double P[3] = {(X[0]+X[1]+X[2])/3.,(Y[0]+Y[1]+Y[2])/3.,(Z[0]+Z[1]+Z[2])/3.};
+	  double v1[3] = {X[0]-X[1],Y[0]-Y[1],Z[0]-Z[1]};
+	  double v2[3] = {X[2]-X[1],Y[2]-Y[1],Z[2]-Z[1]};
+	  double N[3] ;
+	  prodve ( v1, v2 , N );
+	  norme(N);
+	  std::list<GFace*>::iterator it_b = faces.begin();
+	  while (it_b != faces.end())
+	    {
+	      GFace *gf_b = (*it_b);
+	      for (int i_b = 0; i_b< gf_b->triangles.size(); i_b++)
+		{
+		  MTriangle *t_b = gf_b->triangles[i_b];
+		  if (t_b != t)
+		    {
+		      double X_b[3] = {t_b->getVertex(0)->x(),t_b->getVertex(1)->x(),t_b->getVertex(2)->x()};
+		      double Y_b[3] = {t_b->getVertex(0)->y(),t_b->getVertex(1)->y(),t_b->getVertex(2)->y()};
+		      double Z_b[3] = {t_b->getVertex(0)->z(),t_b->getVertex(1)->z(),t_b->getVertex(2)->z()};
+		      int inters = intersect_line_triangle ( X_b,Y_b,Z_b,P,N );
+		      nb_intersect += inters;
+		    }
+		}
+	      ++it_b;
+	    }
+	  Msg (INFO,"Region %d Face %d, %d intersect",gr->tag(),gf->tag(),nb_intersect);
+	  if (nb_intersect >= 0) break; // negative value means intersetcion is not "robust"
+	}
+
+
+      if (nb_intersect % 2 == 1) // odd nb of intersections: the normal points inside the region 
+	{
+	  for (int i = 0; i< gf->triangles.size(); i++)
+	    {
+	      gf->triangles[i]->revert();
+	    }
+	}
+      ++it;
+    }
+}
+
+
+
+void meshGRegion :: operator() (GRegion *gr) 
+{  
+  if(gr->geomType() == GEntity::DiscreteVolume) return;
+
+  // destroy the mesh if it exists
+  deMeshGRegion dem;
+  dem(gr);
+
+  // Send a messsage to the GMSH environment
+  Msg(STATUS2, "Meshing volume %d", gr->tag());
+
+  // orient the triangles of with respect to this region
+  meshNormalsPointOutOfTheRegion (gr); 
+  
+  if(CTX.mesh.algo3d == FRONTAL_NETGEN)
+    {
+#if !defined(HAVE_NETGEN)
+    Msg(GERROR, "Netgen is not compiled in this version of Gmsh");
+#else
+    std::vector<MVertex*> numberedV;
+    Ng_Mesh * _ngmesh = buildNetgenStructure (gr, 0, numberedV);
+    Ng_Meshing_Parameters mp;
+    mp.maxh = 1;
+    mp.fineness = 1;
+    mp.secondorder = 0;
+    NgAddOn_GenerateVolumeMesh(_ngmesh, &mp); // does not optimize
+    TransferVolumeMesh(gr, _ngmesh, numberedV);
+    Ng_DeleteMesh(_ngmesh);
+    Ng_Exit();
+#endif
+    }
+}
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
new file mode 100644
index 0000000000000000000000000000000000000000..b85ccd75e2166d035559ba1ea512d3947f191eba
--- /dev/null
+++ b/Mesh/meshGRegion.h
@@ -0,0 +1,41 @@
+#ifndef _MESH_GREGION_H_
+#define _MESH_GREGION_H_
+
+// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+class GRegion;
+class GModel;
+// Create the mesh of the region
+class meshGRegion 
+{
+public :
+  void operator () ( GRegion * );
+};
+
+void meshAllRegionsAtOnce ( GModel *gm );
+
+// destroy the mesh of the face
+class deMeshGRegion 
+{
+public :
+  void operator () ( GRegion * );
+};
+
+#endif