diff --git a/benchmarks/3d/Falcon/README.txt b/benchmarks/3d/Falcon/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..eb3850f17276095f3e746b3f833018ad3214063a
--- /dev/null
+++ b/benchmarks/3d/Falcon/README.txt
@@ -0,0 +1,14 @@
+// for the symmetric falcon:
+
+gmsh falcon_surface.script - // generates remeshed falcon boundary mesh(es)
+gmsh falcon_volume.script - // generates volume mesh in box around falcon(s)
+gmsh partition.script - // partitions the mesh and saves using new physicals
+gmsh split.script - // splits the partitionned mesh into separate files
+
+// for the full falcon(s)
+
+gmsh fullfalcon_surface.script - // generates remeshed falcon boundary mesh(es)
+gmsh fullfalcon_volume.script - // generates volume mesh in box around falcon(s)
+gmsh partition.script - // partitions the mesh and saves using new physicals
+gmsh fullsplit.script - // splits the partitionned mesh into separate files
+
diff --git a/benchmarks/3d/Falcon/falcon_surface.script b/benchmarks/3d/Falcon/falcon_surface.script
new file mode 100644
index 0000000000000000000000000000000000000000..e1a209f8247bebcbbcda3b7e70bfef7bfe1d2ea7
--- /dev/null
+++ b/benchmarks/3d/Falcon/falcon_surface.script
@@ -0,0 +1,33 @@
+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;
+
+Merge "InitialMeshFalcon.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+100) = ss[i];
+ Physical Surface(i+1) = { i+100 };
+EndFor
+
+Mesh 2;
+Save "tmp.stl";
diff --git a/benchmarks/3d/Falcon/falcon.script b/benchmarks/3d/Falcon/falcon_volume.script
similarity index 50%
rename from benchmarks/3d/Falcon/falcon.script
rename to benchmarks/3d/Falcon/falcon_volume.script
index af884b711366b41c611e56f392c53f731d91eb87..53f9b838ecf35f8856991ed25341c894844fd3fa 100644
--- a/benchmarks/3d/Falcon/falcon.script
+++ b/benchmarks/3d/Falcon/falcon_volume.script
@@ -1,37 +1,4 @@
-lc = 0.3;
-
-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;
-
-Merge "InitialMeshFalcon.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+100) = ss[i];
- Physical Surface(i+1) = { i+100 };
-EndFor
-
-Mesh 2;
-Save "tmp.stl";
-Delete All;
+Include "param.geo";
Merge "tmp.stl";
@@ -44,14 +11,14 @@ xmax = 3;
ymin = -3.0;
ymax = 6.0;
-Point(20001) = {xmin, ymin, z0};
-Point(20002) = {xmax, ymin, z0};
-Point(20003) = {xmax, ymax, z0};
-Point(20004) = {xmin, ymax, z0};
-Point(20005) = {xmin, ymin, zmax};
-Point(20006) = {xmax, ymin, zmax};
-Point(20007) = {xmax, ymax, zmax};
-Point(20008) = {xmin, ymax, zmax};
+Point(20001) = {xmin, ymin, z0, lc2};
+Point(20002) = {xmax, ymin, z0, lc2};
+Point(20003) = {xmax, ymax, z0, lc2};
+Point(20004) = {xmin, ymax, z0, lc2};
+Point(20005) = {xmin, ymin, zmax, lc2};
+Point(20006) = {xmax, ymin, zmax, lc2};
+Point(20007) = {xmax, ymax, zmax, lc2};
+Point(20008) = {xmin, ymax, zmax, lc2};
Line(10) = {20001, 20002};
Line(11) = {20002, 20003};
diff --git a/benchmarks/3d/Falcon/fullfalcon.script b/benchmarks/3d/Falcon/fullfalcon.script
deleted file mode 100644
index a3e47431d6b12a454245f0997a43f3d49a1afae9..0000000000000000000000000000000000000000
--- a/benchmarks/3d/Falcon/fullfalcon.script
+++ /dev/null
@@ -1,129 +0,0 @@
-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";
diff --git a/benchmarks/3d/Falcon/fullfalcon_surface.script b/benchmarks/3d/Falcon/fullfalcon_surface.script
new file mode 100644
index 0000000000000000000000000000000000000000..046560f0dbadda6d48c0b7696abcb1da3dd56e1e
--- /dev/null
+++ b/benchmarks/3d/Falcon/fullfalcon_surface.script
@@ -0,0 +1,60 @@
+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;
diff --git a/benchmarks/3d/Falcon/fullfalcon_volume.script b/benchmarks/3d/Falcon/fullfalcon_volume.script
new file mode 100644
index 0000000000000000000000000000000000000000..b953e41d9fe2605714065b715a428518eec67342
--- /dev/null
+++ b/benchmarks/3d/Falcon/fullfalcon_volume.script
@@ -0,0 +1,69 @@
+
+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, lc2};
+Point(20002) = {xmax, ymin, zmin, lc2};
+Point(20003) = {xmax, ymax, zmin, lc2};
+Point(20004) = {xmin, ymax, zmin, lc2};
+Point(20005) = {xmin, ymin, zmax, lc2};
+Point(20006) = {xmax, ymin, zmax, lc2};
+Point(20007) = {xmax, ymax, zmax, lc2};
+Point(20008) = {xmin, ymax, zmax, lc2};
+
+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(2) = {22:27};
+Physical Volume(4) = {1};
+
+Mesh 3;
+Save "falcon.msh";
diff --git a/benchmarks/3d/Falcon/fullsplit.script b/benchmarks/3d/Falcon/fullsplit.script
new file mode 100644
index 0000000000000000000000000000000000000000..0c03bb9efeea709abf8828c32ef7f65b5b8d9e32
--- /dev/null
+++ b/benchmarks/3d/Falcon/fullsplit.script
@@ -0,0 +1,36 @@
+Include "param.geo";
+Merge "part.msh";
+N_DOM=NUM_SLICES;
+
+vol[] = {}; idx_vol[] = {};
+sur_scat[] = {}; idx_sur_scat[] = {};
+sur_inf[] = {}; idx_sur_inf[] = {};
+sigma_left[] = {}; idx_sigma_left[] = {};
+sigma_right[] = {}; idx_sigma_right[] = {};
+For i In {1:N_DOM}
+ idx_vol[] += #vol[]; vol[] += Physical Volume{4000 + i};
+ idx_sur_scat[] += #sur_scat[]; sur_scat[] += Physical Surface{1000 + i};
+ idx_sur_inf[] += #sur_inf[]; sur_inf[] += Physical Surface{2000 + i};
+ idx_sigma_left[] += #sigma_left[]; sigma_left[] += Physical Surface{3000 + 1000 * (i-2)};
+ idx_sigma_right[] += #sigma_right[]; sigma_right[] += Physical Surface{3000 + 1000 * (i-1)};
+EndFor
+idx_vol[] += #vol[];
+idx_sur_scat[] += #sur_scat[];
+idx_sur_inf[] += #sur_inf[];
+idx_sigma_left[] += #sigma_left[];
+idx_sigma_right[] += #sigma_right[];
+
+Mesh.Binary = 1;
+Geometry.OrientedPhysicals = 0; // partitions have negative tags!!
+
+For i In {1:N_DOM}
+ idom = i-1;
+ Delete Physicals;
+ Physical Volume(4000+i) = vol[{idx_vol[i-1]:idx_vol[i]-1}];
+ Physical Surface(1000+i) = sur_scat[{idx_sur_scat[i-1]:idx_sur_scat[i]-1:1}];
+ Physical Surface(2000+i) = sur_inf[{idx_sur_inf[i-1]:idx_sur_inf[i]-1:1}];
+ Physical Surface(-(4000 + 1000 * (i-2))) = sigma_left[{idx_sigma_left[i-1]:idx_sigma_left[i]-1:1}];
+ Physical Surface(4000 + 1000 * (i-1)) = sigma_right[{idx_sigma_right[i-1]:idx_sigma_right[i]-1:1}];
+ Save Sprintf("falcon_mshcut%g.msh", idom);
+EndFor
+
diff --git a/benchmarks/3d/Falcon/param.geo b/benchmarks/3d/Falcon/param.geo
new file mode 100644
index 0000000000000000000000000000000000000000..4e3a2894aee5f159755d9bd397de871db291f596
--- /dev/null
+++ b/benchmarks/3d/Falcon/param.geo
@@ -0,0 +1,6 @@
+lc = 0.25; // surface
+lc2 = lc * 5; // volume
+dd = 2;
+TWO_PLANES = 1;
+NUM_SLICES = 4;
+
diff --git a/benchmarks/3d/Falcon/partition.script b/benchmarks/3d/Falcon/partition.script
index e1043cfa4c3be613c1687bf4921aaf563b553fd6..049372918f80bdfba7305431bb131fd94190f619 100644
--- a/benchmarks/3d/Falcon/partition.script
+++ b/benchmarks/3d/Falcon/partition.script
@@ -1,6 +1,8 @@
+Include "param.geo";
+
Merge "falcon.msh";
-Plugin(SimplePartition).NumSlices = 16;
+Plugin(SimplePartition).NumSlices = NUM_SLICES;
Plugin(SimplePartition).Run;
Mesh.Binary=1;
diff --git a/benchmarks/3d/Falcon/split.script b/benchmarks/3d/Falcon/split.script
index 2e948846b96c28fa45cf0b4645d7b1ac4e947658..29b41b3f56238ddbc2b7ab4afd56a98b0d103535 100644
--- a/benchmarks/3d/Falcon/split.script
+++ b/benchmarks/3d/Falcon/split.script
@@ -1,6 +1,6 @@
-//Merge "part.msh_new";
+Include "param.geo";
Merge "part.msh";
-N_DOM=16;
+N_DOM=NUM_SLICES;
vol[] = {}; idx_vol[] = {};
sur_scat[] = {}; idx_sur_scat[] = {};