diff --git a/benchmarks/levelset/cube.geo b/benchmarks/levelset/cube.geo
index 45c1076681f5f53c79436abd4f5695d24791aff7..ce096cfa80949ee372b642ab1f0aa34f56aa22d3 100644
--- a/benchmarks/levelset/cube.geo
+++ b/benchmarks/levelset/cube.geo
@@ -1,3 +1,4 @@
+
 lc = 0.1;
 Point(1) = {0.0,0.0,0.0,lc};
 Point(2) = {1,0.0,0.0,lc};
@@ -12,25 +13,11 @@ Plane Surface(6) = {5};
 Extrude {0,0.0,1} {
   Surface{6};
 }
-
-Mesh 3;
-Levelset Plane (1) = {0,-1,0,0.5};
-Levelset Plane (2) = {-1,0,0,0.5};
-Levelset Sphere (3) = {{0.75,0.5,0.5},0.5};
-Levelset Ellipsoid (4) = { {0,0,0}, {1,0,0}, 0.45, 0.25, 0.75 };
-Levelset MathEval (5) = "(x-0.5)^2+(y-0.5)^2+(z-0.5)^2-0.5^2";
-Levelset MathEval (6) = "-((x-0.5)^2+(y-0.5)^2+(z-0.5)^2-0.4^2)";
-Levelset Sphere (7) = {{0.75,0.5,0.5},0.4};
-Levelset Intersection (10) = {1,2,3};
-Levelset Intersection (11) = {5,6};
-Levelset Cut (12) = {3,7};
-
-Levelset Point(13)={{0.5,0,0},{0.5,0.5,0.},{0.5,1,0},{0.5,0,0.5},{0.5,0,0.8},{0.5,0.5,0.5},{0.5,0.5,1},{0.5,1,0.5},{1.5,1,1}};
-
-Levelset CutMeshTri {5};
-
-Print "cube_cut.msh";
-
-
-
+Physical Surface("back") = {6};
+Physical Surface("right") = {15};
+Physical Surface("bottom") = {19};
+Physical Surface("left") = {23};
+Physical Surface("top") = {27};
+Physical Surface("front") = {28};
+Physical Volume("vol")={1};
 
diff --git a/benchmarks/levelset/test.geo b/benchmarks/levelset/test.geo
new file mode 100644
index 0000000000000000000000000000000000000000..ac1635f559ad798836d19e32b05d6e4f575c7bd2
--- /dev/null
+++ b/benchmarks/levelset/test.geo
@@ -0,0 +1,40 @@
+lc = 0.1;
+Point(1) = {0.0,0.0,0.0,lc};
+Point(2) = {1,0.0,0.0,lc};
+Point(3) = {1,1,0.0,lc};
+Point(4) = {0,1,0.0,lc};
+Line(1) = {4,3};
+Line(2) = {3,2};
+Line(3) = {2,1};
+Line(4) = {1,4};
+Line Loop(5) = {2,3,4,1};
+Plane Surface(6) = {5};
+Extrude {0,0.0,1} {
+  Surface{6};
+}
+
+Physical Point (1) = {1};
+Physical Point (2) = {2};
+Physical Point (3) = {3};
+Physical Point (4) = {4};
+Physical Line (1) = {1};
+Physical Line (2) = {2};
+Physical Line (3) = {3};
+Physical Line (4) = {4};
+Physical Surface("back") = {6};
+Physical Surface("right") = {15};
+Physical Surface("bottom") = {19};
+Physical Surface("left") = {23};
+Physical Surface("top") = {27};
+Physical Surface("front") = {28};
+
+Physical Volume("vol") = {1};
+
+Mesh 3;
+
+// Definition of the levelsets
+Levelset MathEval (1) = "x^2+y^2-0.8";
+Levelset Plane (2) = {0,0,1,-0.5};
+Levelset Union (3) = {1,2};
+Levelset CutMesh {3};
+
diff --git a/benchmarks/levelset/test.py b/benchmarks/levelset/test.py
new file mode 100644
index 0000000000000000000000000000000000000000..5b7f394c0b1085d6e5cc1b606eb483c96e5903d5
--- /dev/null
+++ b/benchmarks/levelset/test.py
@@ -0,0 +1,20 @@
+
+from gmshpy import *
+import math
+
+#Define levelset
+ls1 = gLevelsetMathEval('x^2+y^2-0.8',1)
+ls2 = gLevelsetPlane(0,0,1,-0.5,2)
+ls3 = gLevelsetUnion([ls1,ls2])
+		
+#Mesh geometry
+mesh = GModel()
+mesh.load("cube.geo")
+mesh.mesh(3)
+
+#Cut mesh
+model1 = GModel()
+model1 = mesh.buildCutGModel(ls3, True, False) # ls, cutElements?, write triangles?
+
+model1.save("cube2.msh")
+