diff --git a/libsrc/interface/writevtk.cpp b/libsrc/interface/writevtk.cpp
index 4e6f80ea..ff7ba9fe 100644
--- a/libsrc/interface/writevtk.cpp
+++ b/libsrc/interface/writevtk.cpp
@@ -50,64 +50,68 @@ void WriteVtkFormat (const Mesh & mesh,
   outfile << "POINTS " << np << " double\n";
   for (int i=0; i<np; i++)
   {
-	auto & p = mesh.Point(i+1);
+  auto & p = mesh.Point(i+1);
     outfile << p[0] << " " << p[1] << " " << p[2] << "\n";
   }
 
   std::vector<int> types;
+  std::vector<int> domains;
   if (ne > 0)
   {
     unsigned int size = 0;
     for (int i=0; i<ne; i++)
-	  size += mesh.VolumeElement(i+1).GetNV() + 1; // only save "linear" corners
+      size += mesh.VolumeElement(i+1).GetNV() + 1; // only save "linear" corners
 
     outfile << "CELLS " << ne << " " << size << "\n";
     for (int i=0; i<ne; i++)
     {
-	  auto& el = mesh.VolumeElement(i+1);
-	  switch (el.GetType())
-	  {
-	  case TET:
-	  case TET10: // reorder to follow VTK convention & zero based indices
+      auto& el = mesh.VolumeElement(i+1);
+      domains.push_back(el.GetIndex());
+      switch (el.GetType())
+      {
+      case TET:
+      case TET10: // reorder to follow VTK convention & zero based indices
         outfile << 4 << " " << el[0]-1 << " " << el[1]-1 << " " << el[3]-1 << " " << el[2]-1 << "\n";
-		types.push_back(10);
-		break;
-	  case PRISM: // reorder to follow VTK convention & zero based indices
+        types.push_back(10);
+        break;
+      case PRISM: // reorder to follow VTK convention & zero based indices
         outfile << 6 << " "
-				<< el[0]-1 << " " << el[2]-1 << " " << el[1]-1 << " " 
-				<< el[3]-1 << " " << el[5]-1 << " " << el[4]-1 << "\n";
-		types.push_back(13);
-		break;
-	  default:
-	    throw ngcore::Exception("Unexpected element type");
-		break;
-	  }
+          << el[0]-1 << " " << el[2]-1 << " " << el[1]-1 << " " 
+          << el[3]-1 << " " << el[5]-1 << " " << el[4]-1 << "\n";
+        types.push_back(13);
+        break;
+      default:
+        throw ngcore::Exception("Unexpected element type");
+        break;
+      }
     }
   }
   else
   {
     unsigned int size = 0;
-    for (int i=0; i<ne; i++)
-	  size += mesh.SurfaceElement(i+1).GetNV() + 1;
+    for (int i=0; i<nse; i++)
+      size += mesh.SurfaceElement(i+1).GetNV() + 1;
 
     outfile << "CELLS " << nse << " " << size << "\n";
     for (int i=0; i<nse; i++)
     {
-	  auto& el = mesh.SurfaceElement(i+1);
-	  switch (el.GetType())
-	  {
-	  case TRIG:
-	  case TRIG6:
-        outfile << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << "\n";
-		types.push_back(5);
-		break;
-	  case QUAD:
-        outfile << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << " " << el[3]-1 << "\n";
-		types.push_back(9);
-	  default:
-	    throw ngcore::Exception("Unexpected element type");
-		break;
-	  }
+      auto& el = mesh.SurfaceElement(i+1);
+      domains.push_back(el.GetIndex());
+      switch (el.GetType())
+      {
+      case TRIG:
+      case TRIG6:
+        outfile << 3 << " " << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << "\n";
+        types.push_back(5);
+        break;
+      case QUAD:
+        outfile << 4 << " " << el[0]-1 << " " << el[1]-1 << " " << el[2]-1 << " " << el[3]-1 << "\n";
+        types.push_back(9);
+        break;
+      default:
+        throw ngcore::Exception("Unexpected element type");
+      break;
+      }
     }
   }
 
@@ -117,6 +121,14 @@ void WriteVtkFormat (const Mesh & mesh,
     outfile << type_id << "\n";
   }
 
+  outfile << "CELL_DATA " << domains.size() << "\n";
+  outfile << "SCALARS scalars int 1\n";
+  outfile << "LOOKUP_TABLE default\n";
+  for (auto id: domains)
+  {
+    outfile << id << "\n";
+  }
+
   outfile.close();
 }
 }
diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp
index df3447a0..18dfd24b 100644
--- a/nglib/nglib.cpp
+++ b/nglib/nglib.cpp
@@ -591,6 +591,68 @@ namespace nglib
    }
 
 
+   DLL_HEADER Ng_Geometry_2D * Ng_NewGeometry_2D ()
+   {
+      SplineGeometry2d * geom = new SplineGeometry2d();
+      return (Ng_Geometry_2D *)geom;
+   }
+
+   DLL_HEADER void Ng_DeleteGeometry_2D (Ng_Geometry_2D * geom)
+   {
+      if (geom)
+      {
+         SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
+         delete spline_geom;
+         geom = NULL;
+      }
+   }
+
+   DLL_HEADER void Ng_AppendPoint_2D (Ng_Geometry_2D* geom, double * x, double h)
+   {
+      if (geom)
+      {
+         SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
+         Point<2> p(x[0],x[1]);
+         spline_geom->AppendPoint(p, h);
+      }
+   }
+
+   DLL_HEADER void Ng_AppendLineSegment_2D (Ng_Geometry_2D* geom, int n1, int n2,
+      int leftdomain, int rightdomain, double h)
+   {
+      if (geom)
+      {
+         SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
+         // zero-offset!
+         LineSeg<2>* line = new LineSeg<2>(spline_geom->geompoints[n1-1], spline_geom->geompoints[n2-1]);
+         SplineSegExt* seg = new SplineSegExt(*line);
+         seg->leftdom = leftdomain;
+         seg->rightdom = rightdomain;
+         seg->hmax = h;
+         seg->copyfrom = -1;
+         seg->bc = 1;
+         spline_geom->AppendSegment(seg);
+      }
+   }
+
+   DLL_HEADER void Ng_AppendSplinSegment_2D (Ng_Geometry_2D* geom, int n1, int n2, int n3,
+      int leftdomain, int rightdomain, double h)
+   {
+      if (geom)
+      {
+         SplineGeometry2d* spline_geom = (SplineGeometry2d*)geom;
+         // zero-offset!
+         SplineSeg3<2>* line = new SplineSeg3<2>(spline_geom->geompoints[n1-1], spline_geom->geompoints[n2-1], spline_geom->geompoints[n3-1]);
+         SplineSegExt* seg = new SplineSegExt(*line);
+         seg->leftdom = leftdomain;
+         seg->rightdom = rightdomain;
+         seg->hmax = h;
+         seg->copyfrom = -1;
+         seg->bc = 1;
+         spline_geom->AppendSegment(seg);
+      }
+   }
+
    DLL_HEADER Ng_Result Ng_GenerateMesh_2D (Ng_Geometry_2D * geom,
                                             Ng_Mesh ** mesh,
                                             Ng_Meshing_Parameters * mp)
@@ -601,9 +663,6 @@ namespace nglib
 
       shared_ptr<Mesh> m(new Mesh, &NOOP_Deleter);
       MeshFromSpline2D (*(SplineGeometry2d*)geom, m, mparam);
-      // new shared_ptr<Mesh> (m);  // hack to keep mesh m alive 
-      
-      cout << m->GetNSE() << " elements, " << m->GetNP() << " points" << endl;
 
       *mesh = (Ng_Mesh*)m.get();
       return NG_OK;
diff --git a/nglib/nglib.h b/nglib/nglib.h
index 38e60faa..28155ad7 100644
--- a/nglib/nglib.h
+++ b/nglib/nglib.h
@@ -625,6 +625,9 @@ DLL_HEADER void Ng_AppendPoint_2D(Ng_Geometry_2D* geom, double * x, double h);
 DLL_HEADER void Ng_AppendLineSegment_2D(Ng_Geometry_2D* geom, int pi1, int pi2,
 	int leftdomain, int rightdomain, double h);
 
+DLL_HEADER void Ng_AppendSplinSegment_2D(Ng_Geometry_2D* geom, int pi1, int pi2, int pi3,
+	int leftdomain, int rightdomain, double h);
+
 DLL_HEADER void Ng_AddPoint_2D (Ng_Mesh * mesh, double * x);
 DLL_HEADER void Ng_AddBoundarySeg_2D (Ng_Mesh * mesh, int pi1, int pi2, int domain_in = -1, int domain_out = -1);