diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index cb40e80db9f39cc868abeda76984043f657d1a29..9075defe56ad2ca2073ee8c77c819569a8907df3 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -594,6 +594,7 @@ Surface *Create_Surface(int Num, int Typ)
   pS->TransfiniteSmoothing = -1;
   pS->TrsfPoints = List_Create(4, 4, sizeof(Vertex *));
   pS->Generatrices = NULL;
+  pS->GeneratricesByTag = NULL;
   pS->EmbeddedPoints = NULL;
   pS->EmbeddedCurves = NULL;
   pS->Extrude = NULL;
@@ -607,6 +608,7 @@ static void Free_Surface(void *a, void *b)
   if(pS) {
     List_Delete(pS->TrsfPoints);
     List_Delete(pS->Generatrices);
+    List_Delete(pS->GeneratricesByTag);
     List_Delete(pS->EmbeddedCurves);
     List_Delete(pS->EmbeddedPoints);
     delete pS;
@@ -1023,9 +1025,11 @@ static void CopySurface(Surface *s, Surface *ss, bool copyMeshingMethod)
     if(List_Nbr(s->TrsfPoints))
       Msg::Warning("Only automatic transfinite surface specifications can be copied");
   }
-  ss->Generatrices = List_Create(List_Nbr(s->Generatrices), 1, sizeof(Curve *));
+  ss->Generatrices = List_Create(List_Nbr(s->Generatrices) + 1, 1, sizeof(Curve *));
+  ss->GeneratricesByTag = List_Create(List_Nbr(s->GeneratricesByTag) + 1, 1, sizeof(int));
   ss->InSphereCenter = s->InSphereCenter; // FIXME: hack...
   List_Copy(s->Generatrices, ss->Generatrices);
+  List_Copy(s->GeneratricesByTag, ss->GeneratricesByTag);
   End_Surface(ss);
 }
 
@@ -3626,7 +3630,10 @@ void setSurfaceEmbeddedCurves(Surface *s, List_T *curves)
 void setSurfaceGeneratrices(Surface *s, List_T *loops)
 {
   int nbLoop = List_Nbr(loops);
+  List_Delete(s->Generatrices);
   s->Generatrices = List_Create(4, 4, sizeof(Curve *));
+  List_Delete(s->GeneratricesByTag);
+  s->GeneratricesByTag = List_Create(4, 4, sizeof(int));
   for(int i = 0; i < nbLoop; i++) {
     int iLoop;
     List_Read(loops, i, &iLoop);
@@ -3647,10 +3654,9 @@ void setSurfaceGeneratrices(Surface *s, List_T *loops)
           ic *= sign(iLoop);
           if(i != 0) ic *= -1; // hole
           if(!(c = FindCurve(ic))) {
-            Msg::Error("Unknown curve %d", ic);
             List_Delete(s->Generatrices);
             s->Generatrices = NULL;
-            return;
+            break;
           }
           else
             List_Add(s->Generatrices, &c);
@@ -3665,12 +3671,25 @@ void setSurfaceGeneratrices(Surface *s, List_T *loops)
             Msg::Error("Unknown curve %d", ic);
             List_Delete(s->Generatrices);
             s->Generatrices = NULL;
-            return;
+            break;
           }
           else
             List_Add(s->Generatrices, &c);
         }
       }
+      if(!s->Generatrices){
+        for(int j = 0; j < List_Nbr(el->Curves); j++) {
+          List_Read(el->Curves, j, &ic);
+          GEdge *ge = GModel::current()->getEdgeByTag(abs(ic));
+          if(ge) {
+            List_Add(s->GeneratricesByTag, &ic);
+          }
+          else{
+            Msg::Error("Unknown curve %d", ic);
+            return;
+          }
+        }
+      }
     }
   }
 }
diff --git a/Geo/Geo.h b/Geo/Geo.h
index 5de728a88eeb81d966b34a5c6bdffb73d48816b5..77dff77e62d1a56dbe2e297cf9feeb17f329f066 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -152,6 +152,7 @@ class Surface{
   double RecombineAngle;
   int TransfiniteSmoothing;
   List_T *Generatrices;
+  List_T *GeneratricesByTag;
   List_T *EmbeddedCurves;
   List_T *EmbeddedPoints;
   List_T *TrsfPoints;
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index ea18b06527c3d73fb9d71b6af0f483b7d2c6421c..cf0b529beef3bd4fd502a3090f884cafc3d5b669 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -24,34 +24,53 @@ gmshFace::gmshFace(GModel *m, Surface *face)
 
   edgeCounterparts = s->edgeCounterparts;
 
-  std::list<GEdge*> l_wire;
-  GVertex *first = 0;
+  std::vector<GEdge*> eds;
+  std::vector<int> nums;
   for(int i = 0; i < List_Nbr(s->Generatrices); i++){
     Curve *c;
     List_Read(s->Generatrices, i, &c);
     GEdge *e = m->getEdgeByTag(abs(c->Num));
     if(e){
-      GVertex *start = (c->Num > 0) ? e->getBeginVertex() : e->getEndVertex();
-      GVertex *next  = (c->Num > 0) ? e->getEndVertex() : e->getBeginVertex();
-      if (!first) first = start;
-      l_wire.push_back(e);
-      if (next == first){
-        edgeLoops.push_back(GEdgeLoop(l_wire));
-        l_wire.clear();
-        first = 0;
-      }
-
-      l_edges.push_back(e);
-      e->addFace(this);
-      l_dirs.push_back((c->Num > 0) ? 1 : -1);
-      if (List_Nbr(s->Generatrices) == 2){
-        e->meshAttributes.minimumMeshSegments =
-          std::max(e->meshAttributes.minimumMeshSegments, 2);
-      }
+      eds.push_back(e);
+      nums.push_back(c->Num);
     }
     else
       Msg::Error("Unknown curve %d", c->Num);
   }
+  for(int i = 0; i < List_Nbr(s->GeneratricesByTag); i++){
+    int j;
+    List_Read(s->GeneratricesByTag, i, &j);
+    GEdge *e = m->getEdgeByTag(abs(j));
+    if(e){
+      eds.push_back(e);
+      nums.push_back(j);
+    }
+    else
+      Msg::Error("Unknown curve %d", j);
+  }
+
+  std::list<GEdge*> l_wire;
+  GVertex *first = 0;
+  for(unsigned int i = 0; i < eds.size(); i++){
+    GEdge *e = eds[i];
+    int num = nums[i];
+    GVertex *start = (num > 0) ? e->getBeginVertex() : e->getEndVertex();
+    GVertex *next  = (num > 0) ? e->getEndVertex() : e->getBeginVertex();
+    if (!first) first = start;
+    l_wire.push_back(e);
+    if (next == first){
+      edgeLoops.push_back(GEdgeLoop(l_wire));
+      l_wire.clear();
+      first = 0;
+    }
+    l_edges.push_back(e);
+    e->addFace(this);
+    l_dirs.push_back((num > 0) ? 1 : -1);
+    if (List_Nbr(s->Generatrices) == 2){
+      e->meshAttributes.minimumMeshSegments =
+        std::max(e->meshAttributes.minimumMeshSegments, 2);
+    }
+  }
 
   // always compute and store the mean plane for plane surfaces (using
   // the bounding vertices)
@@ -331,7 +350,7 @@ bool gmshFace::containsPoint(const SPoint3 &pt) const
 bool gmshFace::buildSTLTriangulation(bool force)
 {
   return false;
-  
+
   if(va_geom_triangles){
     if(force)
       delete va_geom_triangles;
@@ -355,8 +374,8 @@ bool gmshFace::buildSTLTriangulation(bool force)
     CTX::instance()->mesh.lcFactor = 1;
     CTX::instance()->mesh.order = 1;
     CTX::instance()->mesh.lcIntegrationPrecision = 1.e-3;
-    //  CTX::instance()->mesh.Algorithm = 5;    
-    model()->mesh(2);    
+    //  CTX::instance()->mesh.Algorithm = 5;
+    model()->mesh(2);
     CTX::instance()->mesh = _temp;
     fields->setBackgroundField(fields->get(BGM));
   }
@@ -400,5 +419,5 @@ bool gmshFace::buildSTLTriangulation(bool force)
   }
   va_geom_triangles->finalize();
   return true;
-  
+
 }
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index 499fcca1f3ee40ec8e96ab448f522c6b975d5bd0..fcec42f11be6a8aebe08b96dd90959cfbcdc8ed7 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -39,7 +39,6 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
       Msg::Error("Unknown surface %d", is);
   }
   resetMeshAttributes();
-
 }
 
 void gmshRegion::resetMeshAttributes()
diff --git a/benchmarks/3d/Falcon/CreateVolumeAroundFalcon_FrontalISO.geo b/benchmarks/3d/Falcon/CreateVolumeAroundFalcon_FrontalISO.geo
index e4261b2f195797a667d2e3649951c06a37bd111e..f04278cec39844c2e2ec6f5966b1da22a6824ed3 100644
--- a/benchmarks/3d/Falcon/CreateVolumeAroundFalcon_FrontalISO.geo
+++ b/benchmarks/3d/Falcon/CreateVolumeAroundFalcon_FrontalISO.geo
@@ -1,19 +1,21 @@
-Mesh.Algorithm3D=4; //1-Delaunay, 4-Frontal, 5-Frontal Delaunay 6-Frontal hex, 7-MMG3D, 9-RTree
+//Mesh.Algorithm3D=4; //1-Delaunay, 4-Frontal, 5-Frontal Delaunay 6-Frontal hex, 7-MMG3D, 9-RTree
 
-Mesh.CharacteristicLengthFactor = 1.2/10;
+lc = 0.3;
+Mesh.CharacteristicLengthMin = lc;
+Mesh.CharacteristicLengthMax = lc;
 
-Merge "InitialMeshFalcon.stl";
-//Merge "SurfaceMeshUniform.stl";
+//Merge "InitialMeshFalcon.stl";
+Merge "SurfaceMeshUniform.stl";
 //Merge "SurfaceMeshAnisoCurvature.stl";
 
 CreateTopology;
 
 z0 = 0.272342;
-zmax = 15.0;
-xmin = -27.0;
-xmax = 10.32;
-ymin = -10.0;
-ymax = 10.0;
+zmax = 11;
+xmin = -20;
+xmax = 3;
+ymin = -3.0;
+ymax = 6.0;
 
 Point(20001) = {xmin, ymin, z0};
 Point(20002) = {xmax, ymin, z0};
@@ -54,3 +56,7 @@ Plane Surface(7) = {27};
 Surface Loop(8) = {1,2,3,4,5,6,7};
 Volume(1) = {8};
 
+Physical Surface(1) = {1};
+Physical Surface(2) = {5, 4, 3, 2, 6};
+Physical Surface(3) = {7};
+Physical Volume(4) = {1};
diff --git a/benchmarks/3d/Falcon/SurfaceMeshUniform.geo b/benchmarks/3d/Falcon/SurfaceMeshUniform.geo
index e15d7b85a7d087deaedf5d4cc1aaa3b2dd2d7098..5af988447e5339b20265a754eb7df58e232d2850 100644
--- a/benchmarks/3d/Falcon/SurfaceMeshUniform.geo
+++ b/benchmarks/3d/Falcon/SurfaceMeshUniform.geo
@@ -2,9 +2,11 @@ Mesh.RemeshParametrization=7; //0=harmonic_circle, 1=conformal, 2=rbf, 3=harmoni
 Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split only with metis
 Mesh.Algorithm=6; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
 
+lc = 0.3;
+
 Mesh.CharacteristicLengthFromPoints = 0;
-Mesh.CharacteristicLengthMin = 0.02;
-Mesh.CharacteristicLengthMax = 0.02;
+Mesh.CharacteristicLengthMin = lc;
+Mesh.CharacteristicLengthMax = lc;
 Mesh.LcIntegrationPrecision = 1.e-5;
 Mesh.MinimumCirclePoints = 50;
 Mesh.CharacteristicLengthExtendFromBoundary = 0;
@@ -28,3 +30,4 @@ For i In {0 : #ss[]-1}
   Compound Surface(i+100) = ss[i];
   Physical Surface(i+1) = { i+100 };
 EndFor
+