diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 205d5074..2a990358 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -48,7 +48,7 @@ test_win:
   stage: test
   script:
     - cd %NETGEN_BUILD_DIR%\netgen
-    - ctest -C Release -V
+    - ctest -C Release -V --output-on-failure
     - cd ..
 
 cleanup_win:
@@ -106,7 +106,7 @@ test_ubuntu:
       -e NETGENDIR=/opt/netgen/bin
       -e PYTHONPATH=/opt/netgen/lib/python3/dist-packages
       netgen_${CI_PIPELINE_ID}_installed:${UBUNTU_VERSION}
-      bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V"'
+      bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V --output-on-failure"'
 
 test_ubuntu_mpi:
   <<: *ubuntu
@@ -117,7 +117,7 @@ test_ubuntu_mpi:
       -e NETGENDIR=/opt/netgen/bin
       -e PYTHONPATH=/opt/netgen/lib/python3/dist-packages
       netgen_mpi_${CI_PIPELINE_ID}_installed:${UBUNTU_VERSION}
-      bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V"'
+      bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V --output-on-failure"'
 
 test_build_ngsolve:
   <<: *ubuntu
@@ -205,7 +205,7 @@ test_mac:
   stage: test
   script:
     - cd $BUILD_DIR/netgen
-    - ctest . --output-on-failure
+    - ctest . -V --output-on-failure
 
 cleanup_mac:
   <<: *mac
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index b00a630e..1e111a7d 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -3270,12 +3270,8 @@ namespace netgen
     for (int i = 0; i < surfelements.Size(); i++)
       if (surfelements[i].IsDeleted())
         {
-          surfelements[i] = surfelements.Last();
-          surfelements.DeleteLast();
-          /*
-          surfelements.Delete(i);
+          surfelements.DeleteElement(i);
           i--;
-          */
         }
 
     for (int i = 0; i < segments.Size(); i++)
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 4646e3b4..e72dc282 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -569,6 +569,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     */
     
     .def_property_readonly("_timestamp", &Mesh::GetTimeStamp)
+    .def_property_readonly("ne", [](Mesh& m) { return m.GetNE(); })
     .def("Distribute", [](shared_ptr<Mesh> self, NgMPI_Comm comm) {
 	self->SetCommunicator(comm);
 	if(comm.Size()==1) return self;
diff --git a/tests/pytest/results.py b/tests/pytest/results.py
new file mode 100644
index 00000000..39f82c79
--- /dev/null
+++ b/tests/pytest/results.py
@@ -0,0 +1,34 @@
+number_elements = {}
+number_elements['cylsphere.geo'] = (243,508,706,2665,16573)
+number_elements['cubeandspheres.geo'] = (100,98,98,366,1078)
+number_elements['ellipsoid.geo'] = (624,596,1255,5384,35591)
+number_elements['manyholes2.geo'] = (126463,291489)
+number_elements['sculpture.geo'] = (140,260,477,1316,6395)
+number_elements['ortho.geo'] = (6,6,6,28,180)
+number_elements['ellipticcone.geo'] = (574,1750,4635,12578,67070)
+number_elements['cube.geo'] = (6,6,6,28,178)
+number_elements['twobricks.geo'] = (22,22,42,177,587)
+number_elements['revolution.geo'] = (1237,3779,7958,31475,192314)
+number_elements['circle_on_cube.geo'] = (39,173,586,1860,11364)
+number_elements['sphereincube.geo'] = (178,328,508,1584,12778)
+number_elements['twocubes.geo'] = (22,22,42,177,587)
+number_elements['boundarycondition.geo'] = (22,22,39,165,507)
+number_elements['ellipticcyl.geo'] = (325,1090,2113,7594,52464)
+number_elements['trafo.geo'] = (1358,2384,5081,16891,110333)
+number_elements['boxcyl.geo'] = (146,370,820,3427,17236)
+number_elements['sphere.geo'] = (56,80,126,347,2303)
+number_elements['torus.geo'] = (2518,2715,5286,24144,169082)
+number_elements['shaft.geo'] = (863,1743,2589,10855,59980)
+number_elements['cone.geo'] = (441,706,1173,4298,25551)
+number_elements['cubeandring.geo'] = (226,568,1942,7083,35417)
+number_elements['manyholes.geo'] = (28638,69054)
+number_elements['period.geo'] = (587,1362,3136,10876,64321)
+number_elements['lshape3d.geo'] = (12,12,18,83,317)
+number_elements['cubemsphere.geo'] = (737,1414,4090,16091,105619)
+number_elements['twocyl.geo'] = (151,418,578,1838,12863)
+number_elements['cubemcyl.geo'] = (3102,7714,17803,81423,491425)
+number_elements['matrix.geo'] = (1936,2765,5005,15244,93613)
+number_elements['fichera.geo'] = (18,18,35,211,477)
+number_elements['cylinder.geo'] = (84,272,397,1152,7696)
+number_elements['part1.stl'] = (328,501,1739,4054,76653)
+number_elements['hinge.stl'] = (789,1130,2588,6595,125874)
diff --git a/tests/pytest/test_tutorials.py b/tests/pytest/test_tutorials.py
new file mode 100644
index 00000000..cde78728
--- /dev/null
+++ b/tests/pytest/test_tutorials.py
@@ -0,0 +1,75 @@
+
+import os, pytest
+from netgen.meshing import meshsize
+import netgen.gui
+import netgen.csg as csg
+import netgen.stl as stl
+import netgen.csg as csg
+from results import *
+
+def getFiles(fileEnding):
+    r, d, files = next(os.walk(os.path.join("..","..","tutorials")))
+    return (f for f in files if f.endswith(fileEnding))
+
+
+def getCheckFunc(filename):
+    def func(mesh,i):
+        if filename in number_elements:
+            # number of elements should be in 2% range of expected value
+            assert mesh.ne == pytest.approx(number_elements[filename][i], rel=0.02)
+        return func
+            
+def getResultFunc(filename):
+    def resultFunc(mesh):
+        results = {}
+        results["number_elements"] = mesh.ne
+        return results
+    return resultFunc
+
+def getMeshingparameters(filename):
+    standard = (meshsize.very_coarse, meshsize.coarse, meshsize.moderate, meshsize.fine, meshsize.very_fine)
+    if filename == "shell.geo":
+        return [] # do not test this example cause it needs so long...
+    if filename == "extrusion.geo":
+        return [] # this segfaults right now
+    if filename == "manyholes2.geo" or filename == "manyholes.geo":
+        return standard[:2] # this gets too big for finer meshsizes
+    return standard
+
+# don't test step files as they do not respect all meshing parameters correctly yet.
+_geofiles = [f for f in getFiles(".geo")] + [f for f in getFiles(".stl")] # + [f for f in getFiles(".step")]
+
+def generateMesh(filename, mp):
+    if filename.endswith(".geo"):
+        geo = csg.CSGeometry(os.path.join("..","..","tutorials", filename))
+    elif filename.endswith(".stl"):
+        geo = stl.STLGeometry(os.path.join("..","..","tutorials", filename))
+    elif filename.endswith(".step"):
+        geo = occ.OCCGeometry(os.path.join("..","..","tutorials", filename))
+    return geo.GenerateMesh(mp)
+
+@pytest.mark.parametrize("filename, checkFunc", [(f, getCheckFunc(f)) for f in _geofiles])
+def test_geoFiles(filename, checkFunc):
+    for i, mp in enumerate(getMeshingparameters(filename)):
+        print("load geo", filename)
+        mesh = generateMesh(filename, mp)
+        if checkFunc is not None:
+            checkFunc(mesh,i)
+
+import time
+def generateResultFile():
+    with open("results.py", "w") as f:
+        print("number_elements = {}", file=f)
+        for _file, _func in ((gf, getResultFunc(gf)) for gf in _geofiles):
+            start = time.time()
+            print("write", _file)
+            mps = getMeshingparameters(_file)
+            if not mps:
+                continue
+            results = [_func(generateMesh(_file, mp)) for mp in mps]
+            print("number_elements['{}'] = {}".format(_file, "(" + ",".join((str(r["number_elements"]) for r in results)) + ")"), file=f)
+            print("needed", time.time() - start, "seconds")
+        
+
+if __name__ == "__main__":
+    generateResultFile()