diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 6bdde900fa97005b561b68a558d44fa27785711b..67d08c6b19729c2771a0665070030fa58d76171f 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2281,11 +2281,12 @@ GEntity *GModel::extrude(GEntity *e, std::vector<double> p1, std::vector<double>
   return 0;
 }
 
-GEntity *GModel::extrudeBoundaryLayer(GEntity *e, int nbLayers, double hLayers, int dir, int view)
+std::vector<GEntity*> GModel::extrudeBoundaryLayer(GEntity *e, int nbLayers, double hLayers, int dir, int view)
 {
   if(_factory)
     return _factory->extrudeBoundaryLayer(this, e, nbLayers,hLayers, dir, view);
-  return 0;
+  std::vector<GEntity*> empty;
+  return empty;
 }
 
 GEntity *GModel::addPipe(GEntity *e, std::vector<GEdge *>  edges)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 7a890dc271bbe0dcf7c54faddf98b7d8d831bb46..cc7b3c3374430fa2e7631e7bef2bc8b8e2b802f4 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -450,7 +450,7 @@ class GModel
   GEntity *revolve(GEntity *e, std::vector<double> p1, std::vector<double> p2,
                    double angle);
   GEntity *extrude(GEntity *e, std::vector<double> p1, std::vector<double> p2);
-  GEntity *extrudeBoundaryLayer(GEntity *e, int nbLayers, double hLayers, int dir=1, int view=-1);
+  std::vector<GEntity*> extrudeBoundaryLayer(GEntity *e, int nbLayers, double hLayers, int dir=1, int view=-1);
   GEntity *addPipe(GEntity *e, std::vector<GEdge *>  edges);
 
   void addRuledFaces(std::vector<std::vector<GEdge *> > edges);
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 9a7a5768db0cc871a862dcbe50031c169f232069..35020c6b682f564e4bbd2b73fa0ebd819fda762b 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -285,7 +285,7 @@ std::vector<GFace *> GeoFactory::addRuledFaces(GModel *gm,
   return faces;
 }
 
-GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayer, int dir, int view)
+std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayer, int dir, int view)
 {
 
   ExtrudeParams *ep = new  ExtrudeParams;
@@ -297,6 +297,8 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
   ep->mesh.NbElmLayer.clear(); 
   ep->mesh.NbElmLayer.push_back(nbLayers);
   ep->mesh.ExtrudeMesh = true;
+  if (CTX::instance()->mesh.recombineAll)  ep->mesh.Recombine = true;
+  else ep->mesh.Recombine = false;
   ep->geo.Source = e->tag();
  
   int type  = BOUNDARY_LAYER; 
@@ -322,18 +324,6 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
   else if (e->dim() == 2){
     ((GFace*)e)->meshAttributes.extrude = ep;
     Surface *s = FindSurface(e->tag());
-     // Msg::Info("Geo internal model has:");
-     // List_T *points = Tree2List(gm->getGEOInternals()->Points);
-     // List_T *curves = Tree2List(gm->getGEOInternals()->Curves);
-     // List_T *surfaces = Tree2List(gm->getGEOInternals()->Surfaces);
-     // Msg::Info("%d Vertices", List_Nbr(points));
-     // Msg::Info("%d Edges", List_Nbr(curves));
-     // Msg::Info("%d Faces", List_Nbr(surfaces));
-    // for(int i = 0; i < List_Nbr(surfaces); i++) {
-    //   Surface *s;
-    //   List_Read(surfaces, i, &s);
-    //   printf("surface %d \n", s->Num);
-    // }
     if(!s) {
       printf("surface %d NOT found \n", e->tag());
       exit(1);
@@ -357,33 +347,72 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
   gm->importGEOInternals();
 
   //return the new created entity
-  GEntity *newEnt=0;
-  if (e->dim() == 1){
-   for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++){
-     GEdge *ge = *it;
-    if(ge->getNativeType() == GEntity::GmshModel &&
-       ge->geomType() == GEntity::BoundaryLayerCurve){
-      ExtrudeParams *ep = ge->meshAttributes.extrude;
-      if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY &&
-	std::abs(ep->geo.Source) ==e->tag() )
-	  newEnt = ge;
-      }
-    }
+  int nbout = List_Nbr(list_out);
+  std::vector<GEntity*> extrudedEntities;
+  //GEntity *newEnt =0;
+  if(e->dim()==1){
+    Shape e;
+    Shape s;
+    List_Read(list_out, 0, &e);
+    List_Read(list_out, 1, &s);
+    GEdge *ge = gm->getEdgeByTag(e.Num);
+    GFace *gf = gm->getFaceByTag(s.Num);
+    extrudedEntities.push_back((GEntity*)ge);
+    extrudedEntities.push_back((GEntity*)gf);
+    for (int j=2; j<nbout; j++){
+      Shape el;
+      List_Read(list_out, j, &el);
+      GEdge *gel = gm->getEdgeByTag(el.Num);
+      extrudedEntities.push_back((GEntity*)gel);
+    }  
   }
-  else if (e->dim() ==2){
-    for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); it++){
-      GFace *gf = *it;
-      if(gf->getNativeType() == GEntity::GmshModel &&
-	 gf->geomType() == GEntity::BoundaryLayerSurface){
-	ExtrudeParams *ep = gf->meshAttributes.extrude;
-	if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY 
-	   && std::abs(ep->geo.Source) == e->tag())
-	  newEnt = gf;
-      }
-    }
+  else if(e->dim()==2){
+    Shape s;
+    Shape v;
+    List_Read(list_out, 0, &s);
+    List_Read(list_out, 1, &v);
+    GFace *gf = gm->getFaceByTag(s.Num);
+    GRegion *gr = gm->getRegionByTag(v.Num);
+    extrudedEntities.push_back((GEntity*)gf);
+    extrudedEntities.push_back((GEntity*)gr);
+    for (int j=2; j<nbout; j++){
+      Shape sl;
+      List_Read(list_out, j, &sl);
+      GFace *gfl = gm->getFaceByTag(sl.Num);
+      extrudedEntities.push_back((GEntity*)gfl);
+    }  
   }
 
-  return newEnt;
+  return extrudedEntities;
+
+  // //return the new created entity
+  // GEntity *newEnt=0;
+  // if (e->dim() == 1){
+  //  for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++){
+  //    GEdge *ge = *it;
+  //   if(ge->getNativeType() == GEntity::GmshModel &&
+  //      ge->geomType() == GEntity::BoundaryLayerCurve){
+  //     ExtrudeParams *ep = ge->meshAttributes.extrude;
+  //     if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY &&
+  // 	std::abs(ep->geo.Source) ==e->tag() )
+  // 	  newEnt = ge;
+  //     }
+  //   }
+  // }
+  // else if (e->dim() ==2){
+  //   for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); it++){
+  //     GFace *gf = *it;
+  //     if(gf->getNativeType() == GEntity::GmshModel &&
+  // 	 gf->geomType() == GEntity::BoundaryLayerSurface){
+  // 	ExtrudeParams *ep = gf->meshAttributes.extrude;
+  // 	if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY 
+  // 	   && std::abs(ep->geo.Source) == e->tag())
+  // 	  newEnt = gf;
+  //     }
+  //   }
+  // }
+
+  // return newEnt;
 
 };
 
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index 076c805c1b9cdf35ac949ae6995facd4d2ac9477..3379f622cf5c5bbc960d68e3f6db4ad7408801a2 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -90,10 +90,11 @@ class GModelFactory {
     Msg::Error("extrude not implemented yet");
     return 0;
   }
-  virtual GEntity* extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayers, int dir, int view)
+  virtual std::vector<GEntity*> extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayers, int dir, int view)
   {
     Msg::Error("extrude normals not implemented yet");
-    return 0;
+    std::vector<GEntity*> empty;
+    return empty;
   }
   virtual GEntity *addPipe(GModel *gm, GEntity *base, std::vector<GEdge *> wire)
   {
@@ -182,7 +183,7 @@ class GeoFactory : public GModelFactory {
   GRegion *addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces);
   GEdge *addCircleArc(GModel *gm,GVertex *begin, GVertex *center, GVertex *end);
   std::vector<GFace *> addRuledFaces(GModel *gm, std::vector<std::vector<GEdge *> > edges);
-  GEntity* extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayers, int dir, int view);
+  std::vector<GEntity*> extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayers, int dir, int view);
 };
 
 
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 8b77da198a7072b75e179851959813de19ad445c..82580b932f75840da9cf21d6abf30757e8648969 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2685,20 +2685,6 @@ void ExtrudeShapes(int type, List_T *list_in,
                    ExtrudeParams *e,
                    List_T *list_out)
 {
-
-   // Msg::Info("IN EXTRUDE SHAPE Geo internal model has:");
-   //  List_T *points = Tree2List(GModel::current()->getGEOInternals()->Points);
-   //  List_T *curves = Tree2List(GModel::current()->getGEOInternals()->Curves);
-   //  List_T *surfaces = Tree2List(GModel::current()->getGEOInternals()->Surfaces);
-   //  Msg::Info("%d Vertices", List_Nbr(points));
-   //  Msg::Info("%d Edges", List_Nbr(curves));
-   //  Msg::Info("%d Faces", List_Nbr(surfaces));
-   //  for(int i = 0; i < List_Nbr(surfaces); i++) {
-   //    Surface *s;
-   //    List_Read(surfaces, i, &s);
-   //    printf("surface %d \n", s->Num);
-   //  }
-
  
   for(int i = 0; i < List_Nbr(list_in); i++){
     Shape shape;
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 8eb07acf3eee17b7fcd146612638002e38bcc484..ebaef3a4d9bcbdc7a90be3bc8283a753b51e5903 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -349,10 +349,6 @@ Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
   is_closed = 0;
   is_extruded = 0;
 
-  // callbacks["closeVolume"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::closeVolume, "Create In/Outlet planar faces \n");
-  // callbacks["extrudeWall"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::extrudeWall, "Extrude wall \n");
-  // callbacks["cutMesh"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::cutMesh, "Cut the initial mesh in different mesh partitions using the centerlines \n");
-
   options["closeVolume"] = new FieldOptionInt(is_closed, "Action: Create In/Outlet planar faces");
   options["extrudeWall"] = new FieldOptionInt(is_extruded, "Action: Extrude wall");
   options["cutMesh"] = new FieldOptionInt(is_cut, "Action: Cut the initial mesh in different mesh partitions using the centerlines");
@@ -505,7 +501,8 @@ void Centerline::createBranches(int maxN){
       edges.push_back(newBranch) ;
     }
   }
-  printf("*** Centerline in/outlets =%d branches =%d \n", (int)color_edges.size()+1, (int)edges.size());
+
+  Msg::Info("Centerline: in/outlets =%d branches =%d ", (int)color_edges.size()+1, (int)edges.size());
 
   //create children
   for(unsigned int i = 0; i < edges.size(); ++i) {
@@ -628,7 +625,7 @@ void Centerline::createSplitCompounds(){
   NR = current->getMaxElementaryNumber(3);
 
   // Remesh new faces (Compound Lines and Compound Surfaces)
-  Msg::Info("Centerline: create split compounds:");
+  Msg::Info("Centerline: creates split compounds");
 
   //Parametrize Compound Lines
   for (int i=0; i < NE; i++){
@@ -639,8 +636,6 @@ void Centerline::createSplitCompounds(){
     Msg::Info("Create Compound Line (%d) = %d discrete edge",
               num_gec, pe->tag());
     GEdge *gec = current->addCompoundEdge(e_compound,num_gec);
-    //GEdgeCompound *gec = new GEdgeCompound(current, num_gec, e_compound);
-    //current->add(gec);
   }
 
   // Parametrize Compound surfaces
@@ -653,12 +648,7 @@ void Centerline::createSplitCompounds(){
     Msg::Info("Create Compound Surface (%d) = %d discrete face",
               num_gfc, pf->tag());
 
-    GFace *gfc = current->addCompoundFace(f_compound, 7, 0, num_gfc); //1=conf_spectral 4=convex_circle
-    ////GFaceCompound::typeOfCompound typ = GFaceCompound::CONVEX_CIRCLE; 
-    ////GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_PLANE; 
-    //GFaceCompound::typeOfCompound typ = GFaceCompound::CONFORMAL_SPECTRAL; 
-    //GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,typ, 0);
-    //current->add(gfc);
+    GFace *gfc = current->addCompoundFace(f_compound, 7, 0, num_gfc); //1=conf_spectral 4=convex_circle, 7=conf_fe
 
     gfc->meshAttributes.recombine = recombine;
     gfc->addPhysicalEntity(100);
@@ -771,7 +761,7 @@ void Centerline::createFaces(){
         it != touched.end(); ++it)
       e2e.erase(*it);
   }
-  Msg::Info("Centerline action (cutMesh) has cut surface mesh in %d faces ", (int)faces.size());
+  Msg::Info("Centerline: action (cutMesh) has cut surface mesh in %d faces ", (int)faces.size());
 
   //create discFaces
   for(unsigned int i = 0; i < faces.size(); ++i){
@@ -796,28 +786,7 @@ void Centerline::createFaces(){
 
 }
 
-void Centerline::createClosedVolume(){
-
-  //identify the boundary edges by looping over all discreteFaces
-  std::vector<GEdge*> boundEdges;
-  double dist_inlet = 1.e6;
-  GEdge *gin = NULL;
-  for (int i= 0; i< NF; i++){
-    GFace *gf = current->getFaceByTag(i+1);
-    std::list<GEdge*> l_edges = gf->edges();
-    for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
-      std::vector<GEdge*>::iterator ite = std::find(boundEdges.begin(), boundEdges.end(), *it);
-      if (ite != boundEdges.end()) boundEdges.erase(ite);
-      else boundEdges.push_back(*it);
-      GVertex *gv = (*it)->getBeginVertex();
-      SPoint3 pt(gv->x(), gv->y(), gv->z());
-      double dist = pt.distance(ptin);
-      if(dist < dist_inlet){
-	dist_inlet = dist;
-	gin = *it;
-      }			      
-    }
-  }
+void Centerline::createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges){
  
   current->setFactory("Gmsh");
   std::vector<std::vector<GFace *> > myFaceLoops;
@@ -842,24 +811,26 @@ void Centerline::createClosedVolume(){
     myFaces.push_back(newFace);
   }
 
- Msg::Info("Centerline action (closeVolume) has created %d in/out planar faces ", (int)boundEdges.size());
+  Msg::Info("Centerline: action (closeVolume) has created %d in/out planar faces ", (int)boundEdges.size());
 
   for (int i = 0; i < NF; i++){
-    GFace * gf = current->getFaceByTag(NF+i+1);
-     myFaces.push_back(gf);
+    GFace * gf;
+    if(is_cut) gf = current->getFaceByTag(NF+i+1);
+    else gf = current->getFaceByTag(i+1);
+    myFaces.push_back(gf);
   }
   myFaceLoops.push_back(myFaces);
   GRegion *reg = current->addVolume(myFaceLoops);
   reg->addPhysicalEntity(reg->tag());
   current->setPhysicalName("lumenVolume", 3, reg->tag());
 
-  Msg::Info("Centerline action (closeVolume) has created volume %d ", reg->tag());
+  Msg::Info("Centerline: action (closeVolume) has created volume %d ", reg->tag());
 
 }
 
-void Centerline::extrudeBoundaryLayerWall(){
+void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundEdges){
   
-  Msg::Info("Centerline: extrude boundary layer wall %d %g ", nbElemLayer,  hLayer);
+  Msg::Info("Centerline: extrude boundary layer wall (%d, %g%%R) ", nbElemLayer,  hLayer);
 
   //orient extrude direction outward
   int dir = 0;
@@ -873,13 +844,39 @@ void Centerline::extrudeBoundaryLayerWall(){
   if (dot(ne,nc) < 0) dir = 1;
   if (dir ==1 && hLayer > 0 ) hLayer *= -1.0;
 
+  int shift = 0;
+  if(is_cut) shift = NE;
   for (int i= 0; i< NF; i++){
     GFace *gfc ; 
-    if (is_cut) gfc = current->getFaceByTag(NF+i+1);//at this moment compound is not meshed yet exist yet
+    if (is_cut) gfc = current->getFaceByTag(NF+i+1);
     else gfc = current->getFaceByTag(i+1);
     current->setFactory("Gmsh");
-    current->extrudeBoundaryLayer(gfc, nbElemLayer,  hLayer, dir, -5);
     //view -5 to scale hLayer by radius in BoundaryLayers.cpp
+    std::vector<GEntity*> extrudedE = current->extrudeBoundaryLayer(gfc, nbElemLayer,  hLayer, dir, -5);
+    GFace *eFace = (GFace*) extrudedE[0];
+    eFace->addPhysicalEntity(500);
+    current->setPhysicalName("outerWall", 2, 500);
+    GRegion *eRegion = (GRegion*) extrudedE[1];
+    eRegion->addPhysicalEntity(600);
+    current->setPhysicalName("wallVolume", 3, 600);
+    for (int j= 2; j < extrudedE.size(); j++){
+      GFace *elFace = (GFace*) extrudedE[j];
+      std::list<GEdge*> l_edges = elFace->edges();
+      for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
+	GEdge *myEdge = *it;
+	if (is_cut) myEdge = current->getEdgeByTag((*it)->tag()-NE);
+	if( std::find(boundEdges.begin(), boundEdges.end(), myEdge) != boundEdges.end() ){
+	  if (myEdge==gin){
+	    elFace->addPhysicalEntity(700);
+	    current->setPhysicalName("inletRing", 2, 700);
+	  }
+	  else{
+	    elFace->addPhysicalEntity(800);
+	    current->setPhysicalName("outletRing", 2, 800);
+	  }
+	}
+      } 
+    }
   }
 
 }
@@ -898,35 +895,45 @@ void Centerline::run(){
 
   if (is_cut) cutMesh();
   else{
+    GFace *gf = current->getFaceByTag(1);
+    gf->addPhysicalEntity(100);
+    current->setPhysicalName("wall", 2, 100);
     current->createTopologyFromMesh();
     NV = current->getMaxElementaryNumber(0);
     NE = current->getMaxElementaryNumber(1);
     NF = current->getMaxElementaryNumber(2);
     NR = current->getMaxElementaryNumber(3);  
   }
+  
+  //identify the boundary edges by looping over all discreteFaces
+  std::vector<GEdge*> boundEdges;
+  double dist_inlet = 1.e6;
+  GEdge *gin = NULL;
+  for (int i= 0; i< NF; i++){
+    GFace *gf = current->getFaceByTag(i+1);
+    std::list<GEdge*> l_edges = gf->edges();
+    for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
+      std::vector<GEdge*>::iterator ite = std::find(boundEdges.begin(), boundEdges.end(), *it);
+      if (ite != boundEdges.end()) boundEdges.erase(ite);
+      else boundEdges.push_back(*it);
+      GVertex *gv = (*it)->getBeginVertex();
+      SPoint3 pt(gv->x(), gv->y(), gv->z());
+      double dist = pt.distance(ptin);
+      if(dist < dist_inlet){
+	dist_inlet = dist;
+	gin = *it;
+      }			      
+    }
+  }
 
-  if (is_closed) createClosedVolume();
-  if (is_extruded) extrudeBoundaryLayerWall();
+  if (is_closed)   createClosedVolume(gin, boundEdges);
+  if (is_extruded) extrudeBoundaryLayerWall(gin, boundEdges);
 
 }
 
 void Centerline::cutMesh(){
 
-  is_cut = 1;
-
-  // std::vector<GFace*> currentFaces =  current->bindingsGetFaces();
-  // for (int i=0; i< currentFaces.size(); i++){
-  //   printf("gf =%d \n", currentFaces[i]->tag());
-  //   if(currentFaces[i]->geomType() == GEntity::CompoundSurface) {
-  //     std::list<GFace*> cFaces = ((GFaceCompound*)currentFaces[i])->getCompounds();
-  //     for (std::list<GFace*>::iterator it = cFaces.begin(); it!= cFaces.end(); it++) {
-  //   	printf("comp  = %d\n", (*it)->tag());
-  //     }
-  //     if (cFaces.size() > 0) currentGFC.push_back(currentFaces[i]);
-  //   }
-  // }
-
-  Msg::Info("Splitting surface mesh (%d tris) with centerline %s ", triangles.size(), fileName.c_str());
+  Msg::Info("Centerline: action (cutMesh) splits surface mesh (%d tris) using %s ", triangles.size(), fileName.c_str());
 
   //splitMesh
   for(unsigned int i = 0; i < edges.size(); i++){
@@ -934,13 +941,14 @@ void Centerline::cutMesh(){
     double L = edges[i].length;
     double D = (edges[i].minRad+edges[i].maxRad);
     double AR = L/D;
-    printf("*** Centerline branch %d (AR=%g):  ", edges[i].tag, AR);
-    printf("children (%d) = ", edges[i].children.size());
-    for (int k= 0; k< edges[i].children.size() ; k++) printf("%d ", edges[i].children[k].tag);
-    printf("\n");
+    printf("*** Centerline branch %d (AR=%.1f) \n", edges[i].tag, AR);
+    //if ( edges[i].children.size()) printf("children (%d) = ", edges[i].children.size());
+    //for (int k= 0; k< edges[i].children.size() ; k++) printf("%d ", edges[i].children[k].tag);
+    //printf("\n");
    
-    int nbSplit = (int)floor(AR/2. + 0.5); 
-    if( AR > 3.0){
+    int nbSplit = (int)floor(AR/2 + 0.9); 
+    if( nbSplit > 1 ){
+      printf("->> cut branch in %d parts \n",  nbSplit);
       double li  = L/nbSplit;
       double lc = 0.0;
       for (unsigned int j= 0; j < lines.size(); j++){
@@ -950,7 +958,6 @@ void Centerline::cutMesh(){
 	  MVertex *v2 = lines[j]->getVertex(1);
 	  SVector3 pt(v1->x(), v1->y(), v1->z());
 	  SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
-	  printf("->> cut length %d split \n",  nbSplit);
 	  std::map<MLine*,double>::iterator itr = radiusl.find(lines[j]);
 	  bool cutted = cutByDisk(pt, dir, itr->second);
 	  nbSplit--;
@@ -966,19 +973,19 @@ void Centerline::cutMesh(){
       else v2 = lines[lines.size()-1]->getVertex(0);
       SVector3 pt(v1->x(), v1->y(), v1->z());
       SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
-      printf("-->> cut bifurcation \n");
+      printf("-->> cut branch at bifurcation \n");
       std::map<MLine*,double>::iterator itr = radiusl.find(lines[lines.size()-1]);
       bool cutted = cutByDisk(pt, dir, itr->second);
-      if(!cutted){
-	int l = lines.size()-1-lines.size()/(4*nbSplit);
-	v1 = lines[l]->getVertex(1);
-	v2 = lines[l]->getVertex(0);
-	pt = SVector3(v1->x(), v1->y(), v1->z());
-	dir = SVector3(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
-	printf("-->> cut bifurcation NEW \n");
-	itr = radiusl.find(lines[l]);
-	cutted = cutByDisk(pt, dir, itr->second);
-      }
+      // if(!cutted){
+      // 	int l = lines.size()-1-lines.size()/(4*nbSplit); //chech this!
+      // 	v1 = lines[l]->getVertex(1);
+      // 	v2 = lines[l]->getVertex(0);
+      // 	pt = SVector3(v1->x(), v1->y(), v1->z());
+      // 	dir = SVector3(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
+      // 	printf("-->> cut bifurcation NEW \n");
+      // 	itr = radiusl.find(lines[l]);
+      // 	cutted = cutByDisk(pt, dir, itr->second);
+      // }
     }
  }
 
@@ -1006,7 +1013,7 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
   double d = -a * PT.x() - b * PT.y() - c * PT.z();
   //printf("cut disk (R=%g)= %g %g %g %g \n", maxRad, a, b, c, d);
 
-  int maxStep = 40;
+  int maxStep = 20;
   const double EPS = 0.007;
 
   std::set<MEdge,Less_Edge> allEdges;
@@ -1169,7 +1176,7 @@ void  Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
    //dir_n = normal direction of the disk
    //dir_t = tangential direction of the disk
    SVector3 dir_a = p1-p0;
-   SVector3 dir_n(x-p0.x(), y-p0.y(), z-p0.z()); 
+   SVector3 dir_n(xyz[0]-p0.x(), xyz[1]-p0.y(), xyz[2]-p0.z()); 
    SVector3 dir_t;
    buildOrthoBasis2(dir_a, dir_n, dir_t);
    
@@ -1204,10 +1211,12 @@ void  Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
    return;
 }
 
-/**************************************************/
-/******************Temporary code******************/
-/**************************************************/
+//used so far by Tristan for mesh algo hex
 void Centerline::operator()(double x,double y,double z,SVector3& v1,SVector3& v2,SVector3& v3,GEntity* ge){
+
+  printf("EMI do not go here \n");
+  exit(1);
+
   if(update_needed){
     std::ifstream input;
 	input.open(fileName.c_str());
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index 532d4749a1ff3faa6facdfcb7b1a27eacc8baf2d..14225c6206d08ae843a8091a54c13a6897fe8b0d 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -147,9 +147,9 @@ class Centerline : public Field{
   // Cut the mesh in different parts of small aspect ratio
   void cutMesh();
   //Create In and Outlet Planar Faces
-  void createClosedVolume();
+  void createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges);
   //extrude outer wall
-  void extrudeBoundaryLayerWall();
+  void extrudeBoundaryLayerWall(GEdge *gin, std::vector<GEdge*> boundEdges);
 
   // Cut the tubular structure with a disk
   // perpendicular to the tubular structure
diff --git a/benchmarks/boolean/aneurysmBL.py b/benchmarks/boolean/aneurysmBL.py
index 68b0b95745f1c5efc18e8a38e6a9dedb29b3f64e..c5b15484aef5249d1751ff0f51d1bc34279b4ce1 100644
--- a/benchmarks/boolean/aneurysmBL.py
+++ b/benchmarks/boolean/aneurysmBL.py
@@ -11,11 +11,11 @@ g.load("aneurysm.stl");
 g.createTopologyFromMesh();
 
 face = g.getFaceByTag(1);
-newface  = g.extrudeBoundaryLayer(face, 4, 0.5, 0, -1)
-newface2 = g.extrudeBoundaryLayer(face, 4, -0.5, 1, -1)
+newent  = g.extrudeBoundaryLayer(face, 4, 0.5, 0, -1)
+newent2 = g.extrudeBoundaryLayer(face, 4, -0.5, 1, -1)
 
 print "*** face = %d " % face.tag()
-print "*** new face = %d newface2 = %d " % (newface.tag(), newface2.tag())
+print "*** new face = %d newface2 = %d " % (newent[0].tag(), newent2[0].tag())
 
 g.mesh(2)
 g.save("aneurysmBL.msh")