diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 00000000..fccd5ba0
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "external_dependencies/pybind11"]
+	path = external_dependencies/pybind11
+	url = https://github.com/pybind/pybind11.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e76ca6f2..13e30bb5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -124,7 +124,7 @@ endmacro()
 
 if(WIN32)
     get_WIN32_WINNT(ver)
-    add_definitions(-D_WIN32_WINNT=${ver} -DWNT -DWNT_WINDOW)
+    add_definitions(-D_WIN32_WINNT=${ver} -DWNT -DWNT_WINDOW -DNOMINMAX)
     set(CMAKE_MFC_FLAG 0)
 #     add_definitions(-DNGINTERFACE_EXPORTS)
 #     add_definitions(-DNGLIB_EXPORTS)
@@ -135,7 +135,7 @@ if(WIN32)
 else(WIN32)
     # build shared libraries
     set(NG_LIB_TYPE SHARED)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14")
 endif(WIN32)
 if(APPLE)
 #    set(MACOSX_BUNDLE ON)
@@ -185,6 +185,17 @@ endif (USE_GUI)
 
 #######################################################################
 if (USE_PYTHON)
+    find_path(PYBIND_INCLUDE_DIR pybind11/pybind11.h ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies/pybind11/include)
+    if( NOT PYBIND_INCLUDE_DIR )
+      execute_process(COMMAND git submodule update --init --recursive WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+        find_path(PYBIND_INCLUDE_DIR pybind11/pybind11.h ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies/pybind11/include)
+    endif( NOT PYBIND_INCLUDE_DIR )
+    if( PYBIND_INCLUDE_DIR )
+        message("-- Found Pybind11: ${PYBIND_INCLUDE_DIR}")
+    else( PYBIND_INCLUDE_DIR )
+        message(FATAL_ERROR "Could NOT find pybind11!")
+    endif( PYBIND_INCLUDE_DIR )
+    include_directories(${PYBIND_INCLUDE_DIR})
     set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${PROJECT_SOURCE_DIR}/cmake_modules/python")
     set(PYTHON_VERSION "3" CACHE STRING "")
     set(Python_ADDITIONAL_VERSIONS 3.5)
@@ -198,36 +209,10 @@ if (USE_PYTHON)
         find_package(PythonInterp ${PYTHON_VERSION} EXACT REQUIRED)
         find_package(PythonLibs ${PYTHON_VERSION}  EXACT REQUIRED)
     endif( PYTHON_VERSION STREQUAL "3")
-    set(Boost_FIND_QUIETLY ON)
-
-    foreach( component IN ITEMS python${PYTHON_VERSION_MAJOR} python-py${PYTHON_VERSION_MAJOR}${PYTHON_VERSION_MINOR} python)
-      string(TOUPPER ${component} component_upper)
-      if(NOT Boost_LIBRARIES)
-          find_package(Boost COMPONENTS ${component})
-          set(BOOST_PYTHON_LIB ${Boost_${component_upper}_LIBRARY_RELEASE})
-      endif(NOT Boost_LIBRARIES)
-    endforeach()
-
-    if(NOT Boost_LIBRARIES)
-        message(FATAL_ERROR "Could NOT find Boost-Python" )
-    else(NOT Boost_LIBRARIES)
-        message("-- Found Boost_python: ${BOOST_PYTHON_LIB}")
-    endif(NOT Boost_LIBRARIES)
 
     add_definitions(-DNG_PYTHON)
-    if(WIN32)
-      install(FILES ${BOOST_PYTHON_LIB} DESTINATION lib COMPONENT netgen)
-      install(DIRECTORY ${Boost_INCLUDE_DIR}/boost DESTINATION include COMPONENT netgen)
-    endif(WIN32)
-
-    if(WIN32 AND ${CMAKE_CXX_COMPILER_ID} STREQUAL "Intel")
-        # fix problems with boosts config/auto_link.hpp
-        add_definitions(-DBOOST_LIB_TOOLSET="iw")
-    endif()
-
     include_directories(${PYTHON_INCLUDE_DIRS})
-    include_directories(${Boost_INCLUDE_DIRS})
-    set(PYTHON_LIBS "${PYTHON_LIBRARIES};${Boost_LIBRARIES}")
+    set(PYTHON_LIBS "${PYTHON_LIBRARIES}")
     execute_process(COMMAND ${PYTHON_EXECUTABLE} -c "from distutils.sysconfig import get_python_lib; print(get_python_lib(1,0,''))" OUTPUT_VARIABLE PYTHON_PACKAGES_INSTALL_DIR OUTPUT_STRIP_TRAILING_WHITESPACE)
     message("python install path: ${PYTHON_PACKAGES_INSTALL_DIR}")
 endif (USE_PYTHON)
@@ -288,7 +273,7 @@ if(INSTALL_DEPENDENCIES)
 #    find_library(LIBSTDCXX NAMES stdc++.6)
 #    find_library(LIBGCCS NAMES gcc_s.1)
 #     install( FILES "${Boost_LIBRARIES}" ${LIBGOMP} ${LIBSTDCXX} ${LIBGCCS} DESTINATION ${ng_install_dir_lib} COMPONENT netgen )
-    install( FILES "${Boost_LIBRARIES}" DESTINATION ${ng_install_dir_lib} COMPONENT netgen )
+#     install( FILES "${Boost_LIBRARIES}" DESTINATION ${ng_install_dir_lib} COMPONENT netgen )
     get_filename_component(MY_LIB_DIR ${TK_LIBRARY} DIRECTORY)
     if(APPLE)
         if(TIX_LIBRARY)
@@ -355,8 +340,7 @@ if(UNIX)
         string(SUBSTRING ${temp} 17 -1 UBUNTU_VERSION)
         message("ubuntu version: ${UBUNTU_VERSION}")
 
-        set(BOOST_PACKAGE "libboost-python${Boost_MAJOR_VERSION}.${Boost_MINOR_VERSION}.${Boost_SUBMINOR_VERSION}")
-        set(CPACK_DEBIAN_PACKAGE_DEPENDS "python3, ${BOOST_PACKAGE}, libtk8.5, libtcl8.5, tix, libxmu6")
+        set(CPACK_DEBIAN_PACKAGE_DEPENDS "python3, libtk8.5, libtcl8.5, tix, libxmu6")
         execute_process(COMMAND dpkg --print-architecture OUTPUT_VARIABLE CPACK_DEBIAN_PACKAGE_ARCHITECTURE OUTPUT_STRIP_TRAILING_WHITESPACE)
         set(CPACK_DEBIAN_PACKAGE_MAINTAINER "Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>")
         if(USE_MPI)
diff --git a/external_dependencies/pybind11 b/external_dependencies/pybind11
new file mode 160000
index 00000000..030d10e8
--- /dev/null
+++ b/external_dependencies/pybind11
@@ -0,0 +1 @@
+Subproject commit 030d10e826b87f8cdf0816aa36b9a515fb7d064d
diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp
index 917993a5..fb89cd46 100644
--- a/libsrc/csg/csgeom.cpp
+++ b/libsrc/csg/csgeom.cpp
@@ -1494,7 +1494,7 @@ namespace netgen
     // virtual VisualScene * GetVisualScene (const NetgenGeometry * geom) const;
   };
 
-  extern CSGeometry * ParseCSG (istream & istr);
+  extern CSGeometry * ParseCSG (istream & istr, CSGeometry *instance=nullptr);
 
   NetgenGeometry *  CSGeometryRegister :: Load (string filename) const
   {
diff --git a/libsrc/csg/csgparser.cpp b/libsrc/csg/csgparser.cpp
index c2728bcd..8ac0a6fe 100644
--- a/libsrc/csg/csgparser.cpp
+++ b/libsrc/csg/csgparser.cpp
@@ -822,11 +822,17 @@ namespace netgen
   /*
     Main parsing function for CSG geometry
   */
-  CSGeometry * ParseCSG (istream & istr)
+  CSGeometry * ParseCSG (istream & istr, CSGeometry * instance=nullptr)
   {
     CSGScanner scan(istr);
     
-    geom = new CSGeometry;
+    if (instance)
+    {
+      new (instance) CSGeometry;
+      geom = instance;
+    }
+    else
+      geom = new CSGeometry;
 
     scan.ReadNext();
     if (scan.GetToken() != TOK_RECO)  // keyword 'algebraic3d'
diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp
index d8c65906..e2875bdc 100644
--- a/libsrc/csg/python_csg.cpp
+++ b/libsrc/csg/python_csg.cpp
@@ -11,11 +11,6 @@ namespace netgen
   extern shared_ptr<NetgenGeometry> ng_geometry;
 }
 
-inline void NOOP_Deleter(void *) { ; }
-
-#if defined(_MSC_FULL_VER) && _MSC_FULL_VER == 190024213
-namespace boost { template<> const volatile CSGeometry* get_pointer(const volatile CSGeometry* p) { return p; } }
-#endif
 
 
 // a shadow solid tree using shared pointers.
@@ -164,28 +159,32 @@ inline ostream & operator<< (ostream & ost, const SPSolid & sol)
 
 namespace netgen
 {
-  extern CSGeometry * ParseCSG (istream & istr);
+  extern CSGeometry * ParseCSG (istream & istr, CSGeometry *instance=nullptr);
 }
 
 
 static Transformation<3> global_trafo(Vec<3> (0,0,0));
 
-DLL_HEADER void ExportCSG() 
+DLL_HEADER void ExportCSG(py::module &m)
 {
-  ModuleScope module("csg");
-
-  bp::class_<Point<2>> ("Point2d", bp::init<double,double>()) 
+  py::class_<NGDummyArgument>(m, "NGDummyArgument")
+    .def("__bool__", []( NGDummyArgument &self ) { return false; } )
+    ;
+  
+  py::class_<Point<2>> (m, "Point2d")
+    .def(py::init<double,double>())
     .def ("__str__", &ToString<Point<2>>)
-    .def(bp::self-bp::self)
-    .def(bp::self+Vec<2>())
-    .def(bp::self-Vec<2>())
+    .def(py::self-py::self)
+    .def(py::self+Vec<2>())
+    .def(py::self-Vec<2>())
     ;
 
-  bp::class_<Point<3>> ("Point3d", bp::init<double,double,double>()) 
+  py::class_<Point<3>> (m, "Point3d")
+    .def(py::init<double,double,double>())
     .def ("__str__", &ToString<Point<3>>)
-    .def(bp::self-bp::self)
-    .def(bp::self+Vec<3>())
-    .def(bp::self-Vec<3>())
+    .def(py::self-py::self)
+    .def(py::self+Vec<3>())
+    .def(py::self-Vec<3>())
     ;
 
   bp::def ("Pnt", FunctionPointer
@@ -193,7 +192,13 @@ DLL_HEADER void ExportCSG()
   bp::def ("Pnt", FunctionPointer
            ([](double x, double y) { return Point<2>(x,y); }));
 
-  bp::def ("SetTransformation", FunctionPointer
+           
+  m.def ("Pnt", FunctionPointer
+           ([](double x, double y, double z) { return Point<3>(x,y,z); }));
+  m.def ("Pnt", FunctionPointer
+           ([](double x, double y) { return Point<2>(x,y); }));
+
+  m.def ("SetTransformation", FunctionPointer
            ([](int dir, double angle)
             {
               if (dir > 0)
@@ -201,35 +206,36 @@ DLL_HEADER void ExportCSG()
               else
                 global_trafo = Transformation<3> (Vec<3>(0,0,0));
             }),
-           (bp::arg("dir")=int(0), bp::arg("angle")=0));
+           bp::arg("dir")=int(0), bp::arg("angle")=0);
 
-           
-  
-  bp::class_<Vec<2>> ("Vec2d", bp::init<double,double>()) 
+  py::class_<Vec<2>> (m, "Vec2d")
+    .def(py::init<double,double>())
     .def ("__str__", &ToString<Vec<3>>)
-    .def(bp::self+bp::self)
-    .def(bp::self-bp::self)
-    .def(-bp::self)
-    .def(double()*bp::self)
+    .def(py::self+py::self)
+    .def(py::self-py::self)
+    .def(-py::self)
+    .def(double()*py::self)
     .def("Norm", &Vec<2>::Length)
     ;
 
-  bp::class_<Vec<3>> ("Vec3d", bp::init<double,double,double>()) 
+  py::class_<Vec<3>> (m, "Vec3d")
+    .def(py::init<double,double,double>())
     .def ("__str__", &ToString<Vec<3>>)
-    .def(bp::self+bp::self)
-    .def(bp::self-bp::self)
-    .def(-bp::self)
-    .def(double()*bp::self)
+    .def(py::self+py::self)
+    .def(py::self-py::self)
+    .def(-py::self)
+    .def(double()*py::self)
     .def("Norm", &Vec<3>::Length)
     ;
 
-  bp::def ("Vec", FunctionPointer
-           ([] (double x, double y, double z) { return global_trafo(Vec<3>(x,y,z)); }));
-  bp::def ("Vec", FunctionPointer
+  m.def ("Vec", FunctionPointer
+           ([] (double x, double y, double z) { return global_trafo<3>(x,y,z); }));
+  m.def ("Vec", FunctionPointer
            ([] (double x, double y) { return Vec<2>(x,y); }));
 
 
-  bp::class_<SplineGeometry<2>> ("SplineCurve2d")
+  py::class_<SplineGeometry<2>> (m, "SplineCurve2d")
+    .def(py::init<>())
     .def ("AddPoint", FunctionPointer
           ([] (SplineGeometry<2> & self, double x, double y)
            {
@@ -249,10 +255,7 @@ DLL_HEADER void ExportCSG()
     ;
 
   
-#if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100)
-  bp::register_ptr_to_python<shared_ptr<SPSolid>>();
-#endif
-  bp::class_<SPSolid, shared_ptr<SPSolid>, boost::noncopyable> ("Solid", bp::no_init)
+  py::class_<SPSolid, shared_ptr<SPSolid>> (m, "Solid")
     .def ("__str__", &ToString<SPSolid>)
     .def ("__add__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::UNION, self, other); }))
     .def ("__mul__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::SECTION, self, other); }))
@@ -270,48 +273,48 @@ DLL_HEADER void ExportCSG()
     .def ("mat", FunctionPointer([](shared_ptr<SPSolid> & self, string mat) -> shared_ptr<SPSolid> 
                                  { self->SetMaterial(mat); return self; }))
     .def ("mat", &SPSolid::GetMaterial)
-    .def("col", FunctionPointer([](shared_ptr<SPSolid> & self, bp::list rgb) -> shared_ptr<SPSolid>
+    .def("col", FunctionPointer([](shared_ptr<SPSolid> & self, py::list rgb) -> shared_ptr<SPSolid>
                                 { 
-                                  bp::extract<double> red(rgb[0]);
-                                  bp::extract<double> green(rgb[1]);
-                                  bp::extract<double> blue(rgb[2]);
+                                  py::extract<double> red(rgb[0]);
+                                  py::extract<double> green(rgb[1]);
+                                  py::extract<double> blue(rgb[2]);
                                   self->SetColor(red(),green(),blue());
                                   return self; 
                                 }))
     .def("transp", FunctionPointer([](shared_ptr<SPSolid> & self)->shared_ptr < SPSolid > { self->SetTransparent(); return self; }))
     ;
 
-  bp::def ("Sphere", FunctionPointer([](Point<3> c, double r)
+  m.def ("Sphere", FunctionPointer([](Point<3> c, double r)
                                      {
                                        Sphere * sp = new Sphere (c, r);
                                        Solid * sol = new Solid (sp);
                                        return make_shared<SPSolid> (sol);
                                      }));
-  bp::def ("Plane", FunctionPointer([](Point<3> p, Vec<3> n)
+  m.def ("Plane", FunctionPointer([](Point<3> p, Vec<3> n)
                                     {
                                       Plane * sp = new Plane (p,n);
                                       Solid * sol = new Solid (sp);
                                       return make_shared<SPSolid> (sol);
                                     }));
-  bp::def ("Cone", FunctionPointer([](Point<3> a, Point<3> b, double ra, double rb)
+  m.def ("Cone", FunctionPointer([](Point<3> a, Point<3> b, double ra, double rb)
                                        {
                                          Cone * cyl = new Cone (a, b, ra, rb);
                                          Solid * sol = new Solid (cyl);
                                          return make_shared<SPSolid> (sol);
                                        }));
-  bp::def ("Cylinder", FunctionPointer([](Point<3> a, Point<3> b, double r)
+  m.def ("Cylinder", FunctionPointer([](Point<3> a, Point<3> b, double r)
                                        {
                                          Cylinder * cyl = new Cylinder (a, b, r);
                                          Solid * sol = new Solid (cyl);
                                          return make_shared<SPSolid> (sol);
                                        }));
-  bp::def ("OrthoBrick", FunctionPointer([](Point<3> p1, Point<3> p2)
+  m.def ("OrthoBrick", FunctionPointer([](Point<3> p1, Point<3> p2)
                                          {
                                            OrthoBrick * brick = new OrthoBrick (p1,p2);
                                            Solid * sol = new Solid (brick);
                                            return make_shared<SPSolid> (sol);
                                          }));
-  bp::def ("Revolution", FunctionPointer([](Point<3> p1, Point<3> p2,
+  m.def ("Revolution", FunctionPointer([](Point<3> p1, Point<3> p2,
                                             const SplineGeometry<2> & spline)
                                          {
                                            Revolution * rev = new Revolution (p1, p2, spline);
@@ -319,58 +322,56 @@ DLL_HEADER void ExportCSG()
                                            return make_shared<SPSolid> (sol);
                                          }));
   
-  bp::def ("Or", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
+  m.def ("Or", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
                                  {
                                    return make_shared<SPSolid> (SPSolid::UNION, s1, s2);
                                  }));
-  bp::def ("And", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
+  m.def ("And", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
                                   {
                                     return make_shared<SPSolid> (SPSolid::SECTION, s1, s2);
                                   }));
 
 
-  bp::class_<CSGeometry,shared_ptr<CSGeometry>, boost::noncopyable> ("CSGeometry")
-    .def("__init__", bp::make_constructor (FunctionPointer
-                                           ([](const string & filename)
+  py::class_<CSGeometry,shared_ptr<CSGeometry>> (m, "CSGeometry")
+    .def(py::init<>())
+    .def("__init__", 
+                                           [](CSGeometry *instance, const string & filename)
                                             {
                                               cout << "load geometry";
                                               ifstream ist(filename);
-                                              shared_ptr<CSGeometry> geom(ParseCSG(ist));
-                                              geom -> FindIdenticSurfaces(1e-8 * geom->MaxSize()); 
-                                              return geom;
-                                            })))
-
-    .def("__init__", bp::make_constructor (FunctionPointer
-                                           ([](const bp::list & solidlist)
+                                              ParseCSG(ist, instance);
+                                              instance -> FindIdenticSurfaces(1e-8 * instance->MaxSize()); 
+                                            })
+    .def("__init__",
+                                           [](CSGeometry *instance, const py::list & solidlist)
                                             {
                                               cout << "csg from list";
-                                              auto geom = make_shared<CSGeometry>();
+                                              new (instance) CSGeometry();
                                               for (int i = 0; i < len(solidlist); i++)
                                                 {
-                                                  bp::object obj = solidlist[i];
+                                                  py::object obj = solidlist[i];
                                                   cout << "obj " << i << endl;
 
-                                                  bp::extract<shared_ptr<SPSolid>> solid(solidlist[i]);
+                                                  py::extract<shared_ptr<SPSolid>> solid(solidlist[i]);
                                                   if(solid.check())
                                                     {
                                                       cout << "its a solid" << endl;
-                                                      solid()->AddSurfaces (*geom);
+                                                      solid()->AddSurfaces (*instance);
                                                       solid()->GiveUpOwner();
-                                                      int tlonr = geom->SetTopLevelObject (solid()->GetSolid());
-                                                      geom->GetTopLevelObject(tlonr) -> SetMaterial(solid()->GetMaterial());
+                                                      int tlonr = instance->SetTopLevelObject (solid()->GetSolid());
+                                                      instance->GetTopLevelObject(tlonr) -> SetMaterial(solid()->GetMaterial());
                                                     }
                                                 }
-                                              geom -> FindIdenticSurfaces(1e-8 * geom->MaxSize()); 
-                                              return geom;
-                                            })))
+                                              instance -> FindIdenticSurfaces(1e-8 * instance->MaxSize()); 
+                                            })
 
     .def("Save", FunctionPointer([] (CSGeometry & self, string filename)
                                  {
                                    cout << "save geometry to file " << filename << endl;
                                    self.Save (filename);
                                  }))
-    .def("Add", FunctionPointer
-         ([] (CSGeometry & self, shared_ptr<SPSolid> solid, bp::list bcmod)
+    .def("Add",
+         [] (CSGeometry & self, shared_ptr<SPSolid> solid, py::list bcmod)
           {
             solid->AddSurfaces (self);
             solid->GiveUpOwner();
@@ -380,15 +381,15 @@ DLL_HEADER void ExportCSG()
             self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent());
 
             // bcmod is list of tuples ( solid, bcnr )
-            for (int i = 0; i < bp::len(bcmod); i++)
+            for (int i = 0; i < py::len(bcmod); i++)
               {
-                bp::tuple tup = bp::extract<bp::tuple> (bcmod[i]) ();
-                auto mod_solid = bp::extract<shared_ptr<SPSolid>> (tup[0]) ();
+                py::tuple tup = py::extract<py::tuple> (bcmod[i]) ();
+                auto mod_solid = py::extract<shared_ptr<SPSolid>> (tup[0]) ();
                 int mod_nr = -1;
                 string * bcname = nullptr;
-                bp::object val = tup[1];
-                if (bp::extract<int>(val).check()) mod_nr = bp::extract<int> (val)();
-                if (bp::extract<string>(val).check()) bcname = new string ( bp::extract<string> (val)());
+                py::object val = tup[1];
+                if (py::extract<int>(val).check()) mod_nr = py::extract<int> (val)();
+                if (py::extract<string>(val).check()) bcname = new string ( py::extract<string> (val)());
 
                 Array<int> si;
                 mod_solid -> GetSolid() -> GetSurfaceIndices (si);
@@ -406,9 +407,8 @@ DLL_HEADER void ExportCSG()
                 delete bcname;
               }
             return tlonr;
-
-          }),
-         (bp::arg("self"), bp::arg("solid"), bp::arg("bcmod")=bp::list())
+          },
+          py::arg("solid"), py::arg("bcmod")=py::list()
          )
 
     .def("AddSurface", FunctionPointer
@@ -425,13 +425,13 @@ DLL_HEADER void ExportCSG()
             self.GetTopLevelObject(tlonr) -> SetRGB(solid->GetRed(),solid->GetGreen(),solid->GetBlue());
             self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent());
           }),
-         (bp::arg("self"), bp::arg("surface"), bp::arg("solid"))
+         py::arg("surface"), py::arg("solid")
          )
          
 
     
     .def("CloseSurfaces", FunctionPointer
-         ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, bp::list aslices )
+         ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, py::list aslices )
           {
             Array<int> si1, si2;
             s1->GetSolid()->GetSurfaceIndices (si1);
@@ -443,15 +443,15 @@ DLL_HEADER void ExportCSG()
 
             try
             {
-                int n = bp::len(aslices);
+                int n = py::len(aslices);
                 Array<double> slices(n);
                 for(int i=0; i<n; i++)
                 {
-                    slices[i]= bp::extract<double>(aslices[i])();
+                    slices[i]= py::extract<double>(aslices[i])();
                 }
                 flags.SetFlag("slices", slices);
             }
-            catch( bp::error_already_set const & ) {
+            catch( py::error_already_set const & ) {
                 cout << "caught python error:" << endl;
                 PyErr_Print();
             }
@@ -464,7 +464,7 @@ DLL_HEADER void ExportCSG()
                 domain,
                 flags));
           }),
-         (bp::arg("self"), bp::arg("solid1"), bp::arg("solid2"), bp::arg("slices"))
+         py::arg("solid1"), py::arg("solid2"), py::arg("slices")
          )
     .def("CloseSurfaces", FunctionPointer
          ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, int reflevels)
@@ -484,7 +484,7 @@ DLL_HEADER void ExportCSG()
                 domain,
                 flags));
           }),
-         (bp::arg("self"), bp::arg("solid1"), bp::arg("solid2"), bp::arg("reflevels")=2)
+         py::arg("solid1"), py::arg("solid2"), py::arg("reflevels")=2
          )
     
     .def("PeriodicSurfaces", FunctionPointer
@@ -499,7 +499,7 @@ DLL_HEADER void ExportCSG()
                (self.GetNIdentifications()+1, self, 
                 self.GetSurface (si1[0]), self.GetSurface (si2[0])));
           }),
-         (bp::arg("self"), bp::arg("solid1"), bp::arg("solid2"))
+         py::arg("solid1"), py::arg("solid2")
          )
     
     .def("GetTransparent", FunctionPointer
@@ -507,14 +507,14 @@ DLL_HEADER void ExportCSG()
           {
             return self.GetTopLevelObject(tlonr)->GetTransparent();
           }),
-         (bp::arg("self"), bp::arg("tlonr"))
+         py::arg("tlonr")
          )
     .def("SetTransparent", FunctionPointer
          ([] (CSGeometry & self, int tlonr, bool transparent)
           {
             self.GetTopLevelObject(tlonr)->SetTransparent(transparent);
           }),
-         (bp::arg("self"), bp::arg("tlonr"), bp::arg("transparent"))
+         py::arg("tlonr"), py::arg("transparent")
          )
 
     .def("GetVisible", FunctionPointer
@@ -522,21 +522,21 @@ DLL_HEADER void ExportCSG()
           {
             return self.GetTopLevelObject(tlonr)->GetVisible();
           }),
-         (bp::arg("self"), bp::arg("tlonr"))
+         py::arg("tlonr")
          )
     .def("SetVisible", FunctionPointer
          ([] (CSGeometry & self, int tlonr, bool visible)
           {
             self.GetTopLevelObject(tlonr)->SetVisible(visible);
           }),
-         (bp::arg("self"), bp::arg("tlonr"), bp::arg("visible"))
+         py::arg("tlonr"), py::arg("visible")
          )
     .def("SetBoundingBox", FunctionPointer
          ([] (CSGeometry & self, Point<3> pmin, Point<3> pmax)
           {
             self.SetBoundingBox(Box<3> (pmin, pmax));
           }),
-         (bp::arg("self"), bp::arg("pmin"), bp::arg("pmax"))
+         py::arg("pmin"), py::arg("pmax")
          )
     .def("Draw", FunctionPointer
          ([] (shared_ptr<CSGeometry> self)
@@ -544,13 +544,12 @@ DLL_HEADER void ExportCSG()
              self->FindIdenticSurfaces(1e-6);
              self->CalcTriangleApproximation(0.01, 20);
              ng_geometry = self;
-          }),
-         (bp::arg("self"))
+          })
          )
-    .add_property ("ntlo", &CSGeometry::GetNTopLevelObjects)
+    .def_property_readonly ("ntlo", &CSGeometry::GetNTopLevelObjects)
     ;
 
-  bp::def("GenerateMesh", FunctionPointer
+  m.def("GenerateMesh", FunctionPointer
           ([](shared_ptr<CSGeometry> geo, MeshingParameters & param)
            {
              auto dummy = make_shared<Mesh>();
@@ -570,7 +569,7 @@ DLL_HEADER void ExportCSG()
            }))
     ;
 
-  bp::def("Save", FunctionPointer 
+  m.def("Save", FunctionPointer 
           ([](const Mesh & self, const string & filename, const CSGeometry & geom)
            {
              ostream * outfile;
@@ -588,7 +587,7 @@ DLL_HEADER void ExportCSG()
 
 
 
-  bp::def("ZRefinement", FunctionPointer
+  m.def("ZRefinement", FunctionPointer
           ([](Mesh & mesh, CSGeometry & geom)
           {
             ZRefinementOptions opt;
@@ -598,14 +597,10 @@ DLL_HEADER void ExportCSG()
     ;
 }
 
-
-
-
-
-BOOST_PYTHON_MODULE(libcsg) {
-  ExportCSG();
+PYBIND11_PLUGIN(libcsg) {
+  py::module m("csg", "pybind csg");
+  ExportCSG(m);
+  return m.ptr();
 }
-
-
 #endif
 
diff --git a/libsrc/csg/vscsg.cpp b/libsrc/csg/vscsg.cpp
index 199e9c07..b35b4b10 100644
--- a/libsrc/csg/vscsg.cpp
+++ b/libsrc/csg/vscsg.cpp
@@ -488,29 +488,21 @@ namespace netgen
 
 
 #ifdef NG_PYTHON
-#include <boost/python.hpp>
 #include <../general/ngpython.hpp>
-namespace bp = boost::python;
 
-#if defined(_MSC_FULL_VER) && _MSC_FULL_VER == 190024213
-namespace boost { template<> const volatile netgen::VisualSceneGeometry* get_pointer(const volatile netgen::VisualSceneGeometry* p) { return p; } }
-#endif
-
-DLL_HEADER void ExportCSGVis()
+DLL_HEADER void ExportCSGVis(py::module &m)
 {
 	using namespace netgen;
 
-	ModuleScope module("csgvis");
-
-	bp::class_<VisualSceneGeometry, shared_ptr<VisualSceneGeometry>>
-		("VisualSceneGeometry", bp::no_init)
+	py::class_<VisualSceneGeometry, shared_ptr<VisualSceneGeometry>>
+		(m, "VisualSceneGeometry")
 		.def("Draw", &VisualSceneGeometry::DrawScene)
 		;
 
-    bp::def("SetBackGroundColor", &VisualSceneGeometry::SetBackGroundColor);
+    m.def("SetBackGroundColor", &VisualSceneGeometry::SetBackGroundColor);
 
-	bp::def("VS", FunctionPointer
-		([](CSGeometry & geom)
+	m.def("VS",
+		[](CSGeometry & geom)
 	{
 		geom.FindIdenticSurfaces(1e-6);
 		geom.CalcTriangleApproximation(0.01, 20);
@@ -518,17 +510,19 @@ DLL_HEADER void ExportCSGVis()
 
 		vs->SetGeometry(&geom);
 		return vs;
-	}));
+	});
 
-	bp::def("MouseMove", FunctionPointer
-		([](VisualSceneGeometry &vsgeom, int oldx, int oldy, int newx, int newy, char mode)
+	m.def("MouseMove",
+		[](VisualSceneGeometry &vsgeom, int oldx, int oldy, int newx, int newy, char mode)
 	{
 		vsgeom.MouseMove(oldx, oldy, newx, newy, mode);
-	}));
+	});
 }
-BOOST_PYTHON_MODULE(libcsgvis) 
-{
-	ExportCSGVis();  
+
+PYBIND11_PLUGIN(libcsgvis) {
+  py::module m("csg", "pybind csg");
+  ExportCSGVis(m);
+  return m.ptr();
 }
 #endif
 
diff --git a/libsrc/general/ngpython.hpp b/libsrc/general/ngpython.hpp
index 2ccad390..7c2d98a5 100644
--- a/libsrc/general/ngpython.hpp
+++ b/libsrc/general/ngpython.hpp
@@ -1,36 +1,42 @@
 #ifdef NG_PYTHON
 
-#include <boost/python.hpp>
-namespace bp = boost::python;
+#include <pybind11/pybind11.h>
+#include <pybind11/operators.h>
+namespace py = pybind11;
 #include <iostream>
+#include <sstream>
+
+namespace pybind11 {
+template<typename T>
+bool CheckCast( py::handle obj ) {
+  try{
+    obj.cast<T>();
+    return true;
+  }
+  catch (py::cast_error &e) {
+    return false;
+  }
+}
+
+
+template <typename T>
+struct extract
+{
+  py::handle obj;
+  extract( py::handle aobj ) : obj(aobj) {}
+
+  bool check() { return CheckCast<T>(obj); }
+  T operator()() { return obj.cast<T>(); }
+};
+}
+
+struct NGDummyArgument {};
+
+inline void NOOP_Deleter(void *) { ; }
 
 namespace netgen
 {
 
-	class ModuleScope  {
-		bp::scope *local_scope;
-	public:
-		ModuleScope(const std::string name) : local_scope(nullptr) {
-			std::string nested_name = name;
-			if (bp::scope())
-				nested_name = bp::extract<std::string>(bp::scope().attr("__name__") + "." + name);
-
-			bp::object module(bp::handle<>(bp::borrowed(PyImport_AddModule(nested_name.c_str()))));
-
-			std::cout << "exporting " << nested_name << std::endl;
-			bp::object parent = bp::scope() ? bp::scope() : bp::import("__main__");
-			parent.attr(name.c_str()) = module;
-
-			local_scope = new bp::scope(module);
-		}
-
-		~ModuleScope() {
-			if (local_scope)
-				delete (local_scope);
-		}
-		
-	};
-
   //////////////////////////////////////////////////////////////////////
   // Lambda to function pointer conversion
   template <typename Function>
diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp
index efc46ccf..3d8930d8 100644
--- a/libsrc/geom2d/python_geom2d.cpp
+++ b/libsrc/geom2d/python_geom2d.cpp
@@ -1,41 +1,34 @@
 #ifdef NG_PYTHON
 
-#include <boost/python.hpp>
 #include <../general/ngpython.hpp>
 
 #include <meshing.hpp>
 #include <geometry2d.hpp>
 
 using namespace netgen;
-namespace bp = boost::python;
-
-#if defined(_MSC_FULL_VER) && _MSC_FULL_VER == 190024213
-namespace boost { template<> const volatile SplineGeometry2d* get_pointer(const volatile SplineGeometry2d* p) { return p; } }
-#endif
 
 namespace netgen
 {
   extern std::shared_ptr<NetgenGeometry> ng_geometry;
 }
 
-DLL_HEADER void ExportGeom2d() 
-{
-  ModuleScope module("geom2d");
 
-  bp::class_<SplineGeometry2d, shared_ptr<SplineGeometry2d>, boost::noncopyable>
-    ("SplineGeometry",
+DLL_HEADER void ExportGeom2d(py::module &m) 
+{
+
+  py::class_<SplineGeometry2d, shared_ptr<SplineGeometry2d>>
+    (m, "SplineGeometry",
      "a 2d boundary representation geometry model by lines and splines")
-    .def("__init__", bp::make_constructor 
-         (FunctionPointer
-          ([](const string & filename)
+    .def(py::init<>())
+    .def("__init__",
+          [](SplineGeometry2d *instance, const string & filename)
            {
              cout << "load geometry";
              ifstream ist(filename);
-             auto geom = make_shared<SplineGeometry2d>();
-             geom->Load (filename.c_str());
-             ng_geometry = geom;
-             return geom;
-           })))
+             new (instance) SplineGeometry2d();
+             instance->Load (filename.c_str());
+             ng_geometry = shared_ptr<SplineGeometry2d>(instance, NOOP_Deleter);
+           })
     
 	.def("Load",&SplineGeometry2d::Load)
     .def("AppendPoint", FunctionPointer
@@ -50,17 +43,17 @@ DLL_HEADER void ExportGeom2d()
             self.geompoints.Append(gp);
             return self.geompoints.Size()-1;
 	  }),
-         (bp::arg("self"), bp::arg("x"), bp::arg("y"), bp::arg("maxh") = 1e99, bp::arg("hpref")=false))
-    .def("Append", FunctionPointer([](SplineGeometry2d &self, bp::list segment, int leftdomain, int rightdomain,
-                                      bp::object bc, bp::object copy, double maxh, bool hpref)
+         py::arg("x"), py::arg("y"), py::arg("maxh") = 1e99, py::arg("hpref")=false)
+    .def("Append", FunctionPointer([](SplineGeometry2d &self, py::list segment, int leftdomain, int rightdomain,
+                                      py::object bc, py::object copy, double maxh, bool hpref)
 	  {
-            bp::extract<std::string> segtype(segment[0]);
+            py::extract<std::string> segtype(segment[0]);
             
             SplineSegExt * seg;
             if (segtype().compare("line") == 0)
               {
-                bp::extract<int> point_index1(segment[1]);
-                bp::extract<int> point_index2(segment[2]);
+                py::extract<int> point_index1(segment[1]);
+                py::extract<int> point_index2(segment[2]);
                 //point_index1.check()
                 
                 LineSeg<2> * l = new LineSeg<2>(self.GetPoint(point_index1()), self.GetPoint(point_index2()));
@@ -68,9 +61,9 @@ DLL_HEADER void ExportGeom2d()
               }
             else if (segtype().compare("spline3") == 0)
               {
-                bp::extract<int> point_index1(segment[1]);
-                bp::extract<int> point_index2(segment[2]);
-                bp::extract<int> point_index3(segment[3]);
+                py::extract<int> point_index1(segment[1]);
+                py::extract<int> point_index2(segment[2]);
+                py::extract<int> point_index3(segment[3]);
                 
                 SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(point_index1()), self.GetPoint(point_index2()), self.GetPoint(point_index3()));
                 seg = new SplineSegExt(*seg3);
@@ -86,14 +79,14 @@ DLL_HEADER void ExportGeom2d()
             seg->hpref_right = hpref;
             seg->reffak = 1;
             seg->copyfrom = -1;
-            if (bp::extract<int>(copy).check())
-              seg->copyfrom = bp::extract<int>(copy)()+1;
+            if (py::extract<int>(copy).check())
+              seg->copyfrom = py::extract<int>(copy)()+1;
               
-            if (bp::extract<int>(bc).check())
-              seg->bc = bp::extract<int>(bc)();
-            else if (bp::extract<string>(bc).check())
+            if (py::extract<int>(bc).check())
+              seg->bc = py::extract<int>(bc)();
+            else if (py::extract<string>(bc).check())
               {
-                string bcname = bp::extract<string>(bc)();
+                string bcname = py::extract<string>(bc)();
                 int bcnum = self.GetBCNumber(bcname);
                 if (bcnum == 0)
                   bcnum = self.AddBCName(bcname);
@@ -103,25 +96,25 @@ DLL_HEADER void ExportGeom2d()
               seg->bc = self.GetNSplines()+1;
             self.AppendSegment(seg);
             return self.GetNSplines()-1;
-	  }), (bp::arg("self"), bp::arg("point_indices"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0,
-               bp::arg("bc")=bp::object(), bp::arg("copy")=bp::object(), bp::arg("maxh")=1e99, bp::arg("hpref")=false
-               ))
+	  }), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = 0,
+               py::arg("bc")=NGDummyArgument(), py::arg("copy")=NGDummyArgument(), py::arg("maxh")=1e99, py::arg("hpref")=false
+               )
 
     
-    .def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, bp::list point_indices, int leftdomain, int rightdomain)
+    .def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, py::list point_indices, int leftdomain, int rightdomain)
                                           {
-		  int npts = bp::len(point_indices);
+		  int npts = py::len(point_indices);
 		  SplineSegExt * seg;
-		  //int a = bp::extract<int>(point_indices[0]);
+		  //int a = py::extract<int>(point_indices[0]);
 		  if (npts == 2)
 		  {
-			  LineSeg<2> * l = new LineSeg<2>(self.GetPoint(bp::extract<int>(point_indices[0])), self.GetPoint(bp::extract<int>(point_indices[1])));
+			  LineSeg<2> * l = new LineSeg<2>(self.GetPoint(py::extract<int>(point_indices[0])()), self.GetPoint(py::extract<int>(point_indices[1])()));
 			  seg = new SplineSegExt(*l);
 			  
 		  }
 		  else if (npts == 3)
 		  {
-			  SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(bp::extract<int>(point_indices[0])), self.GetPoint(bp::extract<int>(point_indices[1])), self.GetPoint(bp::extract<int>(point_indices[2])));
+			  SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(py::extract<int>(point_indices[0])()), self.GetPoint(py::extract<int>(point_indices[1])()), self.GetPoint(py::extract<int>(point_indices[2])()));
 			  seg = new SplineSegExt(*seg3);
 
 		  }
@@ -131,7 +124,7 @@ DLL_HEADER void ExportGeom2d()
 		  seg->reffak = 1;
 		  seg->copyfrom = -1;
 		  self.AppendSegment(seg);
-		}), (bp::arg("self"), bp::arg("point_indices"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0) )
+		}), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = 0)
 	//.def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, int point_index1, int point_index2)//, int leftdomain, int rightdomain)
 	//  {
 	//	  LineSeg<2> * l = new LineSeg<2>(self.GetPoint(point_index1), self.GetPoint(point_index2));
@@ -143,7 +136,7 @@ DLL_HEADER void ExportGeom2d()
 	//	  seg->copyfrom = -1;
 
 	//	  self.AppendSegment(seg);
-	//  }))//, (bp::arg("self"), bp::arg("point_index1"), bp::arg("point_index2"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0) )
+	//  }))//, (py::arg("self"), py::arg("point_index1"), py::arg("point_index2"), py::arg("leftdomain") = 1, py::arg("rightdomain") = 0) )
 	//.def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, int point_index1, int point_index2, int point_index3)//, int leftdomain, int rightdomain)
 	//  {
 	//	  SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(point_index1), self.GetPoint(point_index2), self.GetPoint(point_index3));
@@ -154,7 +147,7 @@ DLL_HEADER void ExportGeom2d()
 	//	  seg->reffak = 1;
 	//	  seg->copyfrom = -1;
 	//	  self.AppendSegment(seg);
-	//  }))//, (bp::arg("self"), bp::arg("point_index1"), bp::arg("point_index2"), bp::arg("point_index3"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0 ) )
+	//  }))//, (py::arg("self"), py::arg("point_index1"), py::arg("point_index2"), py::arg("point_index3"), py::arg("leftdomain") = 1, py::arg("rightdomain") = 0 ) )
 
 
     .def("SetMaterial", &SplineGeometry2d::SetMaterial)
@@ -169,22 +162,22 @@ DLL_HEADER void ExportGeom2d()
 		  Box<2> box(self.GetBoundingBox());
 		  double xdist = box.PMax()(0) - box.PMin()(0);
 		  double ydist = box.PMax()(1) - box.PMin()(1);
-		  bp::tuple xlim = bp::make_tuple(box.PMin()(0) - 0.1*xdist, box.PMax()(0) + 0.1*xdist);
-		  bp::tuple ylim = bp::make_tuple(box.PMin()(1) - 0.1*ydist, box.PMax()(1) + 0.1*ydist);
+		  py::tuple xlim = py::make_tuple(box.PMin()(0) - 0.1*xdist, box.PMax()(0) + 0.1*xdist);
+		  py::tuple ylim = py::make_tuple(box.PMin()(1) - 0.1*ydist, box.PMax()(1) + 0.1*ydist);
 
-		  bp::list xpoints, ypoints;
+		  py::list xpoints, ypoints;
 
 		  for (int i = 0; i < self.splines.Size(); i++)
 		  {
-			  bp::list xp, yp;
+			  py::list xp, yp;
 			  if (self.splines[i]->GetType().compare("line")==0)
 			  {
 				  GeomPoint<2> p1 = self.splines[i]->StartPI();
 				  GeomPoint<2> p2 = self.splines[i]->EndPI();
-				  xp.append(p1(0));
-				  xp.append(p2(0));
-				  yp.append(p1(1));
-				  yp.append(p2(1));
+				  xp.append(py::cast(p1(0)));
+				  xp.append(py::cast(p2(0)));
+				  yp.append(py::cast(p1(1)));
+				  yp.append(py::cast(p2(1)));
 			  }
 			  else if (self.splines[i]->GetType().compare("spline3")==0)
 			  {
@@ -194,37 +187,37 @@ DLL_HEADER void ExportGeom2d()
 				  for (int j = 0; j <= n; j++)
 				  {
 					  GeomPoint<2> point = self.splines[i]->GetPoint(j*1./n);
-					  xp.append(point(0));
-					  yp.append(point(1));
+					  xp.append(py::cast(point(0)));
+					  yp.append(py::cast(point(1)));
 				  }
 			  }
 			  else
 			  {
 				  cout << "spline is neither line nor spline3" << endl;
 			  }
-			  xpoints.append(xp);
-			  ypoints.append(yp);
+			  xpoints.append(py::cast(xp));
+			  ypoints.append(py::cast(yp));
 				  
 		  }
-		  return bp::tuple(bp::make_tuple(xlim, ylim, xpoints, ypoints));
+		  return py::tuple(py::make_tuple(xlim, ylim, xpoints, ypoints));
 
 	  }))
 	.def("PointData", FunctionPointer([](SplineGeometry2d &self)
 	  {
-		  bp::list xpoints, ypoints, pointindex;
+		  py::list xpoints, ypoints, pointindex;
 		  
 		  for (int i = 0; i < self.geompoints.Size(); i++)
 		  {
-			  pointindex.append(i);
-			  xpoints.append(self.geompoints[i][0]);
-			  ypoints.append(self.geompoints[i][1]);
+			  pointindex.append(py::cast(i));
+			  xpoints.append(py::cast(self.geompoints[i][0]));
+			  ypoints.append(py::cast(self.geompoints[i][1]));
 		  }
-		  return bp::tuple(bp::make_tuple(xpoints, ypoints, pointindex));
+		  return py::tuple(py::make_tuple(xpoints, ypoints, pointindex));
 		  
 	  }))
 	.def("SegmentData", FunctionPointer([](SplineGeometry2d &self)
 	  {
-		  bp::list leftpoints, rightpoints, leftdom, rightdom;
+		  py::list leftpoints, rightpoints, leftdom, rightdom;
 
 		  for (int i = 0; i < self.splines.Size(); i++)
 		  {
@@ -234,13 +227,13 @@ DLL_HEADER void ExportGeom2d()
 			  normal(0) = normal(1);
 			  normal(1) = -temp;
 
-			  leftdom.append(self.GetSpline(i).leftdom);
-			  rightdom.append(self.GetSpline(i).rightdom);
+			  leftdom.append(py::cast(self.GetSpline(i).leftdom));
+			  rightdom.append(py::cast(self.GetSpline(i).rightdom));
 
-			  rightpoints.append(bp::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0));
-			  leftpoints.append(bp::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0));
+			  rightpoints.append(py::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0));
+			  leftpoints.append(py::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0));
 		  }
-		  return bp::tuple(bp::make_tuple(leftpoints, rightpoints, leftdom, rightdom));
+		  return py::tuple(py::make_tuple(leftpoints, rightpoints, leftdom, rightdom));
 
 	  }))
 	.def("Print", FunctionPointer([](SplineGeometry2d &self)
@@ -273,8 +266,10 @@ DLL_HEADER void ExportGeom2d()
   
 }
 
-BOOST_PYTHON_MODULE(libgeom2d) {
-	ExportGeom2d();
+PYBIND11_PLUGIN(libgeom2d) {
+  py::module m("geom2d", "pybind geom2d");
+  ExportGeom2d(m);
+  return m.ptr();
 }
 
 #endif
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index fc5acc62..8e865f91 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1,7 +1,5 @@
 #ifdef NG_PYTHON
 
-#include <boost/python.hpp>
-#include <boost/python/slice.hpp>
 #include <../general/ngpython.hpp>
 
 #include <mystdlib.h>
@@ -12,7 +10,6 @@
 
 
 using namespace netgen;
-namespace bp = boost::python;
 
 namespace netgen
 {
@@ -21,23 +18,19 @@ namespace netgen
 
 
 template <typename T, int BASE = 0, typename TIND = int>
-void ExportArray ()
+void ExportArray (py::module &m)
 {
   string name = string("Array_") + typeid(T).name();
-  bp::class_<Array<T,BASE,TIND>,boost::noncopyable>(name.c_str())
-    .def ("__len__", &Array<T,BASE,TIND>::Size)
+  py::class_<Array<T,BASE,TIND>>(m, name.c_str())
+    .def ("__len__", [] ( Array<T,BASE,TIND> &self ) { return self.Size(); } )
     .def ("__getitem__", 
           FunctionPointer ([](Array<T,BASE,TIND> & self, TIND i) -> T&
                            {
                              if (i < BASE || i >= BASE+self.Size())
-                               bp::exec("raise IndexError()\n");
+                               throw py::index_error();
                              return self[i];
                            }),
-          bp::return_value_policy<bp::reference_existing_object>())
-
-    .def ("__iter__", 
-          bp::range (FunctionPointer([](Array<T,BASE,TIND> & self) { return &self[BASE]; }),
-                     FunctionPointer([](Array<T,BASE,TIND> & self) { return &self[BASE+self.Size()]; })))
+          py::return_value_policy::reference)
 
     ;
 }
@@ -49,45 +42,46 @@ void TranslateException (const NgException & ex)
 }
 
 
-DLL_HEADER void ExportNetgenMeshing() 
+DLL_HEADER void ExportNetgenMeshing(py::module &m) 
 {
   
-  ModuleScope module("meshing");
-
-  bp::register_exception_translator<NgException>(&TranslateException);
   
-  bp::class_<PointIndex>("PointId", bp::init<int>())
+  py::class_<PointIndex>(m, "PointId")
+    .def(py::init<int>())
     .def("__repr__", &ToString<PointIndex>)
     .def("__str__", &ToString<PointIndex>)
-    .add_property("nr", &PointIndex::operator int)
+    .def_property_readonly("nr", &PointIndex::operator int)
     .def("__eq__" , FunctionPointer( [](PointIndex &self, PointIndex &other)
                   { return static_cast<int>(self)==static_cast<int>(other); }) )
     .def("__hash__" , FunctionPointer( [](PointIndex &self ) { return static_cast<int>(self); }) )
     ;
 
-  bp::class_<ElementIndex>("ElementId3D", bp::init<int>())
+  py::class_<ElementIndex>(m, "ElementId3D")
+    .def(py::init<int>())
     .def("__repr__", &ToString<ElementIndex>)
     .def("__str__", &ToString<ElementIndex>)
-    .add_property("nr", &ElementIndex::operator int)
+    .def_property_readonly("nr", &ElementIndex::operator int)
     .def("__eq__" , FunctionPointer( [](ElementIndex &self, ElementIndex &other)
                   { return static_cast<int>(self)==static_cast<int>(other); }) )
     .def("__hash__" , FunctionPointer( [](ElementIndex &self ) { return static_cast<int>(self); }) )
     ;
 
 
-  bp::class_<SurfaceElementIndex>("ElementId2D", bp::init<int>())
+  py::class_<SurfaceElementIndex>(m, "ElementId2D")
+    .def(py::init<int>())
     .def("__repr__", &ToString<SurfaceElementIndex>)
     .def("__str__", &ToString<SurfaceElementIndex>)
-    .add_property("nr", &SurfaceElementIndex::operator int)
+    .def_property_readonly("nr", &SurfaceElementIndex::operator int)
     .def("__eq__" , FunctionPointer( [](SurfaceElementIndex &self, SurfaceElementIndex &other)
                   { return static_cast<int>(self)==static_cast<int>(other); }) )
     .def("__hash__" , FunctionPointer( [](SurfaceElementIndex &self ) { return static_cast<int>(self); }) )
     ;
 
-  bp::class_<SegmentIndex>("ElementId1D", bp::init<int>())
+  py::class_<SegmentIndex>(m, "ElementId1D")
+    .def(py::init<int>())
     .def("__repr__", &ToString<SegmentIndex>)
     .def("__str__", &ToString<SegmentIndex>)
-    .add_property("nr", &SegmentIndex::operator int)
+    .def_property_readonly("nr", &SegmentIndex::operator int)
     .def("__eq__" , FunctionPointer( [](SegmentIndex &self, SegmentIndex &other)
                   { return static_cast<int>(self)==static_cast<int>(other); }) )
     .def("__hash__" , FunctionPointer( [](SegmentIndex &self ) { return static_cast<int>(self); }) )
@@ -96,176 +90,162 @@ DLL_HEADER void ExportNetgenMeshing()
 
 
   /*  
-  bp::class_<Point<3>> ("Point")
-    .def(bp::init<double,double,double>())
+  py::class_<Point<3>> ("Point")
+    .def(py::init<double,double,double>())
     ;
   */
 
-  bp::class_<MeshPoint /* ,bp::bases<Point<3>> */ >("MeshPoint")
-    .def(bp::init<Point<3>>())
+  py::class_<MeshPoint /* ,py::bases<Point<3>> */ >(m, "MeshPoint")
+    .def(py::init<Point<3>>())
     .def("__str__", &ToString<MeshPoint>)
     .def("__repr__", &ToString<MeshPoint>)
-    .add_property("p", FunctionPointer([](const MeshPoint & self)
-                                        {
-                                          bp::list l;
-                                          l.append ( self[0] );
-                                          l.append ( self[1] );
-                                          l.append ( self[2] );
-                                          return bp::tuple(l);
+    .def_property_readonly("p", FunctionPointer([](const MeshPoint & self)
+                                       {
+                                         py::list l;
+                                         l.append ( py::cast(self[0]) );
+                                         l.append ( py::cast(self[1]) );
+                                         l.append ( py::cast(self[2]) );
+                                         return py::tuple(l);
                                        }))
     .def("__getitem__", FunctionPointer([](const MeshPoint & self, int index) {
 	  if(index<0 || index>2)
-	    bp::exec("raise IndexError()\n");
+              throw py::index_error();
 	  return self[index];
 	}))
     ;
   
-  bp::class_<Element>("Element3D")
-    .def("__init__", bp::make_constructor
-         (FunctionPointer ([](int index, bp::list vertices)
+  py::class_<Element>(m, "Element3D")
+    .def("__init__", [](Element *instance, int index, py::list vertices)
                            {
-                             if (bp::len(vertices) == 4)
+                             if (py::len(vertices) == 4)
                                {
-                                 Element * tmp = new Element(TET);
+                                 new (instance) Element(TET);
                                  for (int i = 0; i < 4; i++)
-                                   (*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
-                                 tmp->SetIndex(index);
-                                 return tmp;
+                                   (*instance)[i] = py::extract<PointIndex>(vertices[i])();
+                                 instance->SetIndex(index);
                                }
-                             if (bp::len(vertices) == 6)
+                             if (py::len(vertices) == 6)
                                {
-                                 Element * tmp = new Element(PRISM);
+                                 new (instance) Element(PRISM);
                                  for (int i = 0; i < 6; i++)
-                                   (*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
-                                 tmp->SetIndex(index);
-                                 return tmp;
+                                   (*instance)[i] = py::extract<PointIndex>(vertices[i])();
+                                 instance->SetIndex(index);
                                }
-                             if (bp::len(vertices) == 8)
+                             if (py::len(vertices) == 8)
                                {
-                                 Element * tmp = new Element(HEX);
+                                 new (instance) Element(HEX);
                                  for (int i = 0; i < 8; i++)
-                                   (*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
-                                 tmp->SetIndex(index);
-                                 return tmp;
+                                   (*instance)[i] = py::extract<PointIndex>(vertices[i])();
+                                 instance->SetIndex(index);
                                }
                              throw NgException ("cannot create element");                             
-                           }),
-          bp::default_call_policies(),        // need it to use arguments
-          (bp::arg("index")=1,bp::arg("vertices"))),
+                           },
+          py::arg("index")=1,py::arg("vertices"),
          "create volume element"
          )
     .def("__repr__", &ToString<Element>)
-    .add_property("index", &Element::GetIndex, &Element::SetIndex)
-    .add_property("vertices", 
-                  FunctionPointer ([](const Element & self) -> bp::list
+    .def_property("index", &Element::GetIndex, &Element::SetIndex)
+    .def_property_readonly("vertices", 
+                  FunctionPointer ([](const Element & self) -> py::list
                                    {
-                                     bp::list li;
+                                     py::list li;
                                      for (int i = 0; i < self.GetNV(); i++)
-                                       li.append (self[i]);
+                                       li.append (py::cast(self[i]));
                                      return li;
                                    }))
     ;
 
-  bp::class_<Element2d>("Element2D")
-    .def("__init__", bp::make_constructor
-         (FunctionPointer ([](int index, bp::list vertices)
+  py::class_<Element2d>(m, "Element2D")
+    .def("__init__", 
+         [](Element2d *instance, int index, py::list vertices)
                            {
-                             if (bp::len(vertices) == 3)
+                             if (py::len(vertices) == 3)
                                {
-                                 Element2d * tmp = new Element2d(TRIG);
+                                 new (instance) Element2d(TRIG);
                                  for (int i = 0; i < 3; i++)
-                                   (*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
-                                 tmp->SetIndex(index);
-                                 return tmp;
+                                   (*instance)[i] = py::extract<PointIndex>(vertices[i])();
+                                 instance->SetIndex(index);
                                }
                              else
                                {
-                                 Element2d * tmp = new Element2d(QUAD);
+                                 new (instance) Element2d(QUAD);
                                  for (int i = 0; i < 4; i++)
-                                   (*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
-                                 tmp->SetIndex(index);
-                                 return tmp;
+                                   (*instance)[i] = py::extract<PointIndex>(vertices[i])();
+                                 instance->SetIndex(index);
                                }
                                
-                           }),
-          bp::default_call_policies(),        // need it to use arguments
-          (bp::arg("index")=1,bp::arg("vertices"))),
+                           },
+          py::arg("index")=1,py::arg("vertices"),
          "create surface element"
          )
-    .add_property("index", &Element2d::GetIndex, &Element2d::SetIndex)
-    .add_property("vertices",
-                  FunctionPointer([](const Element2d & self) -> bp::list
+    .def_property("index", &Element2d::GetIndex, &Element2d::SetIndex)
+    .def_property_readonly("vertices",
+                  FunctionPointer([](const Element2d & self) -> py::list
                                   {
-                                    bp::list li;
+                                    py::list li;
                                     for (int i = 0; i < self.GetNV(); i++)
-                                      li.append(self[i]);
+                                      li.append(py::cast(self[i]));
                                     return li;
                                   }))
     ;
 
-  bp::class_<Segment>("Element1D")
-    .def("__init__", bp::make_constructor
-         (FunctionPointer ([](bp::list vertices, bp::list surfaces, int index)
+  py::class_<Segment>(m, "Element1D")
+    .def("__init__",
+         [](Segment *instance, py::list vertices, py::list surfaces, int index)
                            {
-                             Segment * tmp = new Segment;
+                             new (instance) Segment();
                              for (int i = 0; i < 2; i++)
-                               (*tmp)[i] = bp::extract<PointIndex>(vertices[i]);
-                             tmp -> si = index;
+                               (*instance)[i] = py::extract<PointIndex>(vertices[i])();
+                             instance -> si = index;
                              if (len(surfaces))
                                {
-                                 tmp->surfnr1 = bp::extract<int>(surfaces[0]);
-                                 tmp->surfnr2 = bp::extract<int>(surfaces[1]);
+                                 instance->surfnr1 = py::extract<int>(surfaces[0])();
+                                 instance->surfnr2 = py::extract<int>(surfaces[1])();
                                }
-                             return tmp;
-                           }),
-          bp::default_call_policies(),      
-          (bp::arg("vertices"),
-           bp::arg("surfaces")=bp::list(),
-           bp::arg("index")=1
-           )),
+                           },
+          py::arg("vertices"),
+           py::arg("surfaces")=py::list(),
+           py::arg("index")=1,
          "create segment element"
          )
     .def("__repr__", &ToString<Segment>)
-    .add_property("vertices", 
-                  FunctionPointer ([](const Segment & self) -> bp::list
+    .def_property_readonly("vertices", 
+                  FunctionPointer ([](const Segment & self) -> py::list
                                    {
-                                     bp::list li;
+                                     py::list li;
                                      for (int i = 0; i < 2; i++)
-                                       li.append (self[i]);
+                                       li.append (py::cast(self[i]));
                                      return li;
                                    }))
-    .add_property("surfaces", 
-                  FunctionPointer ([](const Segment & self) -> bp::list
+    .def_property_readonly("surfaces", 
+                  FunctionPointer ([](const Segment & self) -> py::list
                                    {
-                                     bp::list li;
-                                     li.append (self.surfnr1);
-                                     li.append (self.surfnr2);
+                                     py::list li;
+                                     li.append (py::cast(self.surfnr1));
+                                     li.append (py::cast(self.surfnr2));
                                      return li;
                                    }))
     ;
 
 
-  bp::class_<Element0d>("Element0D")
-    .def("__init__", bp::make_constructor
-         (FunctionPointer ([](PointIndex vertex, int index)
+  py::class_<Element0d>(m, "Element0D")
+    .def("__init__",
+         [](Element0d *instance, PointIndex vertex, int index)
                            {
-                             Element0d * tmp = new Element0d;
-                             tmp->pnum = vertex;
-                             tmp->index = index;
-                             return tmp;
-                           }),
-          bp::default_call_policies(),      
-          (bp::arg("vertex"),
-           bp::arg("index")=1
-           )),
+                             new (instance) Element0d;
+                             instance->pnum = vertex;
+                             instance->index = index;
+                           },
+         py::arg("vertex"),
+         py::arg("index")=1,
          "create point element"
          )
     .def("__repr__", &ToString<Element0d>)
-    .add_property("vertices", 
-                  FunctionPointer ([](const Element0d & self) -> bp::list
+    .def_property_readonly("vertices", 
+                  FunctionPointer ([](const Element0d & self) -> py::list
                                    {
-                                     bp::list li;
-                                     li.append (self.pnum);
+                                     py::list li;
+                                     li.append (py::cast(self.pnum));
                                      return li;
                                    }))
     ;
@@ -274,70 +254,61 @@ DLL_HEADER void ExportNetgenMeshing()
   
 
 
-  bp::class_<FaceDescriptor>("FaceDescriptor")
-    .def(bp::init<const FaceDescriptor&>())
-    .def("__init__", bp::make_constructor
-         (FunctionPointer ([](int surfnr, int domin, int domout, int bc)
+  py::class_<FaceDescriptor>(m, "FaceDescriptor")
+    .def(py::init<const FaceDescriptor&>())
+    .def("__init__", 
+         [](FaceDescriptor *instance, int surfnr, int domin, int domout, int bc)
                            {
-                             auto fd = new FaceDescriptor();
-                             fd->SetSurfNr(surfnr);
-                             fd->SetDomainIn(domin);
-                             fd->SetDomainOut(domout);
-                             fd->SetBCProperty(bc);
-                             return fd;
-                           }),
-          bp::default_call_policies(),        // need it to use arguments
-          (bp::arg("surfnr")=1, 
-           bp::arg("domin")=1,
-           bp::arg("domout")=0,
-           bp::arg("bc")=0
-           )),
+                             new (instance) FaceDescriptor();
+                             instance->SetSurfNr(surfnr);
+                             instance->SetDomainIn(domin);
+                             instance->SetDomainOut(domout);
+                             instance->SetBCProperty(bc);
+                           },
+           py::arg("surfnr")=1, 
+           py::arg("domin")=1,
+           py::arg("domout")=0,
+           py::arg("bc")=0,
          "create facedescriptor")
     .def("__str__", &ToString<FaceDescriptor>)
     .def("__repr__", &ToString<FaceDescriptor>)
-    .add_property("surfnr", &FaceDescriptor::SurfNr, &FaceDescriptor::SetSurfNr)
-    .add_property("domin", &FaceDescriptor::DomainIn, &FaceDescriptor::SetDomainIn)
-    .add_property("domout", &FaceDescriptor::DomainOut, &FaceDescriptor::SetDomainOut)
-    .add_property("bc", &FaceDescriptor::BCProperty, &FaceDescriptor::SetBCProperty)
-    .add_property("bcname", FunctionPointer ([](FaceDescriptor & self) -> string { return self.GetBCName(); }))
-    .def("SetSurfaceColor", FunctionPointer ([](FaceDescriptor & self, bp::list color )
+    .def_property("surfnr", &FaceDescriptor::SurfNr, &FaceDescriptor::SetSurfNr)
+    .def_property("domin", &FaceDescriptor::DomainIn, &FaceDescriptor::SetDomainIn)
+    .def_property("domout", &FaceDescriptor::DomainOut, &FaceDescriptor::SetDomainOut)
+    .def_property("bc", &FaceDescriptor::BCProperty, &FaceDescriptor::SetBCProperty)
+    .def_property_readonly("bcname", FunctionPointer ([](FaceDescriptor & self) -> string { return self.GetBCName(); }))
+    .def("SetSurfaceColor", [](FaceDescriptor & self, py::list color )
           {
             Vec3d c;
-            c.X() = bp::extract<double>(color[0])();
-            c.Y() = bp::extract<double>(color[1])();
-            c.Z() = bp::extract<double>(color[2])();
+            c.X() = py::extract<double>(color[0])();
+            c.Y() = py::extract<double>(color[1])();
+            c.Z() = py::extract<double>(color[2])();
             self.SetSurfColour(c);
-          }))
+          })
     ;
 
   
 
-  ExportArray<Element>();
-  ExportArray<Element2d>();
-  ExportArray<Segment>();
-  ExportArray<Element0d>();
-  ExportArray<MeshPoint,PointIndex::BASE,PointIndex>();
-  ExportArray<FaceDescriptor>();
-  ;
-  
-#if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100)  
-  bp::register_ptr_to_python<shared_ptr<Mesh>>();
-#endif
-  bp::class_<Mesh,shared_ptr<Mesh>,boost::noncopyable>("Mesh", bp::no_init)
-    // .def(bp::init<>("create empty mesh"))
+  ExportArray<Element>(m);
+  ExportArray<Element2d>(m);
+  ExportArray<Segment>(m);
+  ExportArray<Element0d>(m);
+  ExportArray<MeshPoint,PointIndex::BASE,PointIndex>(m);
+  ExportArray<FaceDescriptor>(m);
 
-    .def("__init__", bp::make_constructor 
-         (FunctionPointer ([](int dim)
+  py::implicitly_convertible< int, PointIndex>();
+  
+  py::class_<Mesh,shared_ptr<Mesh>>(m, "Mesh")
+    // .def(py::init<>("create empty mesh"))
+
+    .def("__init__",
+         [](Mesh *instance, int dim)
                            {
-                             auto mesh = make_shared<Mesh>();
-                             mesh->SetDimension(dim);
-                             return mesh;
-                           }),
-          bp::default_call_policies(),     // need it to use named arguments
-          (
-           bp::arg("dim")=3
-           )
-          ))
+                             new (instance) Mesh();
+                             instance->SetDimension(dim);
+                           },
+           py::arg("dim")=3
+          )
 
     
     .def("__str__", &ToString<Mesh>)
@@ -379,8 +350,8 @@ DLL_HEADER void ExportNetgenMeshing()
 	  }))
     // static_cast<void(Mesh::*)(const string & name)>(&Mesh::Load))
     .def("Save", static_cast<void(Mesh::*)(const string & name)const>(&Mesh::Save))
-    .def("Export", FunctionPointer
-         ([] (Mesh & self, string filename, string format)
+    .def("Export",
+         [] (Mesh & self, string filename, string format)
           {
             if (WriteUserFormat (format, self, *self.GetGeometry(), filename))
               {
@@ -392,35 +363,35 @@ DLL_HEADER void ExportNetgenMeshing()
                   err += string("'") + name + "'\n";
                 throw NgException (err);
               }
-          }),
-         (bp::arg("self"), bp::arg("filename"), bp::arg("format")))
+          },
+         py::arg("filename"), py::arg("format"))
     
-    .add_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
+    .def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
 
     .def("Elements3D", 
          static_cast<Array<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
-         bp::return_value_policy<bp::reference_existing_object>())
+         py::return_value_policy::reference)
 
     .def("Elements2D", 
          static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements),
-         bp::return_value_policy<bp::reference_existing_object>())
+         py::return_value_policy::reference)
 
     .def("Elements1D", 
          static_cast<Array<Segment>&(Mesh::*)()> (&Mesh::LineSegments),
-         bp::return_value_policy<bp::reference_existing_object>())
+         py::return_value_policy::reference)
 
     .def("Elements0D", FunctionPointer([] (Mesh & self) -> Array<Element0d>&
                                        {
                                          return self.pointelements;
                                        } ),
-         bp::return_value_policy<bp::reference_existing_object>())
+         py::return_value_policy::reference)
 
     .def("Points", 
          static_cast<Mesh::T_POINTS&(Mesh::*)()> (&Mesh::Points),
-         bp::return_value_policy<bp::reference_existing_object>())
+         py::return_value_policy::reference)
 
     .def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor),
-         bp::return_value_policy<bp::reference_existing_object>())
+         py::return_value_policy::reference)
     .def("GetNFaceDescriptors", &Mesh::GetNFD)
     
 
@@ -466,23 +437,22 @@ DLL_HEADER void ExportNetgenMeshing()
     .def ("GetMaterial", FunctionPointer([](Mesh & self, int domnr)
                                          { return string(self.GetMaterial(domnr)); }))
 
-    .def ("GenerateVolumeMesh", FunctionPointer
-          ([](Mesh & self, bp::object pymp)
+    .def ("GenerateVolumeMesh",
+          [](Mesh & self, py::object pymp)
            {
              cout << "generate vol mesh" << endl;
 
              MeshingParameters mp;
-             if (bp::extract<MeshingParameters>(pymp).check())
-               mp = bp::extract<MeshingParameters>(pymp)();
+             if (py::extract<MeshingParameters>(pymp).check())
+               mp = py::extract<MeshingParameters>(pymp)();
              else
                {
                  mp.optsteps3d = 5;
                }
              MeshVolume (mp, self);
              OptimizeVolume (mp, self);
-           }),
-          (bp::arg("self"), bp::arg("mp")=bp::object())
-          )
+           },
+          py::arg("mp")=NGDummyArgument())
 
    .def ("OptimizeVolumeMesh", FunctionPointer
          ([](Mesh & self)
@@ -507,18 +477,19 @@ DLL_HEADER void ExportNetgenMeshing()
              self.SetGeometry(geo);
            }))
 
-    .def ("SetGeometry", FunctionPointer
-          ([](Mesh & self, shared_ptr<SplineGeometry2d> geo)
-           {
-             self.SetGeometry(geo);
-           }))
+    // TODO: fix this dependency on libgeom2d.so
+//     .def ("SetGeometry", FunctionPointer
+//           ([](Mesh & self, shared_ptr<SplineGeometry2d> geo)
+//            {
+//              self.SetGeometry(geo);
+//            }))
 
     .def ("BuildSearchTree", &Mesh::BuildElementSearchTree)
 
     .def ("BoundaryLayer", FunctionPointer 
-          ([](Mesh & self, int bc, bp::list thicknesses, int volnr, bp::list materials)
+          ([](Mesh & self, int bc, py::list thicknesses, int volnr, py::list materials)
            {
-             int n = bp::len(thicknesses);
+             int n = py::len(thicknesses);
              BoundaryLayerParameters blp;
 
              for (int i = 1; i <= self.GetNFD(); i++)
@@ -537,9 +508,9 @@ DLL_HEADER void ExportNetgenMeshing()
              cout << "maxind = " << maxind << endl;
              for ( int i=0; i<n; i++ )
                {
-                 blp.heights.Append( bp::extract<double>(thicknesses[i])()) ;
+                 blp.heights.Append( py::extract<double>(thicknesses[i])()) ;
                  blp.new_matnrs.Append( maxind+1+i );
-                 self.SetMaterial (maxind+1+i, bp::extract<string>(materials[i])().c_str());
+                 self.SetMaterial (maxind+1+i, py::extract<string>(materials[i])().c_str());
                }
              blp.bulk_matnr = volnr;
              GenerateBoundaryLayer (self, blp);
@@ -577,28 +548,26 @@ DLL_HEADER void ExportNetgenMeshing()
   
 
   typedef MeshingParameters MP;
-  bp::class_<MP> ("MeshingParameters", bp::init<>())
-    .def("__init__", bp::make_constructor
-         (FunctionPointer ([](double maxh, bool quad_dominated, int optsteps2d, int optsteps3d)
+  py::class_<MP> (m, "MeshingParameters")
+    .def(py::init<>())
+    .def("__init__",
+         [](MP *instance, double maxh, bool quad_dominated, int optsteps2d, int optsteps3d)
                            {
-                             auto tmp = new MeshingParameters;
-                             tmp->maxh = maxh;
-                             tmp->quad = int(quad_dominated);
-                             tmp->optsteps2d = optsteps2d;
-                             tmp->optsteps3d = optsteps3d;
-                             return tmp;
-                           }),
-          bp::default_call_policies(),        // need it to use arguments
-          (
-           bp::arg("maxh")=1000,
-           bp::arg("quad_dominated")=false,
-           bp::arg("optsteps2d") = 3,
-           bp::arg("optsteps3d") = 3
-           )),
+                             new (instance) MeshingParameters;
+                             instance->maxh = maxh;
+                             instance->quad = int(quad_dominated);
+                             instance->optsteps2d = optsteps2d;
+                             instance->optsteps3d = optsteps3d;
+                           },
+           py::arg("maxh")=1000,
+           py::arg("quad_dominated")=false,
+           py::arg("optsteps2d") = 3,
+           py::arg("optsteps3d") = 3
+           ,
          "create meshing parameters"
           )
     .def("__str__", &ToString<MP>)
-    .add_property("maxh", 
+    .def_property("maxh", 
                   FunctionPointer ([](const MP & mp ) { return mp.maxh; }),
                   FunctionPointer ([](MP & mp, double maxh) { return mp.maxh = maxh; }))
     .def("RestrictH", FunctionPointer
@@ -606,17 +575,17 @@ DLL_HEADER void ExportNetgenMeshing()
           {
             mp.meshsize_points.Append ( MeshingParameters::MeshSizePoint (Point<3> (x,y,z), h));
           }),
-         (bp::arg("x"), bp::arg("y"), bp::arg("z"), bp::arg("h"))
+         py::arg("x"), py::arg("y"), py::arg("z"), py::arg("h")
          )
     ;
 
-  bp::def("SetTestoutFile", FunctionPointer ([] (const string & filename)
+  m.def("SetTestoutFile", FunctionPointer ([] (const string & filename)
                                              {
                                                delete testout;
                                                testout = new ofstream (filename);
                                              }));
 
-  bp::def("SetMessageImportance", FunctionPointer ([] (int importance)
+  m.def("SetMessageImportance", FunctionPointer ([] (int importance)
                                                    {
                                                      int old = printmessage_importance;
                                                      printmessage_importance = importance;
@@ -624,14 +593,11 @@ DLL_HEADER void ExportNetgenMeshing()
                                                    }));
 }
 
-
-
-BOOST_PYTHON_MODULE(libmesh) {
-  ExportNetgenMeshing();
+PYBIND11_PLUGIN(libmesh) {
+  py::module m("mesh", "pybind mesh");
+  ExportNetgenMeshing(m);
+  return m.ptr();
 }
-
-
-
 #endif
 
 
diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp
index bd0e5cb8..2a58e6e6 100644
--- a/libsrc/visualization/vsmesh.cpp
+++ b/libsrc/visualization/vsmesh.cpp
@@ -3529,18 +3529,11 @@ namespace netgen
 
 
 #ifdef NG_PYTHON
-#include <boost/python.hpp>
 #include <../general/ngpython.hpp>
-namespace bp = boost::python;
 
-#if defined(_MSC_FULL_VER) && _MSC_FULL_VER == 190024213
-namespace boost { template<> const volatile netgen::VisualSceneMesh* get_pointer(const volatile netgen::VisualSceneMesh* p) { return p; } }
-#endif
-
-DLL_HEADER void ExportMeshVis()
+DLL_HEADER void ExportMeshVis(py::module &m)
 {
   using namespace netgen;
-  ModuleScope myscope("meshvis");
   vispar.drawcolorbar = true;
   vispar.drawnetgenlogo = true;
   vispar.drawcoordinatecross = true;
@@ -3550,12 +3543,12 @@ DLL_HEADER void ExportMeshVis()
   vispar.drawtets = true;
   vispar.drawprisms = true;
   vispar.drawoutline = true;
-  bp::class_<VisualSceneMesh, shared_ptr<VisualSceneMesh>>
-    ("VisualSceneMesh", bp::no_init)
+  py::class_<VisualSceneMesh, shared_ptr<VisualSceneMesh>>
+    (m, "VisualSceneMesh")
     .def("Draw", &VisualSceneMesh::DrawScene)
     ;
 
-  bp::def("VS", FunctionPointer
+  m.def("VS", FunctionPointer
           ([](shared_ptr<Mesh> mesh)
            {
              auto vs = make_shared<VisualSceneMesh>();
@@ -3564,23 +3557,23 @@ DLL_HEADER void ExportMeshVis()
              return vs;
            }));
 
-  bp::def("MouseMove", FunctionPointer
+  m.def("MouseMove", FunctionPointer
           ([](VisualSceneMesh &vsmesh, int oldx, int oldy, int newx, int 
               newy, char mode)
            {
              vsmesh.MouseMove(oldx, oldy, newx, newy, mode);
            }));
-  bp::def("SelectFace", FunctionPointer
+  m.def("SelectFace", FunctionPointer
       ([] (int facenr) {
        vsmesh.SetSelectedFace(facenr);
        }));
-  bp::def("GetGlobalMesh", FunctionPointer
+  m.def("GetGlobalMesh", FunctionPointer
       ([] () {
        return vsmesh.GetMesh();
        }));
 }
-BOOST_PYTHON_MODULE(libvisual)
-{
-  ExportMeshVis();
-}
+// BOOST_PYTHON_MODULE(libvisual)
+// {
+//   ExportMeshVis();
+// }
 #endif
diff --git a/ng/netgenpy.cpp b/ng/netgenpy.cpp
index 9478dadf..32d9c694 100644
--- a/ng/netgenpy.cpp
+++ b/ng/netgenpy.cpp
@@ -1,7 +1,7 @@
 // a wrapper to load netgen-dll into python
 
 #include <iostream>
-#include <boost/python.hpp>
+#include <../general/ngpython.hpp>
 
 #ifdef WIN32
 #define DLL_HEADER __declspec(dllimport)
@@ -10,19 +10,26 @@
 #endif
 
 
-void DLL_HEADER ExportNetgenMeshing();
-void DLL_HEADER ExportMeshVis();
-void DLL_HEADER ExportCSG();
-void DLL_HEADER ExportCSGVis();
-void DLL_HEADER ExportGeom2d();
+void DLL_HEADER ExportNetgenMeshing(py::module &m);
+void DLL_HEADER ExportMeshVis(py::module &m);
+void DLL_HEADER ExportCSG(py::module &m);
+void DLL_HEADER ExportCSGVis(py::module &m);
+void DLL_HEADER ExportGeom2d(py::module &m);
 
-BOOST_PYTHON_MODULE(libngpy) 
+PYBIND11_PLUGIN(libngpy)
 {
-    ExportCSG();
-    ExportCSGVis();
-    ExportNetgenMeshing();
-    ExportMeshVis();
-    ExportGeom2d();
+    py::module ngpy("libngpy", "pybind netgen module");
+    py::module csg = ngpy.def_submodule("_csg", "pybind csg module");
+    ExportCSG(csg);
+    py::module csgvis = ngpy.def_submodule("csgvis", "pybind csgvis module");
+    ExportCSGVis(csgvis);
+    py::module meshing = ngpy.def_submodule("_meshing", "pybind meshing module");
+    ExportNetgenMeshing(meshing);
+    py::module meshvis = ngpy.def_submodule("meshvis", "pybind meshvis module");
+    ExportMeshVis(meshvis);
+    py::module geom2d = ngpy.def_submodule("_geom2d", "pybind geom2d module");
+    ExportGeom2d(geom2d);
+    return ngpy.ptr();
 }
 
 // Force linking libnglib to libnetgenpy
diff --git a/python/__init__.py b/python/__init__.py
index 1db01972..b854a7b1 100644
--- a/python/__init__.py
+++ b/python/__init__.py
@@ -10,10 +10,14 @@ if __platform.startswith('darwin'):
 
  
 # from libngpy import *
-import libngpy 
+import libngpy
+# from libngpy import *
 
-from . import csg
-from . import meshing
+# import libngpy 
+
+# from . import csg
+# from . import meshing
+# from . import geom2d
 
 del environ
 del path
diff --git a/python/csg.py b/python/csg.py
index cde89fc9..11f54147 100644
--- a/python/csg.py
+++ b/python/csg.py
@@ -1,7 +1,8 @@
-from libngpy.csg import *
+import libngpy
+from libngpy._csg import *
 import libngpy.csgvis as csgvis
 from libngpy.csgvis import MouseMove
-from libngpy.meshing import MeshingParameters
+from libngpy._meshing import MeshingParameters
 
 
 CSGeometry.VS = csgvis.VS
diff --git a/python/geom2d.py b/python/geom2d.py
index debf2983..e4c8d925 100644
--- a/python/geom2d.py
+++ b/python/geom2d.py
@@ -1,5 +1,5 @@
-from libngpy.geom2d import *
-from libngpy.meshing import *
+from libngpy._geom2d import *
+from libngpy._meshing import *
 
 tmp_generate_mesh = SplineGeometry.GenerateMesh
 
diff --git a/python/meshing.py b/python/meshing.py
index 0e3aca83..c9c8f511 100644
--- a/python/meshing.py
+++ b/python/meshing.py
@@ -1 +1 @@
-from libngpy.meshing import *
+from libngpy._meshing import *