diff --git a/benchmarks/3d/Falcon/CreateFullInitialMesh.geo b/benchmarks/3d/Falcon/CreateFullInitialMesh.geo
new file mode 100644
index 0000000000000000000000000000000000000000..c9e9d527c181ece2b49d69da319b9e7c89cd2bb3
--- /dev/null
+++ b/benchmarks/3d/Falcon/CreateFullInitialMesh.geo
@@ -0,0 +1,25 @@
+Merge "InitialMeshFalcon.msh";
+Plugin(NewView).Run;
+Plugin(Transform).A33 = -1;
+Plugin(Transform).Tz = 2 * 0.272342;
+Plugin(Transform).SwapOrientation = 1;
+Plugin(Transform).Run;
+Physical Surface(1) = 1; // keep same num
+For i In {2 : 16}
+Physical Surface(100 + i) = i;
+EndFor
+Save "tmp.msh";
+
+Delete All;
+
+Mesh.SwitchElementTags = 1;
+Merge "tmp.msh";
+
+Mesh.SwitchElementTags = 0;
+Merge "InitialMeshFalcon.msh";
+
+Geometry.Tolerance = 1.e-4;
+Coherence Mesh;
+
+Mesh.SaveAll = 1;
+Save "FullInitialMeshFalcon.msh";
diff --git a/benchmarks/3d/Falcon/FullSurfaceMeshUniform.geo b/benchmarks/3d/Falcon/FullSurfaceMeshUniform.geo
new file mode 100644
index 0000000000000000000000000000000000000000..c6c2549cd4e1697449b2a9d342b045c68fc329d3
--- /dev/null
+++ b/benchmarks/3d/Falcon/FullSurfaceMeshUniform.geo
@@ -0,0 +1,33 @@
+Mesh.RemeshParametrization=7; //0=harmonic_circle, 1=conformal, 2=rbf, 3=harmonic_plane, 4=convex_circle, 5=convex_plane, 6=harmonic square
+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 = lc;
+Mesh.CharacteristicLengthMax = lc;
+Mesh.LcIntegrationPrecision = 1.e-5;
+Mesh.MinimumCirclePoints = 50;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthFromCurvature = 0;
+Mesh.CharacteristicLengthFromPoints = 0;
+
+Merge "FullInitialMeshFalcon.msh";
+
+CreateTopology;
+
+// Make all lines compound:
+ll[] = Line "*";
+For j In {0 : #ll[]-1}
+  Compound Line(1000+j) = ll[j];
+  Physical Line(1000+j) = (1000+j);
+EndFor
+
+// Make all surfaces compound and physical:
+ss[] = Surface "*";
+For i In {0 : #ss[]-1}
+  Compound Surface(i+1000) = ss[i];
+  Physical Surface(i+1) = { i+1000 };
+EndFor
+
diff --git a/benchmarks/3d/Falcon/InitialMeshFalcon.msh b/benchmarks/3d/Falcon/InitialMeshFalcon.msh
index a540ab3d54c05d0d8de493de37e3730411520d99..636ec86fedb0ce4871a4496a4f3171e7638835ab 100644
--- a/benchmarks/3d/Falcon/InitialMeshFalcon.msh
+++ b/benchmarks/3d/Falcon/InitialMeshFalcon.msh
@@ -1296,7 +1296,7 @@ $Nodes
 1291 0.354756 0.244145 0.296104
 1292 0.34252 0.233767 0.300558
 1293 0.354628 0.209661 0.282031
-1294 0.379376 0.211619 0.273384
+1294 0.379376 0.211619 0.2723384
 1295 0.345558 0.201179 0.272907
 1296 0.318762 0.193141 0.2727
 1297 0.32524 0.19931 0.279823
@@ -1309,7 +1309,7 @@ $Nodes
 1304 0.326889 0.220634 0.299753
 1305 0.334657 0.210759 0.290358
 1306 0.365673 0.218769 0.28853
-1307 0.416977 0.224773 0.273925
+1307 0.416977 0.224773 0.2723925
 1308 0.38391 0.218836 0.281475
 1309 0.384859 0.223013 0.28583
 1310 0.40557 0.228232 0.281989
@@ -1317,21 +1317,21 @@ $Nodes
 1312 0.402727 0.236694 0.285627
 1313 0.372881 0.235135 0.294069
 1314 0.382643 0.252874 0.282154
-1315 0.415678 0.250532 0.274051
+1315 0.415678 0.250532 0.2724051
 1316 0.354416 0.258542 0.282622
-1317 0.379211 0.25956 0.273504
-1318 0.349445 0.265407 0.273048
+1317 0.379211 0.25956 0.2723504
+1318 0.349445 0.265407 0.2723048
 1319 0.328055 0.264848 0.280326
-1320 0.325848 0.269618 0.273014
+1320 0.325848 0.269618 0.2723014
 1321 0.304768 0.265833 0.285063
-1322 0.299346 0.274122 0.273039
+1322 0.299346 0.274122 0.2723039
 1323 0.308748 0.256158 0.296493
 1324 0.333901 0.255337 0.291573
 1325 0.364403 0.2508 0.289575
 1326 0.383753 0.248697 0.286816
 1327 0.381965 0.243776 0.290108
 1328 0.404806 0.245625 0.282579
-1329 0.434555 0.238519 0.274184
+1329 0.434555 0.238519 0.2724184
 1330 0.106071 0.130031 0.272436
 1331 0.136716 0.153951 0.300276
 1332 0.127068 0.159146 0.315543
diff --git a/benchmarks/3d/Falcon/falcon.script b/benchmarks/3d/Falcon/falcon.script
index 272b619c52bf19b97b9ca5890af37f64918d1861..061a741cfa607726456c6e29e352ee617e70128d 100644
--- a/benchmarks/3d/Falcon/falcon.script
+++ b/benchmarks/3d/Falcon/falcon.script
@@ -32,10 +32,10 @@ For i In {0 : #ss[]-1}
 EndFor
 
 Mesh 2;
-Save "SurfaceMeshUniform.stl";
+Save "tmp.stl";
 Delete All;
 
-Merge "SurfaceMeshUniform.stl";
+Merge "tmp.stl";
 
 CreateTopology;
 
diff --git a/benchmarks/3d/Falcon/fullfalcon.script b/benchmarks/3d/Falcon/fullfalcon.script
new file mode 100644
index 0000000000000000000000000000000000000000..a3e47431d6b12a454245f0997a43f3d49a1afae9
--- /dev/null
+++ b/benchmarks/3d/Falcon/fullfalcon.script
@@ -0,0 +1,129 @@
+Include "param.geo";
+
+Mesh.RemeshParametrization=7;
+Mesh.RemeshAlgorithm=1;
+Mesh.Algorithm=6;
+
+Mesh.CharacteristicLengthFromPoints = 0;
+Mesh.CharacteristicLengthMin = lc;
+Mesh.CharacteristicLengthMax = lc;
+Mesh.LcIntegrationPrecision = 1.e-5;
+Mesh.MinimumCirclePoints = 50;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthFromCurvature = 0;
+Mesh.CharacteristicLengthFromPoints = 0;
+
+Merge "FullInitialMeshFalcon.msh";
+
+CreateTopology;
+
+// Make all lines compound:
+ll[] = Line "*";
+For j In {0 : #ll[]-1}
+  Compound Line(1000+j) = ll[j];
+  Physical Line(1000+j) = (1000+j);
+EndFor
+
+// Make all surfaces compound and physical:
+ss[] = Surface "*";
+For i In {0 : #ss[]-1}
+  Compound Surface(i+1000) = ss[i];
+  Physical Surface(i+1) = { i+1000 };
+EndFor
+
+Mesh 2;
+Save "SurfaceMeshUniform.stl";
+
+If(TWO_PLANES)
+Plugin(NewView).Run;
+x = -Pi/10;
+y = -Pi/10;
+z = -Pi/10;
+A = Cos(x); B = Sin(x); C = Cos(y); D = Sin(y); E = Cos(z); F = Sin(z);
+AD = A * D; BD = B * D;
+Plugin(Transform).A11 = C*E; 
+Plugin(Transform).A12 = BD*E+A*F;
+Plugin(Transform).A13 =-AD*E+B*F;
+Plugin(Transform).A21 =-C*F;
+Plugin(Transform).A22 =-BD*F+A*E;
+Plugin(Transform).A23 = AD*F+B*E;
+Plugin(Transform).A31 = D;
+Plugin(Transform).A32 =-B*C;
+Plugin(Transform).A33 = A*C;
+Plugin(Transform).Tx = 11;
+Plugin(Transform).Ty = 1;
+Plugin(Transform).Tz = -1;
+Plugin(Transform).Run;
+Save "SurfaceMeshUniform2.stl";
+EndIf
+
+Delete All;
+
+Include "param.geo";
+
+Merge "SurfaceMeshUniform.stl";
+If(TWO_PLANES)
+  Merge "SurfaceMeshUniform2.stl";
+EndIf
+
+zmin = General.MinZ-dd;
+zmax = General.MaxZ+dd;
+xmin = General.MinX-dd;
+xmax = General.MaxX+dd;
+ymin = General.MinY-dd;
+ymax = General.MaxY+dd;
+
+Point(20001) = {xmin, ymin, zmin};
+Point(20002) = {xmax, ymin, zmin};
+Point(20003) = {xmax, ymax, zmin};
+Point(20004) = {xmin, ymax, zmin};
+Point(20005) = {xmin, ymin, zmax};
+Point(20006) = {xmax, ymin, zmax};
+Point(20007) = {xmax, ymax, zmax};
+Point(20008) = {xmin, ymax, zmax};
+
+Line(10) = {20001, 20002};
+Line(11) = {20002, 20003};
+Line(12) = {20003, 20004};
+Line(13) = {20004, 20001};
+Line(14) = {20005, 20006};
+Line(15) = {20006, 20007};
+Line(16) = {20007, 20008};
+Line(17) = {20008, 20005};
+Line(18) = {20001, 20005};
+Line(19) = {20002, 20006};
+Line(20) = {20003, 20007};
+Line(21) = {20004, 20008};
+
+Line Loop(22) = {10,19,-14,-18};
+Line Loop(23) = {11,20,-15,-19};
+Line Loop(24) = {12,21,-16,-20};
+Line Loop(25) = {13,18,-17,-21};
+Line Loop(26) = {14,15,16,17};
+Line Loop(27) = {10,11,12,13};
+
+Plane Surface(22) = {22};
+Plane Surface(23) = {23};
+Plane Surface(24) = {24};
+Plane Surface(25) = {25};
+Plane Surface(26) = {26};
+Plane Surface(27) = {27};
+
+Surface Loop(20) = {22,23,24,25,26,27};
+Surface Loop(1) = {1};
+
+If(!TWO_PLANES)
+  Volume(1) = {20,1};
+  Physical Surface(1) = {1};
+EndIf
+If(TWO_PLANES)
+  Surface Loop(2) = {2};
+  Volume(1) = {20,1,2};
+  Physical Surface(1) = {1,2};
+EndIf
+
+Physical Surface(20) = {22:27};
+Physical Volume(40) = {1};
+
+Mesh 3;
+Save "fullfalcon.msh";