diff --git a/CMakeLists.txt b/CMakeLists.txt
index a40c89c5..db13de81 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -35,8 +35,6 @@ option( USE_SUPERBUILD       "use ccache" ON)
 
 set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_modules")
 
-set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
-
 if(APPLE)
   set(INSTALL_DIR_DEFAULT /Applications/Netgen.app)
 else(APPLE)
@@ -78,6 +76,8 @@ else()
   endif()
 endif()
 
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
 include (${CMAKE_CURRENT_LIST_DIR}/cmake/generate_version_file.cmake)
 set(CPACK_PACKAGE_VERSION "${NETGEN_VERSION}")
 
@@ -191,7 +191,6 @@ check_include_files (dlfcn.h HAVE_DLFCN_H)
 if(HAVE_DLFCN_H)
   add_definitions(-DHAVE_DLFCN_H)
 endif()
-add_definitions(-DNETGEN_VERSION="${NETGEN_VERSION}")
 
 include_directories(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 
diff --git a/cmake/cmake_modules/FindOpenCasCade.cmake b/cmake/cmake_modules/FindOpenCasCade.cmake
index e82f8dac..870e8ed5 100644
--- a/cmake/cmake_modules/FindOpenCasCade.cmake
+++ b/cmake/cmake_modules/FindOpenCasCade.cmake
@@ -27,7 +27,7 @@ endif(WIN32)
 
 if(OCC_LIBRARY AND NOT OCC_LIBRARY_DIR)
     get_filename_component(OCC_LIBRARY_DIR ${OCC_LIBRARY} PATH)
-endif(OCC_LIBRARY)
+endif(OCC_LIBRARY AND NOT OCC_LIBRARY_DIR)
 
 if(OCC_INCLUDE_DIR)
     file(STRINGS ${OCC_INCLUDE_DIR}/Standard_Version.hxx OCC_MAJOR
diff --git a/libsrc/core/array.hpp b/libsrc/core/array.hpp
index ccafef6b..a1161001 100644
--- a/libsrc/core/array.hpp
+++ b/libsrc/core/array.hpp
@@ -1284,7 +1284,7 @@ namespace ngcore
 
   /// bubble sort array
   template <class T, class S>
-  inline void BubbleSort (FlatArray<T> data, FlatArray<S> slave)
+  inline void BubbleSort (FlatArray<T> data, FlatArray<S> index)
   {
     for (size_t i = 0; i < data.Size(); i++)
       for (size_t j = i+1; j < data.Size(); j++)
@@ -1294,9 +1294,9 @@ namespace ngcore
 	    data[i] = data[j];
 	    data[j] = hv;
 
-	    S hvs = slave[i];
-	    slave[i] = slave[j];
-	    slave[j] = hvs;
+	    S hvs = index[i];
+	    index[i] = index[j];
+	    index[j] = hvs;
 	  }
   }
 
diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp
index 4874f0a9..cc3e4464 100644
--- a/libsrc/csg/csgeom.cpp
+++ b/libsrc/csg/csgeom.cpp
@@ -130,6 +130,7 @@ namespace netgen
                Point<3> & newp, EdgePointGeomInfo & newgi) const
   {
     Point<3> hnewp = p1+secpoint*(p2-p1);
+
     //(*testout) << "hnewp " << hnewp << " s1 " << surfi1 << " s2 " << surfi2 << endl;
     if (surfi1 != -1 && surfi2 != -1 && surfi1 != surfi2)
       {
diff --git a/libsrc/general/ngarray.hpp b/libsrc/general/ngarray.hpp
index f3eea9d4..fa160a8b 100644
--- a/libsrc/general/ngarray.hpp
+++ b/libsrc/general/ngarray.hpp
@@ -730,7 +730,7 @@ namespace netgen
 
   /// bubble sort array
   template <class T, class S>
-  inline void BubbleSort (NgFlatArray<T> & data, NgFlatArray<S> & slave)
+  inline void BubbleSort (NgFlatArray<T> & data, NgFlatArray<S> & index)
   {
     for (int i = 0; i < data.Size(); i++)
       for (int j = i+1; j < data.Size(); j++)
@@ -740,16 +740,16 @@ namespace netgen
 	    data[i] = data[j];
 	    data[j] = hv;
 	    
-	    S hvs = slave[i];
-	    slave[i] = slave[j];
-	    slave[j] = hvs;
+	    S hvs = index[i];
+	    index[i] = index[j];
+	    index[j] = hvs;
 	  }
   }
 
 
   template <class T, class S>
   void QuickSortRec (NgFlatArray<T> & data,
-		     NgFlatArray<S> & slave,
+		     NgFlatArray<S> & index,
 		     int left, int right)
   {
     int i = left;
@@ -764,20 +764,20 @@ namespace netgen
 	if (i <= j)
 	  {
             ngcore::Swap (data[i], data[j]);
-            ngcore::Swap (slave[i], slave[j]);
+            ngcore::Swap (index[i], index[j]);
 	    i++; j--;
 	  }
       }
     while (i <= j);
-    if (left < j) QuickSortRec (data, slave, left, j);
-    if (i < right) QuickSortRec (data, slave, i, right);
+    if (left < j) QuickSortRec (data, index, left, j);
+    if (i < right) QuickSortRec (data, index, i, right);
   }
 
   template <class T, class S>
-  void QuickSort (NgFlatArray<T> & data, NgFlatArray<S> & slave)
+  void QuickSort (NgFlatArray<T> & data, NgFlatArray<S> & index)
   {
     if (data.Size() > 1)
-      QuickSortRec (data, slave, 0, data.Size()-1);
+      QuickSortRec (data, index, 0, data.Size()-1);
   }
 
 
diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp
index 09ef0260..c1eadf1a 100644
--- a/libsrc/geom2d/genmesh2d.cpp
+++ b/libsrc/geom2d/genmesh2d.cpp
@@ -527,11 +527,16 @@ namespace netgen
 
 
           for (PointIndex pix = nextpi[c1], ix = 0; pix != c2; pix = nextpi[pix], ix++)
+          {
+            Point<3> px = (*mesh)[pix];
             for (PointIndex piy = nextpi[c2], iy = 0; piy != c3; piy = nextpi[piy], iy++)
               {
-                Point<3> p = (*mesh)[pix] + ( (*mesh)[piy] - (*mesh)[c2] );
-                pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (p , 1, FIXEDPOINT);
+                double lam = Dist((*mesh)[piy],(*mesh)[c2]) / Dist((*mesh)[c3],(*mesh)[c2]);
+                auto pix1 = pts[(nex+1)*ney+ix+1];
+                auto pnew = px + lam*((*mesh)[pix1]-px);
+                pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (pnew, 1, FIXEDPOINT);
               }
+          }
 
           for (int i = 0; i < ney; i++)
             for (int j = 0; j < nex; j++)
@@ -545,6 +550,10 @@ namespace netgen
 
                 mesh -> AddSurfaceElement (el);
               }
+          char* material;
+          geometry.GetMaterial(domnr, material);
+          if(material)
+            mesh->SetMaterial(domnr, material);
         }
 
 
diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp
index ff4c459c..f5841b19 100644
--- a/libsrc/geom2d/geometry2d.hpp
+++ b/libsrc/geom2d/geometry2d.hpp
@@ -133,7 +133,7 @@ namespace netgen
     NgArray<char*> materials;
     NgArray<double> maxh;
     NgArray<bool> quadmeshing;
-    NgArray<bool> tensormeshing;
+    Array<bool> tensormeshing;
     NgArray<int> layer;
     NgArray<string*> bcnames;
     double elto0 = 1.0;
@@ -216,9 +216,20 @@ namespace netgen
     }
     bool GetDomainTensorMeshing ( int domnr ) 
     { 
-      if ( tensormeshing.Size() ) return tensormeshing[domnr-1]; 
+      if ( tensormeshing.Size()>=domnr ) return tensormeshing[domnr-1];
       else return false;
     }
+    void SetDomainTensorMeshing ( int domnr, bool tm )
+    {
+      if ( tensormeshing.Size()<domnr )
+      {
+        auto oldsize = tensormeshing.Size();
+        tensormeshing.SetSize(domnr);
+        for(auto i : IntRange(oldsize, domnr-1))
+          tensormeshing[i] = false;
+      }
+      tensormeshing[domnr-1] = tm;
+    }
     int GetDomainLayer ( int domnr ) 
     { 
       if ( layer.Size() ) return layer[domnr-1]; 
diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp
index a890579b..f591e74e 100644
--- a/libsrc/geom2d/python_geom2d.cpp
+++ b/libsrc/geom2d/python_geom2d.cpp
@@ -392,6 +392,7 @@ DLL_HEADER void ExportGeom2d(py::module &m)
                 }, py::arg("mp") = nullptr,
       py::call_guard<py::gil_scoped_release>(),
       meshingparameter_description.c_str())
+    .def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing)
     ;
   
 }
diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp
index 6d6d1b85..20fccd2a 100644
--- a/libsrc/include/nginterface_v2_impl.hpp
+++ b/libsrc/include/nginterface_v2_impl.hpp
@@ -111,7 +111,13 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const
   ret.faces.num = 0;
   ret.faces.ptr = NULL;
 
-  if (mesh->GetDimension() == 2)
+  if (mesh->GetDimension() == 3)
+    {
+      ret.facets.num = 0;
+      ret.facets.base = 0;
+      ret.facets.ptr = nullptr;
+    }
+  else if (mesh->GetDimension() == 2)
     {
       ret.facets.num = 1;
       ret.facets.base = 0;
diff --git a/libsrc/interface/readtetmesh.cpp b/libsrc/interface/readtetmesh.cpp
index e061c722..46607c59 100644
--- a/libsrc/interface/readtetmesh.cpp
+++ b/libsrc/interface/readtetmesh.cpp
@@ -154,7 +154,7 @@ namespace netgen
             break;
 
           case 7:
-            // NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), PID:
+            // NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), PID:
             {
               cout << "read nodes" << endl;
               for(int i=0; i<nnodes; i++)
@@ -175,7 +175,7 @@ namespace netgen
             break;
 
           case 9:
-            // MasterNodeID, SlaveNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
+            // MasterNodeID, MinionNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
             for(int i=0; i<nperiodicmasternodes; i++)
               {
                 for(int j=0; j<2; j++)
@@ -191,7 +191,7 @@ namespace netgen
             break;
 
           case 11:
-            // MasterNodeID, 3-SlaveNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2)
+            // MasterNodeID, 3-MinionNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2)
             for(int i=0; i<ncornerperiodicmasternodes; i++)
               {
                 for(int j=0; j<4; j++)
@@ -208,7 +208,7 @@ namespace netgen
             break;
 
           case 13:
-            //MasterNodeID, 7-SlaveNodeID's, TranslCodes
+            //MasterNodeID, 7-MinionNodeID's, TranslCodes
             for(int i=0; i<ncubicperiodicmasternodes; i++)
               {
                 for(int j=0; j<8; j++)
@@ -220,7 +220,7 @@ namespace netgen
             break;
 
           case 14:
-            // EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), PID
+            // EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), PID
             cout << "read edges" << endl;
             // nullstarted = false;
             segmentdata.SetSize(nedges);
@@ -243,7 +243,7 @@ namespace netgen
             break;
 
           case 16:
-            // MasterEdgeID, SlaveEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
+            // MasterEdgeID, MinionEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2)
             for(int i=0; i<nperiodicmasteredges; i++)
               in >> dummyint >> dummyint >> dummyint;
             break;
@@ -254,7 +254,7 @@ namespace netgen
             break;
 
           case 18:
-            // MasterEdgeID, 3 SlaveEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2)
+            // MasterEdgeID, 3 MinionEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2)
             for(int i=0; i<ncornerperiodicmasteredges; i++)
               {
                 in >> dummyint;
@@ -266,7 +266,7 @@ namespace netgen
             break;
 
           case 19:
-            // FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PSlave), PID
+            // FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PMinion), PID
             {
               //Segment seg;
               int segnum_ng[3];
@@ -343,7 +343,7 @@ namespace netgen
             break;
 
           case 21:
-            // MasterFaceID, SlaveFaceID, TranslCode (1=dS1 2=dS2)
+            // MasterFaceID, MinionFaceID, TranslCode (1=dS1 2=dS2)
             {
               Vec<3> randomvec(-1.32834,3.82399,0.5429151);
               int maxtransl = -1;
diff --git a/libsrc/interface/writeabaqus.cpp b/libsrc/interface/writeabaqus.cpp
index c1582d0f..12a1eddd 100644
--- a/libsrc/interface/writeabaqus.cpp
+++ b/libsrc/interface/writeabaqus.cpp
@@ -268,25 +268,25 @@ void WriteAbaqusFormat (const Mesh & mesh,
 
         cout << "masternode = " << masternode << " = "
             << mesh.Point(masternode) << endl;
-        NgArray<int> slaves(3);
+        NgArray<int> minions(3);
         for (i = 1; i <= 3; i++)
         {
             mesh.GetIdentifications().GetPairs (i, pairs);
             for (j = 1; j <= pairs.Size(); j++)
             {
                 if (pairs.Get(j).I1() == masternode)
-                    slaves.Elem(i) = pairs.Get(j).I2();
+                    minions.Elem(i) = pairs.Get(j).I2();
             }
-            cout << "slave(" << i << ") = " << slaves.Get(i)
-                << " = " << mesh.Point(slaves.Get(i)) << endl;
+            cout << "minion(" << i << ") = " << minions.Get(i)
+                << " = " << mesh.Point(minions.Get(i)) << endl;
         }
 
 
         outfile << "**\n"
             << "*NSET,NSET=CTENODS\n"
-            << slaves.Get(1) << ", " 
-            << slaves.Get(2) << ", " 
-            << slaves.Get(3) << endl;
+            << minions.Get(1) << ", " 
+            << minions.Get(2) << ", " 
+            << minions.Get(3) << endl;
 
 
         outfile << "**\n"
@@ -300,7 +300,7 @@ void WriteAbaqusFormat (const Mesh & mesh,
             << "*BOUNDARY, OP=NEW\n";
         for (j = 1; j <= 3; j++)
         {
-            Vec3d v(mesh.Point(masternode), mesh.Point(slaves.Get(j)));
+            Vec3d v(mesh.Point(masternode), mesh.Point(minions.Get(j)));
             double vlen = v.Length();
             int dir = 0;
             if (fabs (v.X()) > 0.9 * vlen) dir = 2;
@@ -308,7 +308,7 @@ void WriteAbaqusFormat (const Mesh & mesh,
             if (fabs (v.Z()) > 0.9 * vlen) dir = 1;
             if (!dir)
                 cout << "ERROR: Problem with rigid body constraints" << endl;
-            outfile << slaves.Get(j) << ", " << dir << ",,    0.\n";
+            outfile << minions.Get(j) << ", " << dir << ",,    0.\n";
         }
 
         outfile << "**\n"
@@ -333,14 +333,13 @@ void WriteAbaqusFormat (const Mesh & mesh,
                         mpc << "4" << "\n";
                         mpc << pairs.Get(j).I2() << "," << k << ", -1.0, ";
                         mpc << pairs.Get(j).I1() << "," << k << ", 1.0, ";
-                        mpc << slaves.Get(i) << "," << k << ", 1.0, ";
+                        mpc << minions.Get(i) << "," << k << ", 1.0, ";
                         mpc << masternode << "," << k << ", -1.0 \n";
                     }
                 }
         }
     }
 
-
     cout << "done" << endl;
 }
 
diff --git a/libsrc/interface/writefeap.cpp b/libsrc/interface/writefeap.cpp
index dc2574a9..84f71eef 100644
--- a/libsrc/interface/writefeap.cpp
+++ b/libsrc/interface/writefeap.cpp
@@ -144,11 +144,11 @@ void WriteFEAPFormat (const Mesh & mesh,
 
   // BEGIN CONTACT OUTPUT
   /*      
-	  int masterindex, slaveindex;
+	  int masterindex, minionindex;
 	  cout << "Master Surface index = ";
 	  cin >> masterindex;
-	  cout << "Slave Surface index  = ";
-	  cin >> slaveindex;
+	  cout << "Minion Surface index  = ";
+	  cin >> minionindex;
 
 
 	  // CONTACT SURFACE 1
@@ -196,7 +196,7 @@ void WriteFEAPFormat (const Mesh & mesh,
 	  Element2d sel = mesh.SurfaceElement(i);
 	  if (invertsurf)
 	  sel.Invert();
-	  if (mesh.GetFaceDescriptor(sel.GetIndex ()).BCProperty() == slaveindex)
+	  if (mesh.GetFaceDescriptor(sel.GetIndex ()).BCProperty() == minionindex)
 	  {
 	  zz++;
 	  outfile.width(14);
diff --git a/libsrc/interface/writetet.cpp b/libsrc/interface/writetet.cpp
index 3bc0ba5a..2dbfa440 100644
--- a/libsrc/interface/writetet.cpp
+++ b/libsrc/interface/writetet.cpp
@@ -427,7 +427,7 @@ namespace netgen
 	    << numedges << " " 
 	    << numnodes << endl << endl;
 
-    outfile << "// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), "<< uidpid <<":\n" \
+    outfile << "// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), "<< uidpid <<":\n" \
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
        
 
@@ -515,7 +515,7 @@ namespace netgen
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n" \
 	    << n2 << "\n"			       \
 	    << "\n" \
-	    << "// MasterNodeID, SlaveNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n" \
+	    << "// MasterNodeID, MinionNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n" \
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
     for(int i=0; i<id_groups.Size(); i++)
       {
@@ -538,7 +538,7 @@ namespace netgen
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n" \
 	    << n4 << "\n"				      \
 	    << "\n" \
-	    << "// MasterNodeID, 3-SlaveNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2):\n" \
+	    << "// MasterNodeID, 3-MinionNodeID's, 3-TranslCodes (1=dS1 2=dS2 3=dS1+dS2):\n" \
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
 
 
@@ -565,7 +565,7 @@ namespace netgen
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n" \
 	    << n8 << "\n"				     \
 	    << "\n" \
-	    << "// MasterNodeID, 7-SlaveNodeID's, TranslCodes:\n" \
+	    << "// MasterNodeID, 7-MinionNodeID's, TranslCodes:\n" \
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
     for(int i=0; i<id_groups.Size(); i++)
       {
@@ -586,7 +586,7 @@ namespace netgen
     
 
 
-    outfile << "// EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), "<<uidpid<<":\n" \
+    outfile << "// EdgeID, NodeID0, NodeID1, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), "<<uidpid<<":\n" \
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
 
     
@@ -760,7 +760,7 @@ namespace netgen
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"\
 	    << n2 << "\n"			      \
 	    << "\n"\
-	    << "// MasterEdgeID, SlaveEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
+	    << "// MasterEdgeID, MinionEdgeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
     for(int i=0; i<id_groups.Size(); i++)
       {
@@ -782,7 +782,7 @@ namespace netgen
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"\
 	    << n4 << "\n"				     \
 	    << "\n"\
-	    << "// MasterEdgeID, 3 SlaveEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
+	    << "// MasterEdgeID, 3 MinionEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n"\
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
     for(int i=0; i<id_groups.Size(); i++)
       {
@@ -801,7 +801,7 @@ namespace netgen
     outfile << endl;
 
 
-    outfile << "// FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PSlave), "<<uidpid<<":\n" \
+    outfile << "// FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PMinion), "<<uidpid<<":\n" \
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
 
     
@@ -921,7 +921,7 @@ namespace netgen
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"\
 	    << n2 << "\n"			      \
 	    << "\n"\
-	    << "// MasterFaceID, SlaveFaceID, TranslCode (1=dS1 2=dS2):\n"\
+	    << "// MasterFaceID, MinionFaceID, TranslCode (1=dS1 2=dS2):\n"\
 	    << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";
     for(int i=0; i<id_groups.Size(); i++)
       {
diff --git a/libsrc/meshing/CMakeLists.txt b/libsrc/meshing/CMakeLists.txt
index 47b99e79..9bf45a89 100644
--- a/libsrc/meshing/CMakeLists.txt
+++ b/libsrc/meshing/CMakeLists.txt
@@ -12,7 +12,7 @@ add_library(mesh ${NG_LIB_TYPE}
         smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp		    
         topology.cpp triarls.cpp validate.cpp bcfunctions.cpp   
         parallelmesh.cpp  paralleltop.cpp  paralleltop.hpp basegeom.cpp     
-        python_mesh.cpp hexarls.cpp
+        python_mesh.cpp hexarls.cpp surfacegeom.cpp
         ../../ng/onetcl.cpp
         ${mesh_object_libs}
         )
@@ -37,6 +37,6 @@ install(FILES
   localh.hpp meshclass.hpp meshfunc.hpp meshing2.hpp meshing3.hpp
   meshing.hpp meshtool.hpp meshtype.hpp msghandler.hpp paralleltop.hpp
   ruler2.hpp ruler3.hpp specials.hpp topology.hpp validate.hpp
-  python_mesh.hpp
+  python_mesh.hpp surfacegeom.hpp
   DESTINATION ${NG_INSTALL_DIR_INCLUDE}/meshing COMPONENT netgen_devel
 )
diff --git a/libsrc/meshing/bcfunctions.cpp b/libsrc/meshing/bcfunctions.cpp
index d430a3a7..a7214954 100644
--- a/libsrc/meshing/bcfunctions.cpp
+++ b/libsrc/meshing/bcfunctions.cpp
@@ -322,7 +322,7 @@ namespace netgen
       colours_sorted.SetSize(all_colours.Size()+1);
       faces_sorted = 0;
       
-      // Slave NgArray to identify the colours the faces were assigned to, 
+      // Index NgArray to identify the colours the faces were assigned to, 
       // after the bubble sort routine to sort the automatic boundary 
       // identifiers according to the number of surface mesh elements 
       // of a given colour
diff --git a/libsrc/meshing/bcfunctions.hpp b/libsrc/meshing/bcfunctions.hpp
index 55aca571..967dd0d8 100644
--- a/libsrc/meshing/bcfunctions.hpp
+++ b/libsrc/meshing/bcfunctions.hpp
@@ -26,7 +26,7 @@ namespace netgen
           - Use colour index 0 (zero) for all faces with no colour defined
           - Calculate the number of faces of the surface mesh for each colour
           - Sort the number of surface elements in ascending order, with the 
-            colour indices as a slave
+            colour indices as a index
           - Use the indices of the sorted array as the BC property number
 
           Example: If there are 3 colours, present in the mesh and the number 
diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index ae7c7a30..3ce5891a 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -159,10 +159,13 @@ namespace netgen
                   auto& old_fd = mesh.GetFaceDescriptor(si);
                   int domout = blp.outside ? old_fd.DomainOut() : blp.new_matnrs[layer-1];
                   int domin = blp.outside ? blp.new_matnrs[layer-1] : old_fd.DomainIn();
-                  FaceDescriptor fd(max_surface_index-1,
+                  // -1 surf nr is so that curving does not do anything
+                  FaceDescriptor fd(-1,
                                     domin, domout, -1);
                   fd.SetBCProperty(max_surface_index);
                   mesh.AddFaceDescriptor(fd);
+                  mesh.SetBCName(max_surface_index-1,
+                                 "mapped_" + old_fd.GetBCName());
                   return max_surface_index;
                 }
               return last_layer_surface_index_map[si];
@@ -269,55 +272,53 @@ namespace netgen
           if(blp.grow_edges)
             for(SegmentIndex sei = 0; sei < nseg; sei++)
               {
-                PointIndex seg_p1 = mesh[sei][0];
-                PointIndex seg_p2 = mesh[sei][1];
+                auto& segi = mesh[sei];
 
                 // Only go in if the segment is still active, and if both its
                 // surface index is part of the "hit-list"
                 if(segsel.Test(sei))
                   {
-                    if(blp.surfid.Contains(mesh[sei].si))
-                  {
-                    // clear the bit to indicate that this segment has been processed
-                    segsel.Clear(sei);
-
-                    // Find matching segment pair on other surface
-                    for(SegmentIndex sej = 0; sej < nseg; sej++)
+                    if(blp.surfid.Contains(segi.si))
                       {
-                        PointIndex segpair_p1 = mesh[sej][1];
-                        PointIndex segpair_p2 = mesh[sej][0];
+                        // clear the bit to indicate that this segment has been processed
+                        segsel.Clear(sei);
 
-                        // Find the segment pair on the neighbouring surface element
-                        // Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0]
-                        if(segsel.Test(sej) && ((segpair_p1 == seg_p1) && (segpair_p2 == seg_p2)))
+                        // Find matching segment pair on other surface
+                        for(SegmentIndex sej = 0; sej < nseg; sej++)
                           {
-                            // clear bit to indicate that processing of this segment is done
-                            segsel.Clear(sej);
+                            auto& segj = mesh[sej];
+                            // Find the segment pair on the neighbouring surface element
+                            // Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0]
+                            if(segsel.Test(sej) && ((segi[0] == segj[1]) && (segi[1] == segj[0])))
+                              {
+                                // clear bit to indicate that processing of this segment is done
+                                segsel.Clear(sej);
 
-                            // Only worry about those surfaces which are not in the
-                            // boundary layer list
-                            if(!blp.surfid.Contains(mesh[sej].si))
+                                // if segj is not in surfel list we nned to add quads
+                            if(!blp.surfid.Contains(segj.si))
                               {
                                 SurfaceElementIndex pnt_commelem;
                                 SetInvalid(pnt_commelem);
 
-                                auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segpair_p1);
-                                auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2);
+                                auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segj[0]);
+                                auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segj[1]);
 
                                 for(auto pnt1_sei : pnt1_elems)
-                                  if(mesh[pnt1_sei].GetIndex() == mesh[sej].si)
+                                  if(mesh[pnt1_sei].GetIndex() == segj.si)
                                     for(auto pnt2_sei : pnt2_elems)
                                       if(pnt1_sei == pnt2_sei)
                                         pnt_commelem = pnt1_sei;
 
                                 if(IsInvalid(pnt_commelem))
-                                  throw Exception("Couldn't find element on other side for " + ToString(segpair_p1) + " to " + ToString(segpair_p2));
+                                  throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1]));
 
                                 const auto& commsel = mesh[pnt_commelem];
                                 auto surfelem_vect = GetSurfaceNormal(mesh, commsel);
                                 if(blp.outside)
                                   surfelem_vect *= -1;
                                 Element2d sel(QUAD);
+                                auto seg_p1 = segi[0];
+                                auto seg_p2 = segi[1];
                                 if(blp.outside)
                                   Swap(seg_p1, seg_p2);
                                 sel[0] = seg_p1;
@@ -330,7 +331,7 @@ namespace netgen
                                   {
                                     domains_to_surf_index[domains] = ++max_surface_index;
                                     domains_to_surf_index[make_tuple(max_surface_index, get<1>(domains), get<2>(domains))] = max_surface_index;
-                                    FaceDescriptor fd(max_surface_index-1,
+                                    FaceDescriptor fd(-1,
                                                       get<1>(domains),
                                                       get<2>(domains),
                                                       -1);
@@ -365,46 +366,37 @@ namespace netgen
                                 seg_2.edgenr = pi_to_edgenr[points];
                                 seg_2.si = new_index;
                                 mesh.AddSegment(seg_2);
-                                mesh[sej].si = new_index;
                               }
 
                             // in last layer insert new segments
                             if(layer == blp.heights.Size())
                               {
-                                Segment s1 = mesh[sei];
-                                Segment s2 = mesh[sej];
+                                Segment s1 = segi;
+                                Segment s2 = segj;
                                 s1.edgenr = ++max_edge_nr;
                                 s2.edgenr = max_edge_nr;
-                                bool create_it = true;
-                                if(blp.surfid.Contains(mesh[sej].si))
-                                  {
-                                    if(last_layer_surface_index_map.find(s1.si) != last_layer_surface_index_map.end() &&
-                                       last_layer_surface_index_map.find(s2.si) != last_layer_surface_index_map.end())
-                                      // edge already mapped
-                                      create_it = false;
-                                    s2.si = map_surface_index(s2.si);
-                                  }
+                                if(blp.surfid.Contains(segj.si))
+                                  s2.si = map_surface_index(segj.si);
                                 else
                                   {
-                                    s2.si = domains_to_surf_index[make_tuple(s2.si,
-                                                                             blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
+                                    auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())];
+                                    if(blp.outside)
+                                      s2.si = side_surf;
+                                    else
+                                      segj.si = side_surf;
                                   }
                                 s1.si = map_surface_index(s1.si);
-                                if(create_it)
-                                  {
-                                    mesh.AddSegment(s1);
-                                    mesh.AddSegment(s2);
-                                  }
+                                s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1;
+                                mesh.AddSegment(s1);
+                                mesh.AddSegment(s2);
+                              }
+                            // segi[0] = mapto[segi[0]] not working somehow?
+                            mesh[sei][0] = mapto[segi[0]];
+                            mesh[sei][1] = mapto[segi[1]];
+                            mesh[sej][0] = mapto[segj[0]];
+                            mesh[sej][1] = mapto[segj[1]];
                               }
-
-                            // remap the segments to the new points
-                            mesh[sei][0] = mapto[mesh[sei][0]];
-                            mesh[sei][1] = mapto[mesh[sei][1]];
-                            mesh[sej][1] = mapto[mesh[sej][1]];
-                            mesh[sej][0] = mapto[mesh[sej][0]];
-
                           }
-                      }
                   }
                 else
                   {
@@ -455,9 +447,9 @@ namespace netgen
             {
               for(SurfaceElementIndex si = 0; si < nse; si++)
                 {
-                  if(blp.surfid.Contains(mesh[si].GetIndex()))
+                  const auto& sel = mesh[si];
+                  if(blp.surfid.Contains(sel.GetIndex()))
                     {
-                      const auto& sel = mesh[si];
                       Element2d newel = sel;
                       newel.SetIndex(map_surface_index(sel.GetIndex()));
                       mesh.AddSurfaceElement(newel);
@@ -546,36 +538,28 @@ namespace netgen
           for(SurfaceElementIndex sei : Range(nse))
             {
               auto& sel = mesh[sei];
-              bool to_move = blp.surfid.Contains(sel.GetIndex());
-              if(blp.domains.Size())
+              if(!blp.surfid.Contains(sel.GetIndex()))
                 {
-                  if(blp.outside)
-                    to_move |= blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()];
-                  else
-                    to_move |= !blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()];
-                }
-
-              if(to_move)
-                {
-                  for(auto& pnum : sel.PNums())
-                    // Check (Doublecheck) if the corresponding point has a
-                    // copy available for remapping
-                    if(mapto[pnum].IsValid())
-                      // Map the surface elements to the new points
-                      pnum = mapto[pnum];
+                  const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex());
+                  if(blp.outside &&
+                     (!blp.domains[fd.DomainIn()] && !blp.domains[fd.DomainOut()]))
+                    continue;
+                  if(!blp.outside &&
+                     (blp.domains[fd.DomainIn()] || blp.domains[fd.DomainOut()]))
+                    continue;
                 }
+              for(auto& pnum : sel.PNums())
+                if(mapto[pnum].IsValid())
+                  pnum = mapto[pnum];
             }
 
           for(ElementIndex ei : Range(ne))
             {
               auto& el = mesh[ei];
-              bool to_move = blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()];
-              if(blp.domains.Size() == 0 || to_move)
+              // only move the elements on the correct side
+              if(blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()])
                 for(auto& pnum : el.PNums())
-                  // Check (Doublecheck) if the corresponding point has a
-                  // copy available for remapping
                   if(mapto[pnum].IsValid())
-                    // Map the volume elements to the new points
                     pnum = mapto[pnum];
             }
 
diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp
index 6523c056..7ea1086b 100644
--- a/libsrc/meshing/curvedelems.cpp
+++ b/libsrc/meshing/curvedelems.cpp
@@ -748,21 +748,19 @@ namespace netgen
 
 	      for (int i2 = 0; i2 < edgenrs.Size(); i2++)
 		{
-		  // PointIndex pi1 = el[edges[i2][0]];
-		  // PointIndex pi2 = el[edges[i2][1]];
-
-		  // bool swap = pi1 > pi2;
-		
-		  // Point<3> p1 = mesh[pi1];
-		  // Point<3> p2 = mesh[pi2];
-		
-		  // int order1 = edgeorder[edgenrs[i2]];
-		  // int ndof = max (0, order1-1);
-
-		  surfnr[edgenrs[i2]] = mesh.GetFaceDescriptor(el.GetIndex()).SurfNr();
-		  gi0[edgenrs[i2]] = el.GeomInfoPi(edges[i2][0]+1);
-		  gi1[edgenrs[i2]] = el.GeomInfoPi(edges[i2][1]+1);
-		}
+		  auto enr = edgenrs[i2];
+                  surfnr[enr] = mesh.GetFaceDescriptor(el.GetIndex()).SurfNr();
+                  if (el[edges[i2][0]] < el[edges[i2][1]])
+                    {
+                      gi0[enr] = el.GeomInfoPi(edges[i2][0]+1);
+                      gi1[enr] = el.GeomInfoPi(edges[i2][1]+1);
+                    }
+                  else
+                    {
+                      gi1[enr] = el.GeomInfoPi(edges[i2][0]+1);
+                      gi0[enr] = el.GeomInfoPi(edges[i2][1]+1);
+                    }
+                }
 	    }
 
 
@@ -1303,6 +1301,8 @@ namespace netgen
                       SurfaceElementIndex sei = top.GetFace2SurfaceElement (f+1)-1;
 		      if (sei != SurfaceElementIndex(-1)) {
 			PointGeomInfo gi = mesh[sei].GeomInfoPi(1);
+                        gi.u = 1.0/3.0*(mesh[sei].GeomInfoPi(1).u+mesh[sei].GeomInfoPi(2).u+mesh[sei].GeomInfoPi(3).u);
+                        gi.v = 1.0/3.0*(mesh[sei].GeomInfoPi(1).v+mesh[sei].GeomInfoPi(2).v+mesh[sei].GeomInfoPi(3).v);
 			geo.ProjectPointGI(surfnr[facenr], pp, gi);
 		      }
 		      else
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index b38aeb88..76a22ffc 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -6541,7 +6541,7 @@ namespace netgen
       return defaultstring;
 
     if (bcnr < 0 || bcnr >= bcnames.Size())
-      throw NgException ("illegal bc-number");
+      throw RangeException("Illegal bc number ", bcnr, 0, bcnames.Size());
 
     if ( bcnames[bcnr] )
       return *bcnames[bcnr];
diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp
index 88955a36..84f5d399 100644
--- a/libsrc/meshing/meshfunc2d.cpp
+++ b/libsrc/meshing/meshfunc2d.cpp
@@ -20,6 +20,19 @@ namespace netgen
       }
     mesh.Compress();
 
+    bool optimize_swap_separate_faces = false;
+    if(!mp.quad)
+    {
+      bool mixed = false;
+      ParallelFor( Range(mesh.GetNSE()), [&] (auto i) NETGEN_LAMBDA_INLINE
+          {
+            if (mesh[SurfaceElementIndex(i)].GetNP() != 3)
+                mixed = true;
+          });
+      if(mixed)
+        optimize_swap_separate_faces = true;
+    }
+
     const char * optstr = mp.optimize2d.c_str();
     int optsteps = mp.optsteps2d;
 
@@ -31,16 +44,42 @@ namespace netgen
 	    {
 	    case 's': 
 	      {  // topological swap
-		MeshOptimize2d meshopt(mesh);
+                MeshOptimize2d meshopt(mesh);
                 meshopt.SetMetricWeight (mp.elsizeweight);
-		meshopt.EdgeSwapping (0);
+
+                if(optimize_swap_separate_faces)
+                {
+                  for(auto i : Range(1, mesh.GetNFD()+1))
+                  {
+                    meshopt.SetFaceIndex(i);
+                    meshopt.EdgeSwapping (0);
+                  }
+                }
+                else
+                {
+                  meshopt.SetFaceIndex(0);
+                  meshopt.EdgeSwapping (0);
+                }
 		break;
 	      }
 	    case 'S': 
 	      {  // metric swap
 		MeshOptimize2d meshopt(mesh);
                 meshopt.SetMetricWeight (mp.elsizeweight);
-		meshopt.EdgeSwapping (1);
+
+                if(optimize_swap_separate_faces)
+                {
+                  for(auto i : Range(1, mesh.GetNFD()+1))
+                  {
+                    meshopt.SetFaceIndex(i);
+                    meshopt.EdgeSwapping (1);
+                  }
+                }
+                else
+                {
+                  meshopt.SetFaceIndex(0);
+                  meshopt.EdgeSwapping (1);
+                }
 		break;
 	      }
 	    case 'm': 
diff --git a/libsrc/meshing/meshing.hpp b/libsrc/meshing/meshing.hpp
index 35cb9add..253f2f93 100644
--- a/libsrc/meshing/meshing.hpp
+++ b/libsrc/meshing/meshing.hpp
@@ -57,10 +57,12 @@ namespace netgen
 #include "hprefinement.hpp"
 #include "boundarylayer.hpp"
 #include "specials.hpp"
+
 }
 
 #include "validate.hpp"
 #include "basegeom.hpp"
+#include "surfacegeom.hpp"
 
 #ifdef PARALLEL
 #include "paralleltop.hpp"
diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp
index 99a5f981..fac0b2ef 100644
--- a/libsrc/meshing/parallelmesh.cpp
+++ b/libsrc/meshing/parallelmesh.cpp
@@ -798,7 +798,7 @@ namespace netgen
 
 
 
-  // slaves receive the mesh from the master
+  // workers receive the mesh from the master
   void Mesh :: ReceiveParallelMesh ( )
   {
     int timer = NgProfiler::CreateTimer ("ReceiveParallelMesh");
@@ -1058,7 +1058,7 @@ namespace netgen
 
   
   
-  // distribute the mesh to the slave processors
+  // distribute the mesh to the worker processors
   // call it only for the master !
   void Mesh :: Distribute ()
   {
@@ -1358,7 +1358,7 @@ namespace netgen
 
 
 
-  // distribute the mesh to the slave processors
+  // distribute the mesh to the worker processors
   // call it only for the master !
   void Mesh :: Distribute (NgArray<int> & volume_weights , NgArray<int>  & surface_weights, NgArray<int>  & segment_weights)
   {
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index fdfb2c40..d2d57f89 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -949,8 +949,22 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
                {
                  regex pattern(*get_if<string>(&boundary));
                  for(int i = 1; i<=self.GetNFD(); i++)
-                   if(regex_match(self.GetFaceDescriptor(i).GetBCName(), pattern))
-                     blp.surfid.Append(i);
+                   {
+                     auto& fd = self.GetFaceDescriptor(i);
+                     if(regex_match(fd.GetBCName(), pattern))
+                       {
+                         auto dom_pattern = get_if<string>(&domain);
+                         // only add if adjacent to domain
+                         if(dom_pattern)
+                           {
+                             regex pattern(*dom_pattern);
+                             if(regex_match(self.GetMaterial(fd.DomainIn()), pattern) || (fd.DomainOut() > 0 ? regex_match(self.GetMaterial(fd.DomainOut()), pattern) : false))
+                               blp.surfid.Append(i);
+                           }
+                         else
+                           blp.surfid.Append(i);
+                       }
+                     }
                }
 
              if(double* pthickness = get_if<double>(&thickness); pthickness)
@@ -1132,6 +1146,45 @@ grow_edges : bool = False
 
 
   m.def("ReadCGNSFile", &ReadCGNSFile, py::arg("filename"), py::arg("base")=1, "Read mesh and solution vectors from CGNS file");
+
+    py::class_<SurfaceGeometry, NetgenGeometry, shared_ptr<SurfaceGeometry>> (m, "SurfaceGeometry")
+    .def(py::init<>())
+    .def(py::init([](py::object pyfunc)
+                  {
+                    std::function<Vec<3> (Point<2>)> func = [pyfunc](Point<2> p)
+                                                    {
+                                                      py::gil_scoped_acquire aq;
+                                                      py::tuple pyres = py::extract<py::tuple>(pyfunc(p[0],p[1],0.0)) ();
+                                                      return Vec<3>(py::extract<double>(pyres[0])(),py::extract<double>(pyres[1])(),py::extract<double>(pyres[2])());
+                                                    };
+                    auto geo = make_shared<SurfaceGeometry>(func);
+                    return geo;
+                  }), py::arg("mapping"))
+    .def(NGSPickle<SurfaceGeometry>())
+    .def("GenerateMesh", [](shared_ptr<SurfaceGeometry> geo,
+                            bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames)
+           {
+             if (py::len(py_bbbpts) != py::len(py_bbbnames))
+               throw Exception("In SurfaceGeometry::GenerateMesh bbbpts and bbbnames do not have same lengths.");
+             Array<Point<3>> bbbpts(py::len(py_bbbpts));
+             Array<string> bbbname(py::len(py_bbbpts));
+             for(int i = 0; i<py::len(py_bbbpts);i++)
+		 {
+                   py::tuple pnt = py::extract<py::tuple>(py_bbbpts[i])();
+                   bbbpts[i] = Point<3>(py::extract<double>(pnt[0])(),py::extract<double>(pnt[1])(),py::extract<double>(pnt[2])());
+                   bbbname[i] = py::extract<string>(py_bbbnames[i])();
+                 }
+             auto mesh = make_shared<Mesh>();
+             SetGlobalMesh (mesh);
+             mesh->SetGeometry(geo);
+	     ng_geometry = geo;
+             auto result = geo->GenerateMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname);
+             if(result != 0)
+               throw Exception("SurfaceGeometry: Meshing failed!");
+             return mesh;
+           }, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list())
+      ;
+    ;
 }
 
 PYBIND11_MODULE(libmesh, m) {
diff --git a/libsrc/meshing/surfacegeom.cpp b/libsrc/meshing/surfacegeom.cpp
new file mode 100644
index 00000000..40d272a8
--- /dev/null
+++ b/libsrc/meshing/surfacegeom.cpp
@@ -0,0 +1,419 @@
+/* *************************************************************************/
+/* File:   surfacegeom.cpp                                                 */
+/* Author: Michael Neunteufel                                              */
+/* Date:   Jun. 2020                                                       */
+/* *************************************************************************/
+
+#include <meshing.hpp>
+
+namespace netgen
+{
+  SurfaceGeometry :: SurfaceGeometry()
+  {
+    //identity
+    func = [](Point<2> p) { return Vec<3>(p[0],p[1],0.0); };
+  }
+
+  SurfaceGeometry :: SurfaceGeometry(function<Vec<3>(Point<2>)> _func) : func(_func)
+  {
+    ;
+  }
+
+  SurfaceGeometry :: SurfaceGeometry(const SurfaceGeometry& geom) : func(geom.func), eps(geom.eps)
+  {
+    ;
+  }
+  
+  void SurfaceGeometry :: CalcHesse(double u, double v, Vec<3>& f_uu, Vec<3>& f_vv, Vec<3>& f_uv) const
+  {
+    Point<2> p = Point<2>(u,v);
+    double pr  = p[0]+eps;
+    double pl  = p[0]-eps;
+    double prr = p[0]+2*eps;
+    double pll = p[0]-2*eps;
+    
+    auto dr = GetTangentVectors( pr, v );
+    auto dl = GetTangentVectors( pl, v );
+    auto drr = GetTangentVectors( prr, v );
+    auto dll = GetTangentVectors( pll, v );
+    
+    f_uu = (1.0/(12.0*eps)) * (8.0*dr[0]-8.0*dl[0]-drr[0]+dll[0]);
+    f_uv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]);
+
+    pr  = p[1]+eps;
+    pl  = p[1]-eps;
+    prr = p[1]+2*eps;
+    pll = p[1]-2*eps;
+
+    dr = GetTangentVectors(u, pr);
+    dl = GetTangentVectors(u, pl);
+    drr = GetTangentVectors(u, prr);
+    dll = GetTangentVectors(u, pll);
+    
+    f_vv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]);
+  }
+  
+  Array<Vec<3>> SurfaceGeometry :: GetTangentVectors(double u, double v) const
+  {
+    Array<Vec<3>> tang(2);
+    
+    Point<2> pru  = Point<2>(u+eps,v);
+    Point<2> plu  = Point<2>(u-eps,v);
+    Point<2> prru = Point<2>(u+2*eps,v);
+    Point<2> pllu = Point<2>(u-2*eps,v);
+    
+    Point<2> prv  = Point<2>(u,v+eps);
+    Point<2> plv  = Point<2>(u,v-eps);
+    Point<2> prrv = Point<2>(u,v+2*eps);
+    Point<2> pllv = Point<2>(u,v-2*eps);
+
+    
+    tang[0] = 1/(12.0*eps)*( 8.0*func(pru) - 8.0*func(plu) - func(prru) + func(pllu) );
+    tang[1] = 1/(12.0*eps)*( 8.0*func(prv) - 8.0*func(plv) - func(prrv) + func(pllv) );
+    
+    return tang;
+  }
+
+  Vec<3> SurfaceGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const
+  {
+    Array<Vec<3>> tang = GetTangentVectors(gi->u, gi->v);
+    auto normal = Cross(tang[0], tang[1]);
+    return Cross(tang[0], tang[1]);
+  }
+
+
+  PointGeomInfo SurfaceGeometry :: ProjectPoint(int surfind, Point<3> & p) const
+  {
+    throw Exception("In SurfaceGeometry::ProjectPoint");
+  }
+  
+  void SurfaceGeometry :: ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
+                                            EdgePointGeomInfo* gi) const
+  {
+    if (gi == nullptr)
+      throw Exception("In SurfaceGeometry::ProjectPointEdge: gi is nullptr");
+    throw Exception("In SurfaceGeometry::ProjectPointEdge: not implemented");
+  }
+    
+  bool  SurfaceGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const
+  {
+    Array<Vec<3>> tangs;
+    Vec<3> diff, f_uu, f_vv, f_uv;
+    Vec<2> r, dx;
+    double norm_r, det, energy=0.0, new_energy=0.0, alpha=2.0,u=0.0,v=0.0;
+    Mat<2,2> mat, inv;
+    int num=0, maxit=20;
+    double damping=0.2;
+
+
+    //Solve minimization problem
+    //   argmin_(u,v) 0.5*\| f(u,v)-p\|^2
+    //via Neton's method:
+    //  F(u,v) = ( (f(u,v)-p)*f_u(u,v), (f(u,v)-p)*f_v(u,v))^T = (0,0)^T
+    //Stiffness matrix
+    //  F'(u,v) = ( f_u*f_u + (f-p)*f_uu, f_v*f_u + (f-p)*f_uv, f_v*f_u + (f-p)*f_uv, f_v*f_v + (f-p)*f_vv )
+    do
+      {
+        num++;
+        tangs = GetTangentVectors(gi.u, gi.v);
+        diff = func(Point<2>(gi.u, gi.v)) - Vec<3>(p);
+        energy = diff.Length2();
+        r = Vec<2>( diff*tangs[0], diff*tangs[1] );
+        norm_r = r.Length2();
+
+        CalcHesse(gi.u, gi.v, f_uu, f_vv, f_uv);
+        
+       
+        mat(0,0) = tangs[0]*tangs[0] + diff*f_uu;
+        mat(1,0) = mat(0,1) = tangs[0]*tangs[1]+diff*f_uv;
+        mat(1,1) = tangs[1]*tangs[1]+diff*f_vv;
+    
+        CalcInverse(mat,inv);
+
+        dx = inv*r;
+
+        //Linesearch 
+        alpha = 2.0;
+        do
+          {
+            alpha /= 2.0;
+            u = gi.u - min(1.0,alpha*damping*num)*dx[0];
+            v = gi.v - min(1.0,alpha*damping*num)*dx[1];
+
+            diff = func(Point<2>(u, v)) - Vec<3>(p);
+            new_energy = diff.Length2();
+          }
+        while (alpha > 1e-10 && new_energy > energy+1e-14);
+        if (alpha <= 1e-10)
+          throw Exception("In SurfaceGeometry::ProjectPointGI: Linesearch min alpha reached!");
+        gi.u = u;
+        gi.v = v;
+        
+
+      }
+    while ( norm_r > 1e-12 && num < maxit);
+
+    //Stay in reference domain [0,1]^2
+    if (gi.u < 0 || gi.u > 1 || gi.v < 0 || gi.v > 1)
+      {
+        cout << "Warning: Projected point outside [0,1]^2: u=" << gi.u << ",v=" << gi.v <<". Setting back." << endl;
+        gi.u = min(max(gi.u,0.0),1.0);
+        gi.v = min(max(gi.v,0.0),1.0);
+      }
+
+    p = Point<3>(func(Point<2>(gi.u,gi.v)));
+
+    if (num == maxit)
+      {
+        //cout << "In SurfaceGeometry::ProjectPointGI: Newton did not converge" << endl;
+        throw Exception("In SurfaceGeometry::ProjectPointGI: Newton did not converge");
+      }
+    return true;
+  }
+
+  bool  SurfaceGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const
+  {
+    throw Exception("In SurfaceGeometry::CalcPointGeomInfo: not implemented");
+    return false;
+  }
+
+  void  SurfaceGeometry :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint, int surfi1, int surfi2, const EdgePointGeomInfo & ap1, const EdgePointGeomInfo & ap2, Point<3> & newp, EdgePointGeomInfo & newgi) const
+  {
+    newp = p1+secpoint*(p2-p1);
+
+    PointGeomInfo pgi;
+    pgi.u = ap1.u+secpoint*(ap2.u-ap1.u);
+    pgi.v = ap1.v+secpoint*(ap2.v-ap1.v);
+
+    ProjectPointGI(surfi1, newp, pgi);
+
+    newgi.u = pgi.u;
+    newgi.v = pgi.v;
+    newgi.edgenr = ap1.edgenr;
+    newgi.body = -1;
+    newgi.dist = -1.0;
+  }
+  
+  void SurfaceGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
+                                       int surfi, 
+                                       const PointGeomInfo & gi1, 
+                                       const PointGeomInfo & gi2,
+                                       Point<3> & newp, PointGeomInfo & newgi) const
+  {
+    newp = p1+secpoint*(p2-p1);
+    
+    newgi.u = gi1.u+secpoint*(gi2.u-gi1.u);
+    newgi.v = gi1.v+secpoint*(gi2.v-gi1.v);  
+    newgi.trignum = -1;
+
+    ProjectPointGI(surfi, newp, newgi);
+  }
+
+  int SurfaceGeometry :: GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames)
+  {
+    mesh->SetDimension(3);
+
+    Array<bool> found(bbbpts.Size());
+    found = false;
+    Array<PointIndex> indbbbpts(bbbpts.Size());
+
+    
+    Array<PointIndex> pids;
+    Array<PointGeomInfo> pgis;
+    for(int i=0; i <= ny; i++)
+      for(int j=0; j <= nx; j++)
+        {
+          PointGeomInfo pgi;
+          pgi.trignum = -1;
+          pgi.u = double(j)/nx;
+          pgi.v = double(i)/ny;
+
+          Point<3> pnt = Point<3>(func(Point<2>(pgi.u,pgi.v)));
+          pids.Append(mesh->AddPoint(pnt));
+          pgis.Append(pgi);
+          
+          for (int k = 0; k < bbbpts.Size(); k++)
+            {
+              auto diff = pnt - bbbpts[k];
+              if(diff.Length2() < 1e-14)
+                {
+                  found[k] = true;
+                  indbbbpts[k] = pids[pids.Size()-1];
+                }
+            }
+        }
+
+    for (bool f : found)
+      if (!f)
+        throw Exception("In SurfaceGeometry :: GenerateMesh: bbbpts not resolved in mesh.");
+
+    FaceDescriptor fd;
+    fd.SetSurfNr(1);
+    fd.SetDomainIn(1);
+    fd.SetDomainOut(0);
+    fd.SetBCProperty(1);
+    mesh->AddFaceDescriptor(fd);
+
+
+    for(int i=0; i < ny; i++)
+      {
+        for(int j=0; j < nx; j++)
+          {
+            int base = i * (nx+1) + j;
+            if (quads)
+              {
+                int pnum[4] = {base,base+1,base+nx+2,base+nx+1};
+                Element2d el = Element2d(QUAD);
+                for (int i = 0; i < 4; i++)
+                  {
+                    el[i] = pids[pnum[i]];
+                    el.GeomInfoPi(i+1) = pgis[pnum[i]];
+                  }
+                el.SetIndex(1);
+            
+                mesh->AddSurfaceElement(el);
+              }
+            else
+              {
+                Array<int> pnum1(3);
+                Array<int> pnum2(3);
+                if (flip_triangles)
+                  {
+                    pnum1[0] = base;
+                    pnum1[1] = base+1;
+                    pnum1[2] = base+nx+2;
+                    pnum2[0] = base;
+                    pnum2[1] = base+nx+2;
+                    pnum2[2] = base+nx+1;
+                  }
+                else
+                  {
+                    pnum1[0] = base;
+                    pnum1[1] = base+1;
+                    pnum1[2] = base+nx+1;
+                    pnum2[0] = base+1;
+                    pnum2[1] = base+nx+2;
+                    pnum2[2] = base+nx+1;
+                  }
+
+                Element2d el = Element2d(TRIG);
+                for (int i = 0; i < 3; i++)
+                  {
+                    el[i] = pids[pnum1[i]];
+                    el.GeomInfoPi(i+1) = pgis[pnum1[i]];
+                  }
+                el.SetIndex(1);
+            
+                mesh->AddSurfaceElement(el);
+                for (int i = 0; i < 3; i++)
+                  {
+                    el[i] = pids[pnum2[i]];
+                    el.GeomInfoPi(i+1) = pgis[pnum2[i]];
+                  }
+                mesh->AddSurfaceElement(el);
+              }
+          }
+      }
+
+    Segment seg;
+    seg.si = 1;
+    seg.edgenr = 1;
+    seg.epgeominfo[0].edgenr = 1;
+    seg.epgeominfo[1].edgenr = 1;
+    // needed for codim2 in 3d
+    seg.edgenr = 1;
+    for(int i=0; i < nx; i++)
+      {
+        seg[0] = pids[i];
+        seg[1] = pids[i+1];
+        
+        seg.geominfo[0] = pgis[i];
+        seg.geominfo[1] = pgis[i+1];
+        seg.epgeominfo[0].u = pgis[i].u;
+        seg.epgeominfo[0].v = pgis[i].v;
+        seg.epgeominfo[0].edgenr = seg.edgenr;
+        seg.epgeominfo[1].u = pgis[i+1].u;
+        seg.epgeominfo[1].v = pgis[i+1].v;
+        seg.epgeominfo[1].edgenr = seg.edgenr;
+        
+        mesh->AddSegment(seg);
+      }
+
+    seg.si = 2;
+    seg.edgenr = 2;
+    for(int i=0; i<ny; i++)
+      {
+        seg[0] = pids[i*(nx+1)+nx];
+        seg[1] = pids[(i+1)*(nx+1)+nx];
+
+        seg.geominfo[0] = pgis[i*(nx+1)+nx];
+        seg.geominfo[1] = pgis[(i+1)*(nx+1)+nx];
+        seg.epgeominfo[0].u = pgis[i*(nx+1)+nx].u;
+        seg.epgeominfo[0].v = pgis[i*(nx+1)+nx].v;
+        seg.epgeominfo[0].edgenr = seg.edgenr;
+        seg.epgeominfo[1].u = pgis[(i+1)*(nx+1)+nx].u;
+        seg.epgeominfo[1].v = pgis[(i+1)*(nx+1)+nx].v;
+        seg.epgeominfo[1].edgenr = seg.edgenr;
+
+        mesh->AddSegment(seg);
+      }
+
+    seg.si = 3;
+    seg.edgenr = 3;
+    for(int i=0; i<nx; i++)
+      {
+        seg[0] = pids[ny*(nx+1)+i+1];
+        seg[1] = pids[ny*(nx+1)+i];
+
+        seg.geominfo[0] = pgis[ny*(nx+1)+i+1];
+        seg.geominfo[1] = pgis[ny*(nx+1)+i];
+        seg.epgeominfo[0].u = pgis[ny*(nx+1)+i+1].u;
+        seg.epgeominfo[0].v = pgis[ny*(nx+1)+i+1].v;
+        seg.epgeominfo[0].edgenr = seg.edgenr;
+        seg.epgeominfo[1].u = pgis[ny*(nx+1)+i].u;
+        seg.epgeominfo[1].v = pgis[ny*(nx+1)+i].v;
+        seg.epgeominfo[1].edgenr = seg.edgenr;
+        
+        mesh->AddSegment(seg);
+      }
+
+    seg.si = 4;
+    seg.edgenr = 4;
+    for(int i=0; i<ny; i++)
+      {
+        seg[0] = pids[(i+1)*(nx+1)];
+        seg[1] = pids[i*(nx+1)];
+
+        seg.geominfo[0] = pgis[(i+1)*(nx+1)];
+        seg.geominfo[1] = pgis[i*(nx+1)];
+        seg.epgeominfo[0].u = pgis[(i+1)*(nx+1)].u;
+        seg.epgeominfo[0].v = pgis[(i+1)*(nx+1)].v;
+        seg.epgeominfo[0].edgenr = seg.edgenr;
+        seg.epgeominfo[1].u = pgis[i*(nx+1)].u;
+        seg.epgeominfo[1].v = pgis[i*(nx+1)].v;
+        seg.epgeominfo[1].edgenr = seg.edgenr;
+
+        mesh->AddSegment(seg);
+      }
+
+    mesh->SetCD2Name(1, "bottom");
+    mesh->SetCD2Name(2, "right");
+    mesh->SetCD2Name(3, "top");
+    mesh->SetCD2Name(4, "left");
+
+    for (int i = 0; i < bbbpts.Size(); i++)
+      {
+        Element0d el;
+        el.pnum = indbbbpts[i];
+        el.index = i+1;
+        mesh->pointelements.Append(el);
+        mesh->SetCD3Name(i+1, bbbnames[i]);
+      }
+
+    mesh->Compress();
+    mesh->UpdateTopology();
+
+    return 0;
+  }
+  
+};
diff --git a/libsrc/meshing/surfacegeom.hpp b/libsrc/meshing/surfacegeom.hpp
new file mode 100644
index 00000000..25ca73c1
--- /dev/null
+++ b/libsrc/meshing/surfacegeom.hpp
@@ -0,0 +1,70 @@
+#ifndef FILE_SURFACEGEOM
+#define FILE_SURFACEGEOM
+
+/* *************************************************************************/
+/* File:   surfacegeom.hpp                                                 */
+/* Author: Michael Neunteufel                                              */
+/* Date:   Jun. 2020                                                       */
+/* *************************************************************************/
+
+
+#include <functional>
+
+
+namespace netgen
+{
+
+  class DLL_HEADER SurfaceGeometry : public NetgenGeometry
+  {
+    function<Vec<3>(Point<2>)> func;
+    double eps=1e-4;
+
+  private:
+
+    void CalcHesse(double u, double v, Vec<3>& f_uu, Vec<3>& f_vv, Vec<3>& f_uv) const;
+  public:
+
+    SurfaceGeometry();
+    SurfaceGeometry(function<Vec<3>(Point<2>)> func);
+    SurfaceGeometry(const SurfaceGeometry& geom);
+    SurfaceGeometry& operator =(const SurfaceGeometry& geom)
+    {
+      func = geom.func;
+      eps = geom.eps;
+      return *this;
+    }
+
+    Array<Vec<3>> GetTangentVectors(double u, double v) const;
+
+    virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override;
+
+    virtual PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override;
+    
+    virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p,
+                           EdgePointGeomInfo* gi = nullptr) const override;
+    
+    virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override;
+
+    virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override;
+
+    virtual void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint,
+                          int surfi1, int surfi2, 
+                          const EdgePointGeomInfo & ap1, 
+                          const EdgePointGeomInfo & ap2,
+                          Point<3> & newp, EdgePointGeomInfo & newgi) const override;
+    
+    virtual void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint,
+                      int surfi, 
+                      const PointGeomInfo & gi1, 
+                      const PointGeomInfo & gi2,
+                      Point<3> & newp, PointGeomInfo & newgi) const override;
+
+    int GenerateMesh(shared_ptr<Mesh> & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array<Point<3>>& bbbpts, const Array<string>& bbbnames);
+
+  };
+  
+}
+   
+
+
+#endif //SURFACEGEOM
diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp
index 037022f9..91dfb2cf 100644
--- a/libsrc/occ/occgeom.cpp
+++ b/libsrc/occ/occgeom.cpp
@@ -169,8 +169,91 @@ namespace netgen
 
    }
 
+  void OCCGeometry :: GlueGeometry()
+  {
+    PrintMessage(1, "OCC Glue Geometry");
+    /*
+      // 
+    BRep_Builder builder;
+    TopoDS_Shape my_fuse;
+    int cnt = 0;
+    for (TopExp_Explorer exp_solid(shape, TopAbs_SOLID); exp_solid.More(); exp_solid.Next())
+      {
+        cout << "cnt = " << cnt << endl;
+	if (cnt == 0)
+	  my_fuse = exp_solid.Current();
+	else
+          // my_fuse = BRepAlgoAPI_Fuse (my_fuse, exp_solid.Current());
+          my_fuse = QANewModTopOpe_Glue::QANewModTopOpe_Glue(my_fuse, exp_solid.Current());
+	cnt++;
+      }
+    cout << "remove" << endl;
+    // for (int i = 1; i <= somap.Size(); i++)
+    // builder.Remove (shape, somap(i));
+    cout << "now add" << endl;
+    // builder.Add (shape, my_fuse);
+    shape = my_fuse;
+    cout << "build fmap" << endl;
+    BuildFMap();
+    */
 
 
+    // from 
+    // https://www.opencascade.com/doc/occt-7.4.0/overview/html/occt_user_guides__boolean_operations.html
+    BOPAlgo_Builder aBuilder;
+
+    // Setting arguments
+    TopTools_ListOfShape aLSObjects; 
+    for (TopExp_Explorer exp_solid(shape, TopAbs_SOLID); exp_solid.More(); exp_solid.Next())
+      aLSObjects.Append (exp_solid.Current());
+    aBuilder.SetArguments(aLSObjects);
+
+    // Setting options for GF
+    // Set parallel processing mode (default is false)
+    // Standard_Boolean bRunParallel = Standard_True;
+    // aBuilder.SetRunParallel(bRunParallel);
+    
+    // Set Fuzzy value (default is Precision::Confusion())
+    // Standard_Real aFuzzyValue = 1.e-5;
+    // aBuilder.SetFuzzyValue(aFuzzyValue);
+    
+    // Set safe processing mode (default is false)
+    // Standard_Boolean bSafeMode = Standard_True;
+    // aBuilder.SetNonDestructive(bSafeMode);
+    
+    // Set Gluing mode for coinciding arguments (default is off)
+    // BOPAlgo_GlueEnum aGlue = BOPAlgo_GlueShift;
+    // aBuilder.SetGlue(aGlue);
+    
+    // Disabling/Enabling the check for inverted solids (default is true)
+    // Standard Boolean bCheckInverted = Standard_False;
+    // aBuilder.SetCheckInverted(bCheckInverted);
+    
+    // Set OBB usage (default is false)
+    // Standard_Boolean bUseOBB = Standard_True;
+    // aBuilder.SetUseOBB(buseobb);
+    
+    // Perform the operation
+    aBuilder.Perform();
+    // Check for the errors
+#if OCC_VERSION_HEX >= 0x070200
+    if (aBuilder.HasErrors())
+      {
+        cout << "builder has errors" << endl;
+        return;
+      }
+    // Check for the warnings
+    if (aBuilder.HasWarnings())
+      {
+        // treatment of the warnings
+        ;
+      }
+#endif
+    // result of the operation
+    shape = aBuilder.Shape();
+    BuildFMap();
+  }
+
    void OCCGeometry :: HealGeometry ()
    {
       int nrc = 0, nrcs = 0,
diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp
index 372979ad..447d212f 100644
--- a/libsrc/occ/occgeom.hpp
+++ b/libsrc/occ/occgeom.hpp
@@ -78,7 +78,7 @@
 #include "Bnd_Box.hxx"
 #include "ShapeAnalysis.hxx"
 #include "ShapeBuild_ReShape.hxx"
-
+#include "BOPAlgo_Builder.hxx"
 
 // Philippose - 29/01/2009
 // OpenCascade XDE Support
@@ -343,6 +343,7 @@ namespace netgen
     void MakeSolid();
 
     void HealGeometry();
+    void GlueGeometry();
 
     // Philippose - 15/01/2009
     // Sets the maximum mesh size for a given face
diff --git a/libsrc/occ/python_occ.cpp b/libsrc/occ/python_occ.cpp
index 837ecdb0..a7508c31 100644
--- a/libsrc/occ/python_occ.cpp
+++ b/libsrc/occ/python_occ.cpp
@@ -66,6 +66,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                   }), py::arg("filename"),
         "Load OCC geometry from step, brep or iges file")
     .def(NGSPickle<OCCGeometry>())
+    .def("Glue", &OCCGeometry::GlueGeometry)
     .def("Heal",[](OCCGeometry & self, double tolerance, bool fixsmalledges, bool fixspotstripfaces, bool sewfaces, bool makesolids, bool splitpartitions)
          {
            self.tolerance = tolerance;
diff --git a/tests/pytest/test_splinegeo_tensordomainmeshing.py b/tests/pytest/test_splinegeo_tensordomainmeshing.py
new file mode 100644
index 00000000..0c00ae7d
--- /dev/null
+++ b/tests/pytest/test_splinegeo_tensordomainmeshing.py
@@ -0,0 +1,22 @@
+from netgen.geom2d import *
+
+def test_tensordomainmeshing():
+    geo = SplineGeometry()
+    w = 10
+    h = 0.01
+
+    p = [ (0, 0), (w, 0), (w, h), (0, h) ]
+    p = [geo.AppendPoint(*px) for px in p]
+
+    l0 = geo.Append ( ["line", p[0], p[1]], leftdomain=1, rightdomain=0 )
+    l1 = geo.Append ( ["line", p[1], p[2]], leftdomain=1, rightdomain=0)
+    geo.Append ( ["line", p[3], p[2]], leftdomain=0, rightdomain=1, copy=l0 )
+    geo.Append ( ["line", p[0], p[3]], leftdomain=0, rightdomain=1, copy=l1 )
+
+    geo._SetDomainTensorMeshing(1, True)
+
+    mesh = geo.GenerateMesh(maxh=1)
+
+    for el in mesh.Elements2D():
+        print(el.vertices)
+        assert len(el.vertices) == 4