diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 19d4c8ec..ae7e0095 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -138,7 +138,12 @@ test_win64:
 .template_test_linux: &test_linux
   stage: test
   script:
-    - docker run netgen_${CI_BUILD_REF_NAME}_installed:${UBUNTU_VERSION} bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V"'
+    - >-
+      docker run
+      -e NETGENDIR=/opt/netgen/bin
+      -e PYTHONPATH=/opt/netgen/lib/python3/dist-packages
+      netgen_${CI_BUILD_REF_NAME}_installed:${UBUNTU_VERSION}
+      bash -c 'cd /root/build/netgen && make test_netgen ARGS="-V"'
 
 test_ubuntu_1510:
   <<: *ubuntu_1510
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f2bbd4e2..1befe02f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -258,6 +258,11 @@ execute_process(COMMAND hdiutil create -volname Netgen -srcfolder ${INSTALL_DIR}
 ")
   add_custom_target(bundle COMMAND ${CMAKE_COMMAND} "-P" "${CMAKE_CURRENT_BINARY_DIR}/fixup.cmake")
 #######################################################################
+# CTest
+enable_testing()
+include(CTest)
+
+#######################################################################
 
 add_subdirectory(libsrc)
 add_subdirectory(ng)
@@ -269,13 +274,9 @@ add_subdirectory(py_tutorials)
 add_subdirectory(doc)
 add_subdirectory(windows)
 add_subdirectory(nglib)
+add_subdirectory(tests)
 
 
-#######################################################################
-# CTest
-enable_testing()
-include(CTest)
-
 #######################################################################
 # Debian packager
 
diff --git a/libsrc/csg/csg.hpp b/libsrc/csg/csg.hpp
index d524903f..ec3ad1e1 100644
--- a/libsrc/csg/csg.hpp
+++ b/libsrc/csg/csg.hpp
@@ -10,7 +10,7 @@
 #include <myadt.hpp>
 #include <gprim.hpp>
 #include <meshing.hpp>
-
+#include <vector>
 // #include <geometry2d.hpp>
 #include "../gprim/spline.hpp"
 #include "../gprim/splinegeometry.hpp"
@@ -21,12 +21,12 @@
 #include "solid.hpp"
 #include "identify.hpp"
 #include "singularref.hpp"
+#include "splinesurface.hpp"
 #include "csgeom.hpp"
 #include "csgparser.hpp"
 
 #include "triapprox.hpp"
 #include "algprim.hpp"
-#include "splinesurface.hpp"
 #include "brick.hpp"
 #include "spline3d.hpp"
 #include "manifold.hpp"
diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp
index 0cb2c8c2..892d754f 100644
--- a/libsrc/csg/csgeom.cpp
+++ b/libsrc/csg/csgeom.cpp
@@ -345,7 +345,25 @@ namespace netgen
 	const ExtrusionFace * ef = dynamic_cast< const ExtrusionFace * > (GetSurface(i));
 	const RevolutionFace * rf = dynamic_cast< const RevolutionFace * > (GetSurface(i));
 	const DummySurface * dummyf = dynamic_cast< const DummySurface * > (GetSurface(i));
+        const SplineSurface * splines = dynamic_cast<const SplineSurface *> (GetSurface(i));
 
+        if (splines)
+          {
+            splines->GetBase()->GetPrimitiveData(classname,coeffs);
+            out << classname << " " << coeffs.Size() << "\n";
+            for (int j=0; j<coeffs.Size(); j++)
+              out << coeffs[j] << " ";
+            out << "\n";
+            for (auto cut : *(splines->GetCuts()))
+              {
+                cut->GetPrimitiveData(classname,coeffs);
+                out << classname << " " << coeffs.Size() << "\n";
+                for (int j=0; j<coeffs.Size(); j++)
+                  out << coeffs[j] << " ";
+                out << "\n";
+              }
+            return;
+          }
 
 	if(sp)
 	  {
@@ -819,6 +837,7 @@ namespace netgen
     int inv;
     int nsurf = GetNSurf();
 
+
     isidenticto.SetSize(nsurf);
     for (int i = 0; i < nsurf; i++)
       isidenticto[i] = i;
diff --git a/libsrc/csg/csgeom.hpp b/libsrc/csg/csgeom.hpp
index 66e3e7cb..3106326e 100644
--- a/libsrc/csg/csgeom.hpp
+++ b/libsrc/csg/csgeom.hpp
@@ -143,6 +143,9 @@ namespace netgen
     /// filename of inputfile
     string filename;
 
+    /// store splinesurfaces, such that added ones do not get deleted before geometry does
+    std::vector<shared_ptr<SplineSurface>> spline_surfaces;
+
   public:
     CSGeometry ();
     CSGeometry (const string & afilename);
@@ -315,7 +318,9 @@ namespace netgen
 
     virtual int GenerateMesh (shared_ptr<Mesh> & mesh, MeshingParameters & mparam);
 
-    virtual const Refinement & GetRefinement () const; 
+    virtual const Refinement & GetRefinement () const;
+
+    void AddSplineSurface (shared_ptr<SplineSurface> ss) { spline_surfaces.push_back(ss); }
   };
 
 
diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp
index 68cac826..30ac3b85 100644
--- a/libsrc/csg/edgeflw.cpp
+++ b/libsrc/csg/edgeflw.cpp
@@ -491,7 +491,7 @@ namespace netgen
 	    if(splinesurface)
 	      {
 		auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p);
-		mesh.SetCD2Name(refedges[i].edgenr,*name);
+		mesh.SetCD2Name(refedges[i].edgenr,name);
 	      }
 	    else
 	      {
@@ -499,7 +499,7 @@ namespace netgen
 	    if(splinesurface2)
 	      {
 		auto name = splinesurface2->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p);
-		mesh.SetCD2Name(refedges[i].edgenr,*name);
+		mesh.SetCD2Name(refedges[i].edgenr,name);
 	      }
 		
 	      }
diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp
index b8c4ffbf..b61a9934 100644
--- a/libsrc/csg/python_csg.cpp
+++ b/libsrc/csg/python_csg.cpp
@@ -207,11 +207,11 @@ DLL_HEADER void ExportCSG(py::module &m)
     ;
 
   py::class_<SplineSurface, shared_ptr<SplineSurface>> (m, "SplineSurface",
-			 "A surface for co dim 2 integrals on the splines")
+                        "A surface for co dim 2 integrals on the splines")
     .def("__init__", FunctionPointer  ([](SplineSurface* instance, shared_ptr<SPSolid> base, py::list cuts)
 	     {
 	       auto primitive = dynamic_cast<OneSurfacePrimitive*> (base->GetSolid()->GetPrimitive());
-	       auto acuts = new Array<OneSurfacePrimitive*>();
+	       auto acuts = make_shared<std::vector<shared_ptr<OneSurfacePrimitive>>>();
 	       for(int i = 0; i<py::len(cuts);i++)
 		 {
 		   py::extract<shared_ptr<SPSolid>> sps(cuts[i]);
@@ -219,13 +219,14 @@ DLL_HEADER void ExportCSG(py::module &m)
 		     throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!");
 		   auto sp = dynamic_cast<OneSurfacePrimitive*>(sps()->GetSolid()->GetPrimitive());
 		   if(sp)
-		     acuts->Append(sp);
+		     acuts->push_back(shared_ptr<OneSurfacePrimitive>(sp));
 		   else
 		     throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!");
 		 }
 	       if(!primitive)
 		 throw NgException("Base is not a SurfacePrimitive in constructor of SplineSurface!");
-	       new (instance) SplineSurface(primitive,acuts);
+	       new (instance) SplineSurface(shared_ptr<OneSurfacePrimitive>(primitive),acuts);
+               py::object obj = py::cast(instance);
 	     }),py::arg("base"), py::arg("cuts")=py::list())
     .def("AddPoint", FunctionPointer
 	 ([] (SplineSurface & self, double x, double y, double z, bool hpref)
@@ -237,8 +238,8 @@ DLL_HEADER void ExportCSG(py::module &m)
     .def("AddSegment", FunctionPointer
 	 ([] (SplineSurface & self, int i1, int i2, string bcname, double maxh)
 	  {
-	    auto str = new string(bcname);
-	    self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str,maxh);
+            auto seg = make_shared<LineSeg<3>>(self.GetPoint(i1),self.GetPoint(i2));
+	    self.AppendSegment(seg,bcname,maxh);
 	  }),
 	 py::arg("pnt1"),py::arg("pnt2"),py::arg("bcname")="default", py::arg("maxh")=-1.)
     ;
@@ -433,14 +434,15 @@ DLL_HEADER void ExportCSG(py::module &m)
 	 ([] (CSGeometry & self, shared_ptr<SplineSurface> surf)
 	  {
 	    auto cuttings = surf->CreateCuttingSurfaces();
-	    auto spsol = make_shared<SPSolid>(new Solid(&*surf));
+	    auto spsol = make_shared<SPSolid>(new Solid(surf.get()));
 	    for(auto cut : (*cuttings)){
-	      spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(cut)));
+	      spsol = make_shared<SPSolid>(SPSolid::SECTION,spsol,make_shared<SPSolid>(new Solid(cut.get())));
 	    }
 	    spsol->AddSurfaces(self);
-	    int tlonr = self.SetTopLevelObject(spsol->GetSolid(), &*surf);
+	    int tlonr = self.SetTopLevelObject(spsol->GetSolid(), surf.get());
 	    for(auto p : surf->GetPoints())
 		self.AddUserPoint(p);
+            self.AddSplineSurface(surf);
 	  }),
 	  py::arg("SplineSurface"))
 
diff --git a/libsrc/csg/splinesurface.cpp b/libsrc/csg/splinesurface.cpp
index d66c8bff..65a0f988 100644
--- a/libsrc/csg/splinesurface.cpp
+++ b/libsrc/csg/splinesurface.cpp
@@ -6,22 +6,22 @@ namespace netgen
 void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const bool hpref)
 {
   auto pp = Point<3>(p);
-  geompoints.Append(GeomPoint<3>(pp,reffac));
-  geompoints.Last().hpref = hpref;
+  geompoints.push_back(GeomPoint<3>(pp,reffac));
+  geompoints.back().hpref = hpref;
 }
   
-  void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh)
+  void SplineSurface :: AppendSegment(shared_ptr<SplineSeg<3>> sp, string & bcname, double amaxh)
   {
-    splines.Append(spline);
-    bcnames.Append(bcname);
+    splines.push_back(sp);
+    bcnames.push_back(bcname);
     maxh.Append(amaxh);
   }
 
-  string* SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const
+  string SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const
   {
     
     double eps = 1e-5;
-    for(int i=0; i<splines.Size(); i++)
+    for(int i=0; i<splines.size(); i++)
       {
 	auto pp1 = Point<3>(splines[i]->GetPoint(0));
 	Project(pp1);
@@ -32,36 +32,44 @@ void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const
 	    return bcnames[i];
 	  }
       }
-    return new string("default");
+    return "default";
   }
 
-  Array<OneSurfacePrimitive*>* SplineSurface :: CreateCuttingSurfaces() const
+  const shared_ptr<std::vector<shared_ptr<OneSurfacePrimitive>>> SplineSurface :: CreateCuttingSurfaces()
   {
-    auto cuttings = new Array<OneSurfacePrimitive*>();
+    if(all_cuts)
+      return all_cuts;
+    auto cuttings = make_shared<std::vector<shared_ptr<OneSurfacePrimitive>>>();
     for (auto cut : *cuts)
-      cuttings->Append(cut);
-    for(int i = 0; i<splines.Size(); i++)
+      cuttings->push_back(cut);
+    for(int i = 0; i<splines.size(); i++)
       {
 	auto spline = splines[i];
-	auto lineseg = dynamic_cast<LineSeg<3>*>(spline);
-	auto p1 = Point<3>(spline->GetPoint(0));
-	Project(p1);
-	auto p2 = Point<3>(spline->GetPoint(1));
-	Project(p2);
-	auto vec = Vec<3>(p2)-Vec<3>(p1);
-	auto plane = new Plane(p1,-Cross(vec,baseprimitive->GetNormalVector(p1)));
-	if(maxh[i]>0)
-	  {
-	  plane->SetMaxH(maxh[i]);
-	  }
-	cuttings->Append(plane);
+	auto lineseg = dynamic_cast<LineSeg<3>*>(spline.get());
+        if(lineseg)
+          {
+            auto p1 = Point<3>(spline->GetPoint(0));
+            Project(p1);
+            auto p2 = Point<3>(spline->GetPoint(1));
+            Project(p2);
+            auto vec = Vec<3>(p2)-Vec<3>(p1);
+            auto plane = make_shared<Plane>(p1,-Cross(vec,baseprimitive->GetNormalVector(p1)));
+            if(maxh[i]>0)
+              {
+                plane->SetMaxH(maxh[i]);
+              }
+            cuttings->push_back(plane);
+          }
+        else
+          throw NgException("Spline type not implemented for SplineSurface!");
       }
+    all_cuts = cuttings;
     return cuttings;
   }
   
   void SplineSurface :: Print(ostream & str) const
 {
-  str << "SplineSurface " << endl;
+  str << "SplineSurface with base " << *baseprimitive << endl;
 }
 
 }
diff --git a/libsrc/csg/splinesurface.hpp b/libsrc/csg/splinesurface.hpp
index 1f68ae73..0bc5494c 100644
--- a/libsrc/csg/splinesurface.hpp
+++ b/libsrc/csg/splinesurface.hpp
@@ -7,37 +7,37 @@ namespace netgen
   class SplineSurface  : public  OneSurfacePrimitive
   {
   protected:
-    Array<GeomPoint<3>> geompoints;
-    Array<SplineSeg<3>*> splines;
-    Array<string*> bcnames;
+    std::vector<GeomPoint<3>> geompoints;
+    std::vector<shared_ptr<SplineSeg<3>>> splines;
+    std::vector<string> bcnames;
     Array<double> maxh;
-    OneSurfacePrimitive* baseprimitive;
-    Array<OneSurfacePrimitive*>* cuts;
+    shared_ptr<OneSurfacePrimitive> baseprimitive;
+    shared_ptr<std::vector<shared_ptr<OneSurfacePrimitive>>> cuts;
+    shared_ptr<std::vector<shared_ptr<OneSurfacePrimitive>>> all_cuts;
     
   public:
-    SplineSurface(OneSurfacePrimitive* abaseprimitive, Array<OneSurfacePrimitive*>* acuts) :
+    SplineSurface(shared_ptr<OneSurfacePrimitive> abaseprimitive, shared_ptr<std::vector<shared_ptr<OneSurfacePrimitive>>> acuts) :
       OneSurfacePrimitive(), baseprimitive(abaseprimitive), cuts(acuts)
     { ; }
     virtual ~SplineSurface() { ; }
     
-    const Array<SplineSeg<3>*> & GetSplines() const { return splines; }
-    int GetNSplines() const { return splines.Size(); }
-    const Array<GeomPoint<3>>& GetPoints() const { return geompoints; }
+    const auto & GetSplines() const { return splines; }
+    int GetNSplines() const { return splines.size(); }
+    const std::vector<GeomPoint<3>>& GetPoints() const { return geompoints; }
     string GetSplineType(const int i) const { return splines[i]->GetType(); }
     SplineSeg<3> & GetSpline(const int i) { return *splines[i]; }
     const SplineSeg<3> & GetSpline(const int i) const { return *splines[i]; }
-    int GetNP() const { return geompoints.Size(); }
+    int GetNP() const { return geompoints.size(); }
     const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; }
-    string* GetBCName(int i) const { return bcnames[i]; }
-    string* GetBCNameOf(Point<3> p1, Point<3> p2) const;
+    string GetBCName(int i) const { return bcnames[i]; }
+    string GetBCNameOf(Point<3> p1, Point<3> p2) const;
     
     DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false);
-    void AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh = -1);
-
-    OneSurfacePrimitive* GetBase() const { return baseprimitive; } 
-
-    Array<OneSurfacePrimitive*>* CreateCuttingSurfaces() const;
+    void AppendSegment(shared_ptr<SplineSeg<3>> spline, string & bcname, double amaxh = -1);
 
+    const shared_ptr<std::vector<shared_ptr<OneSurfacePrimitive>>> CreateCuttingSurfaces();
+    const shared_ptr<std::vector<shared_ptr<OneSurfacePrimitive>>> GetCuts() const { return all_cuts; }
+    const shared_ptr<OneSurfacePrimitive> GetBase() const { return baseprimitive; }
     
     virtual void Project (Point<3> & p3d) const { baseprimitive->Project(p3d); }
     virtual double CalcFunctionValue (const Point<3> & point) const
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index 9af43005..21585384 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -75,7 +75,7 @@ namespace netgen
 	 (*testout) << "Meshing subdomain " << k << endl;
 	 
 	 mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
-	 
+
 	 mesh3d.CalcSurfacesOfNode();
 	 mesh3d.FindOpenElements(k);
 	 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
new file mode 100644
index 00000000..f898217f
--- /dev/null
+++ b/tests/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(pytest)
diff --git a/tests/docker_15.10 b/tests/docker_15.10
index e159dfa8..a88f5673 100644
--- a/tests/docker_15.10
+++ b/tests/docker_15.10
@@ -1,4 +1,4 @@
 FROM ubuntu:15.10
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
-RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev libboost-python-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache
+RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy
 ADD . /root/src/netgen
diff --git a/tests/docker_16.04 b/tests/docker_16.04
index 2ab970e7..a1fe001e 100644
--- a/tests/docker_16.04
+++ b/tests/docker_16.04
@@ -1,4 +1,4 @@
 FROM ubuntu:16.04
 MAINTAINER Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>
-RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev libboost-python-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache
+RUN apt-get update && apt-get -y install python3 libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev ccache python3-pytest python3-numpy
 ADD . /root/src/netgen
diff --git a/tests/pytest/CMakeLists.txt b/tests/pytest/CMakeLists.txt
new file mode 100644
index 00000000..9997622a
--- /dev/null
+++ b/tests/pytest/CMakeLists.txt
@@ -0,0 +1,4 @@
+if(USE_PYTHON)
+    add_test(NAME pytest COMMAND ${PYTHON_EXECUTABLE} -m pytest WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+    add_custom_target(pytest ${PYTHON_EXECUTABLE} -m pytest WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+endif(USE_PYTHON)
diff --git a/tests/pytest/test_savemesh.py b/tests/pytest/test_savemesh.py
new file mode 100644
index 00000000..6b7b7e2c
--- /dev/null
+++ b/tests/pytest/test_savemesh.py
@@ -0,0 +1,52 @@
+
+from netgen.csg import *
+from netgen import meshing
+import filecmp
+import difflib
+from math import sqrt, cos, sin
+
+def CreateQuad():
+    base = Plane(Pnt(0,0,0),Vec(0,0,1))
+    surface = SplineSurface(base)
+    pts = [(-0.2,-0.2,0),(-0.2,0.2,0),(0.2,0.2,0),(0.2,-0.2,0)]
+    geopts = [surface.AddPoint(*p) for p in pts]
+    for p1,p2,bc in [(0,1,"wire"), (1, 2,"contact"),(2,3,"wire"),(3,0,"wire")]:
+        surface.AddSegment(geopts[p1],geopts[p2],bc)
+    return surface
+
+Cross = lambda a,b: [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-b[0]*a[1]]
+
+
+def CreateGeo():
+    geo = CSGeometry()
+    air = OrthoBrick(Pnt(-1,-1,-1),Pnt(1,1,1))
+    geo.Add(air.mat("air"))
+    surface = CreateQuad()
+    geo.AddSplineSurface(surface)
+    return geo
+
+def test_BBNDsave():
+    mesh = CreateGeo().GenerateMesh(maxh=0.4,perfstepsend = meshing.MeshingStep.MESHSURFACE)
+    for i in range(2):
+        mesh.GenerateVolumeMesh(mp = MeshingParameters(only3D_domain=i+1,maxh=0.4))
+    mesh.SetGeometry(None)
+    mesh.Save("test.vol")
+    mesh2 = meshing.Mesh()
+    mesh2.Load("test.vol")
+    mesh2.Save("test2.vol")
+    with open("test.vol","r") as f:
+        first = f.readlines()
+    with open("test2.vol","r") as f:
+        second = f.readlines()
+    # exclude the face colours section (because they aren't in the same order)
+    for i,line in enumerate(first):
+        if line[0:12] == "face_colours":
+            first = first[0:i]
+            second = second[0:i]
+            break
+    diff = difflib.context_diff(first,second)
+    print("File diff:")
+    l = list(diff)
+    print(*l)
+    assert len(l)==0
+