diff --git a/Geo/GModel.h b/Geo/GModel.h
index 4a3a54e67a1fd40375a6ba3bebb9b68fe8265264..5a78a4b939daa6e97ee40aa3d96685ed2768d14a 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -331,6 +331,7 @@ class GModel
   int writeOCCBREP(const std::string &name);
   int importOCCShape(const void *shape);
   int applyOCCMeshConstraints(const void *constraints);
+  void addOCCFillet(std::vector<int> & , double &);
   void addShape(std::string name, std::vector<double> &p, std::string op);
 
   // Gmsh mesh file format
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 27ccb2c16edb13559c12998595cfe26601b646df..84f77cfe803ffaaaced931983cbf22972134d117 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -545,15 +545,23 @@ void GModel::addShape(std::string name, std::vector<double> &p,
       _occ_internals->Box(SPoint3(p[0],p[1],p[2]),
                           SPoint3(p[3],p[4],p[5]),o);
     }
+    else if (name=="Fillet"){
+      std::vector<int> edges;
+      for (int i=0;i<p.size()-1;i++)edges.push_back((int)p[i]);
+      addOCCFillet(edges,p[p.size()-1]);
+    }
     else{
-      // we should that at the end, a test now !!
-      _occ_internals->buildLists();
-      _occ_internals->buildGModel(this);
+      return;
     }
+    destroy();
+    _occ_internals->buildLists();
+    _occ_internals->buildGModel(this);
   }
   catch(Standard_Failure &err){
     Msg::Error("%s", err.GetMessageString());
   }
+
+
 }
 
 TopoDS_Shape  GlueFaces (const TopoDS_Shape& theShape,
@@ -604,6 +612,51 @@ TopoDS_Shape  GlueFaces (const TopoDS_Shape& theShape,
 //   return aRes;
 }
 
+void OCC_Internals::Fillet(std::vector<TopoDS_Edge> & edgesToFillet,
+			   double Radius){
+
+  // create a tool for fillet
+  BRepFilletAPI_MakeFillet fill (shape);
+  for (int i=0;i<edgesToFillet.size();++i){
+    fill.Add(edgesToFillet[i]);
+  }
+  for (int i = 1; i <= fill.NbContours(); i++) {
+    fill.SetRadius(Radius, i, 1);
+  }
+  fill.Build();
+  if (!fill.IsDone()) {
+    Msg::Error("Fillet can't be computed on the given shape with the given radius");
+    return;
+  }
+  shape = fill.Shape();
+
+  if (shape.IsNull()) return;
+
+  // Check shape validity
+  BRepCheck_Analyzer ana (shape, false);
+  if (!ana.IsValid()) {
+    Msg::Error("Fillet algorithm have produced an invalid shape result");
+  }
+
+
+}
+
+
+void GModel::addOCCFillet(std::vector<int> & edgesToFillet,
+			  double &Radius)
+{
+  std::vector<TopoDS_Edge> toto;
+
+  for (int i=0;i<edgesToFillet.size();++i){
+    GEdge *ge = getEdgeByTag(edgesToFillet[i]);
+    if (ge->getNativeType() == GEntity::OpenCascadeModel){
+      toto.push_back(*(TopoDS_Edge*)ge->getNativePtr());
+    }  
+  }
+  _occ_internals->Fillet(toto,Radius);
+}
+
+
 void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool, const BooleanOperator &op)
 {
   if (tool.IsNull()) return;
@@ -822,7 +875,7 @@ void OCC_Internals::Torus(const SPoint3 &p, const SVector3 &d, double R1, double
   BRepPrimAPI_MakeTorus MC (anAxes, R1, R2);
   MC.Build();
   if (!MC.IsDone()) {
-    Msg::Error("Cylinder can't be computed from the given parameters");
+    Msg::Error("Torus can't be computed from the given parameters");
     return;
   }
   TopoDS_Shape aShape = MC.Shape();
@@ -839,7 +892,7 @@ void OCC_Internals::Torus(const SPoint3 &p, const SVector3 &d, double R1, double
   BRepPrimAPI_MakeTorus MC(anAxes, R1, R2, angle);
   MC.Build();
   if (!MC.IsDone()) {
-    Msg::Error("Cylinder can't be computed from the given parameters");
+    Msg::Error("Torus can't be computed from the given parameters");
     return;
   }
   TopoDS_Shape aShape = MC.Shape();
@@ -1145,4 +1198,9 @@ void GModel::addShape(std::string name, std::vector<double> &p,
              "Boolean Operators On Solids");
 }
 
+void GModel:: addOCCFillet(std::vector<int> & , double &){
+  Msg::Error("Gmsh must be compiled with OpenCascade support to apply "
+             "the Fillet operator");
+}
+
 #endif
diff --git a/Geo/GModelIO_OCC.h b/Geo/GModelIO_OCC.h
index c66c7863441d8e06ba13ce8583fdf1b0b6adb08a..563b7dea7de3878dcb954eebd1892a6f4eba9c7a 100644
--- a/Geo/GModelIO_OCC.h
+++ b/Geo/GModelIO_OCC.h
@@ -58,6 +58,7 @@ class OCC_Internals {
 	     const BooleanOperator &op);
   void Torus(const SPoint3 &bottom_center, const SVector3 &dir, double R1, double R2, 
 	     double angle,  const BooleanOperator &op);
+  void Fillet(std::vector<TopoDS_Edge> &shapes, double radius);
   void applyBooleanOperator(TopoDS_Shape tool, const BooleanOperator &op);
 };
 
diff --git a/Geo/OCCIncludes.h b/Geo/OCCIncludes.h
index 2430eae64420c18fd62018d913b6d7150d0bc4f1..067b126f7a2ad14b4e5a5306058b31f64245beac 100644
--- a/Geo/OCCIncludes.h
+++ b/Geo/OCCIncludes.h
@@ -114,6 +114,7 @@ using std::iostream;
 #include "BRepAlgoAPI_Cut.hxx"
 #include "BRepAlgoAPI_Section.hxx"
 #include "BRepAlgoAPI_Fuse.hxx"
+#include "BRepFilletAPI_MakeFillet.hxx"
 #endif
 
 #if defined(WIN32)
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 6f7c9f56cc79fadc83810deb05fd46ae6f0b1add..91031ed1dabbb5d20ccfedee7067e22cf5427c81 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -409,10 +409,23 @@ static void Mesh2D(GModel *m)
   // and curve meshes) is global as it depends on a smooth normal
   // field generated from the surface mesh of the source surfaces
   if(!Mesh2DWithBoundaryLayers(m)){
-    std::for_each(m->firstFace(), m->lastFace(), meshGFace());        
+
+    std::vector<GFace*> blob;
+    for(GModel::fiter it = m->firstFace() ; it!=m->lastFace(); ++it){
+      blob.push_back(*it);
+    }    
+
+
+    meshGFace mesher;
+#pragma omp parallel for schedule(dynamic) 
+    for(int i=0; i<blob.size();i++){
+      mesher(blob[i]);
+    }    
+#pragma omp barrier
+
+      //    std::for_each(m->firstFace(), m->lastFace(), meshGFace());        
     int nIter = 0;
     while(1){
-      meshGFace mesher;
       int nbPending = 0;
       for(GModel::fiter it = m->firstFace() ; it!=m->lastFace(); ++it){
         if ((*it)->meshStatistics.status == GFace::PENDING){
diff --git a/benchmarks/boolean/FILLET.geo b/benchmarks/boolean/FILLET.geo
new file mode 100644
index 0000000000000000000000000000000000000000..34b6ee31d8cc529830e09c917487167d59423687
--- /dev/null
+++ b/benchmarks/boolean/FILLET.geo
@@ -0,0 +1,14 @@
+L = 100;
+H = 30;
+Z = 10;
+OCCShape("Box",{0,0,0,L,H,Z},"none");
+R = 4;
+OCCShape("Fillet",{1:12,R},"none");
+OCCShape("Cone",{0*L/2,H/2,-Z,0,0,1,.3*R,2*R,3*Z},"Fuse");
+OCCShape("Fillet",{1,R},"none");
+OCCShape("Fillet",{14,R/8},"none");
+OCCShape("Cone",{0.99*L/2,H/2,-Z,0,0,1,.3*R,2*R,3*Z},"Fuse");
+OCCShape("Fillet",{1,R},"none");
+OCCShape("Fillet",{83,R/8},"none");
+OCCShape("Cone",{0.99*L,H/2,-Z,0,0,1,3*R,.2*R,3*Z},"Cut");
+Compound Surface(100000000) = {31, 40, 43, 13};
diff --git a/benchmarks/boolean/JAW.geo b/benchmarks/boolean/JAW.geo
new file mode 100644
index 0000000000000000000000000000000000000000..4f6e7befaf8ce509d4836908b15fa11c54a65bf7
--- /dev/null
+++ b/benchmarks/boolean/JAW.geo
@@ -0,0 +1,16 @@
+L = 100;
+H = 30;
+Z = 10;
+OCCShape("Box",{0,0,0,L,H,Z},"none");
+R = 10;
+X = 5;
+
+For I In {0:5}
+ OCCShape("Cylinder",{2*I*X,H/2,-3*Z,0,0,1,R,6*Z},"Cut");
+EndFor
+
+OCCShape("Sphere",{H-X,H/2,Z/2,R},"Fuse");
+OCCShape("Torus",{L,H/2,Z/2,0,0,1,2*R,R/2},"Fuse");
+
+
+OCCShape("End",{},"none");
diff --git a/benchmarks/boolean/TORUS.geo b/benchmarks/boolean/TORUS.geo
new file mode 100644
index 0000000000000000000000000000000000000000..6ec01cb0be6d0e2799597a1c4dfaf3043d99dfcf
--- /dev/null
+++ b/benchmarks/boolean/TORUS.geo
@@ -0,0 +1,2 @@
+OCCShape("Torus",{0,0,0,0,0,1,10,9.9},"None");
+OCCShape("End",{},"none");