diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp
index 5c43e23e..c85480da 100644
--- a/libsrc/csg/python_csg.cpp
+++ b/libsrc/csg/python_csg.cpp
@@ -441,6 +441,9 @@ DLL_HEADER void ExportCSG(py::module &m)
 	    }
 	    spsol->AddSurfaces(self);
 	    int tlonr = self.SetTopLevelObject(spsol->GetSolid(), surf.get());
+	    self.GetTopLevelObject(tlonr) -> SetBCProp(surf->GetBase()->GetBCProperty());
+	    self.GetTopLevelObject(tlonr) -> SetBCName(surf->GetBase()->GetBCName());
+	    self.GetTopLevelObject(tlonr) -> SetMaxH(surf->GetBase()->GetMaxH());
 	    for(auto p : surf->GetPoints())
 		self.AddUserPoint(p);
             self.AddSplineSurface(surf);
diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp
index d95f71e1..550ce1a8 100644
--- a/libsrc/interface/nginterface_v2.cpp
+++ b/libsrc/interface/nginterface_v2.cpp
@@ -757,7 +757,7 @@ namespace netgen
     double hdxdxi[4][3];
     for (int j = 0; j<4;j++)
       hxi[j][0] = ((double*)&(xi[0]))[j];
-    MultiElementTransformation<1,3> (elnr, 4, &hxi[0][0], 1, &hx[0][0], 2, &hdxdxi[0][0],4);
+    MultiElementTransformation<1,3> (elnr, 4, &hxi[0][0], 1, &hx[0][0], 3, &hdxdxi[0][0],3);
     for(int j=0; j<4; j++)
       for(int k=0; k<3; k++)
 	((double*)&(x[k]))[j] = hx[j][k];
diff --git a/py_tutorials/merge.py b/py_tutorials/merge.py
index 931d095a..abacd38e 100644
--- a/py_tutorials/merge.py
+++ b/py_tutorials/merge.py
@@ -3,28 +3,35 @@ from netgen.csg import *
 
 from ngsolve import ngsglobals
 ngsglobals.msg_level = 2
-
+# generate brick and mesh it
 geo1 = CSGeometry()
 geo1.Add (OrthoBrick( Pnt(0,0,0), Pnt(1,1,1) ))
 m1 = geo1.GenerateMesh (maxh=0.1)
 m1.Refine()
-m1.Refine()
-
 
+# generate sphere and mesh it
 geo2 = CSGeometry()
 geo2.Add (Sphere (Pnt(0.5,0.5,0.5), 0.1))
 m2 = geo2.GenerateMesh (maxh=0.05)
 m2.Refine()
 m2.Refine()
 
-print ("********************")
-print ("** start merging  **")
-print ("********************")
+print ("***************************")
+print ("** merging suface meshes **")
+print ("***************************")
 
+# create an empty mesh
 mesh = Mesh()
-fd_outside = mesh.Add (FaceDescriptor(surfnr=1,domin=1))
-fd_inside = mesh.Add (FaceDescriptor(surfnr=2,domin=2,domout=1))
 
+# a face-descriptor stores properties associated with a set of surface elements
+# bc .. boundary condition marker,
+# domin/domout .. domain-number in front/back of surface elements (0 = void),
+# surfnr .. number of the surface described by the face-descriptor
+
+fd_outside = mesh.Add (FaceDescriptor(bc=1,domin=1,surfnr=1))
+fd_inside = mesh.Add (FaceDescriptor(bc=2,domin=2,domout=1,surfnr=2))
+# copy all boundary points from first mesh to new mesh.
+# pmap1 maps point-numbers from old to new mesh
 
 pmap1 = { }
 for e in m1.Elements2D():
@@ -32,11 +39,17 @@ for e in m1.Elements2D():
         if (v not in pmap1):
             pmap1[v] = mesh.Add (m1[v])
 
+
+# copy surface elements from first mesh to new mesh
+# we have to map point-numbers:
+
 for e in m1.Elements2D():
     mesh.Add (Element2D (fd_outside, [pmap1[v] for v in e.vertices]))
 
 
 
+# same for the second mesh:
+
 pmap2 = { }
 for e in m2.Elements2D():
     for v in e.vertices: