diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp
index ec8759b6..b99ea2d5 100644
--- a/libsrc/csg/python_csg.cpp
+++ b/libsrc/csg/python_csg.cpp
@@ -623,6 +623,68 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
           })
          )
     .def_property_readonly ("ntlo", &CSGeometry::GetNTopLevelObjects)
+    .def("_visualizationData", [](shared_ptr<CSGeometry> csg_geo)
+         {
+           std::vector<float> vertices;
+           std::vector<int> trigs;
+           std::vector<float> normals;
+           std::vector<float> min = {std::numeric_limits<float>::max(),
+                                     std::numeric_limits<float>::max(),
+                                     std::numeric_limits<float>::max()};
+           std::vector<float> max = {std::numeric_limits<float>::lowest(),
+                                     std::numeric_limits<float>::lowest(),
+                                     std::numeric_limits<float>::lowest()};
+           std::vector<string> surfnames;
+           for (int i = 0; i < csg_geo->GetNSurf(); i++)
+             {
+               auto surf = csg_geo->GetSurface(i);
+               surfnames.push_back(surf->GetBCName());
+             }
+           csg_geo->FindIdenticSurfaces(1e-6);
+           csg_geo->CalcTriangleApproximation(0.01,100);
+           auto nto = csg_geo->GetNTopLevelObjects();
+           size_t np = 0;
+           size_t ntrig = 0;
+           for (int i = 0; i < nto; i++){
+             np += csg_geo->GetTriApprox(i)->GetNP();
+             ntrig += csg_geo->GetTriApprox(i)->GetNT();
+           }
+           vertices.reserve(np*3);
+           trigs.reserve(ntrig*4);
+           normals.reserve(np*3);
+           int offset_points = 0;
+           for (int i = 0; i < nto; i++)
+             {
+               auto triapprox = csg_geo->GetTriApprox(i);
+               for (int j = 0; j < triapprox->GetNP(); j++)
+                 for(int k = 0; k < 3; k++) {
+                   float val = triapprox->GetPoint(j)[k];
+                   vertices.push_back(val);
+                   min[k] = min2(min[k], val);
+                   max[k] = max2(max[k],val);
+                   normals.push_back(triapprox->GetNormal(j)[k]);
+                 }
+               for (int j = 0; j < triapprox->GetNT(); j++)
+                 {
+                   for(int k = 0; k < 3; k++)
+                     trigs.push_back(triapprox->GetTriangle(j)[k]+offset_points);
+                   trigs.push_back(triapprox->GetTriangle(j).SurfaceIndex());
+                 }
+               offset_points += triapprox->GetNP();
+             }
+           py::gil_scoped_acquire ac;
+           py::dict res;
+           py::list snames;
+           for(auto name : surfnames)
+             snames.append(py::cast(name));
+           res["vertices"] = MoveToNumpy(vertices);
+           res["triangles"] = MoveToNumpy(trigs);
+           res["normals"] = MoveToNumpy(normals);
+           res["surfnames"] = snames;
+           res["min"] = MoveToNumpy(min);
+           res["max"] = MoveToNumpy(max);
+           return res;
+         }, py::call_guard<py::gil_scoped_release>())
     ;
 
   m.def("GenerateMesh", FunctionPointer
diff --git a/libsrc/general/ngpython.hpp b/libsrc/general/ngpython.hpp
index c3931215..396c8f06 100644
--- a/libsrc/general/ngpython.hpp
+++ b/libsrc/general/ngpython.hpp
@@ -14,11 +14,21 @@ namespace pybind11 {
 // END EVIL HACK
 #include <pybind11/pybind11.h>
 #include <pybind11/operators.h>
+#include <pybind11/numpy.h>
 namespace py = pybind11;
 #include <iostream>
 #include <sstream>
 
 
+template <typename T>
+py::array MoveToNumpy(std::vector<T>& vec)
+{
+  auto newvec = new std::vector<T>();
+  std::swap(*newvec, vec);
+  auto capsule = py::capsule(newvec, [](void *v) { delete reinterpret_cast<std::vector<T>*>(v); });
+  return py::array(newvec->size(), newvec->data(), capsule);
+}
+
 namespace PYBIND11_NAMESPACE {
 template<typename T>
 bool CheckCast( py::handle obj ) {
diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp
index 5e570af7..f44b525e 100644
--- a/libsrc/geom2d/python_geom2d.cpp
+++ b/libsrc/geom2d/python_geom2d.cpp
@@ -198,6 +198,79 @@ DLL_HEADER void ExportGeom2d(py::module &m)
 		  return py::tuple(py::make_tuple(xlim, ylim, xpoints, ypoints));
 
 	  }))
+    .def("_visualizationData", [](SplineGeometry2d &self)
+         {
+           Box<2> box(self.GetBoundingBox());
+           double xdist = box.PMax()(0) - box.PMin()(0);
+           double ydist = box.PMax()(1) - box.PMin()(1);
+           py::dict data;
+           py::dict segment_data;
+           auto min_val = py::make_tuple(box.PMin()(0), box.PMin()(1),0);
+           auto max_val = py::make_tuple(box.PMax()(1),box.PMax()(1),0);
+           py::list vertices;
+           py::list domains;
+           py::list segment_points;
+           py::list segment_normals;
+           py::list leftdom;
+           py::list rightdom;
+           int max_bcnr = 0;
+           for(int i = 0; i < self.splines.Size(); i++)
+             {
+               std::vector<netgen::GeomPoint<2>> lst;
+               if (self.splines[i]->GetType().compare("line") == 0)
+                 lst = { self.splines[i]->StartPI(), self.splines[i]->EndPI() };
+               else if(self.splines[i]->GetType().compare("spline3") == 0)
+                 {
+                   double len = self.splines[i]->Length();
+                   int n = floor(len/(0.05*min(xdist,ydist)));
+                   lst.push_back(self.splines[i]->StartPI());
+                   for (int j = 1; j < n; j++){
+                     lst.push_back(self.splines[i]->GetPoint(j*1./n));
+                     lst.push_back(self.splines[i]->GetPoint(j*1./n));
+                   }
+                   lst.push_back(self.splines[i]->EndPI());
+                   }
+               else
+                 {
+                   throw NgException("Spline is neither line nor spline3");
+                 }
+               for (auto point : lst)
+                 {
+                   for(auto val : {point(0), point(1), 0.})
+                      vertices.append(val);
+                    int bcnr = self.GetSpline(i).bc;
+                    max_bcnr = max2(max_bcnr, bcnr);
+                    domains.append(bcnr);
+                    domains.append(self.GetSpline(i).leftdom);
+                    domains.append(self.GetSpline(i).rightdom);
+                 }
+
+               // segment data
+               auto pnt = self.splines[i]->GetPoint(0.5);
+               segment_points.append(py::make_tuple(pnt(0),pnt(1)));
+               auto normal = self.GetSpline(i).GetTangent(0.5);
+               std::swap(normal(0),normal(1));
+               normal(1) *= -1;
+               normal *= 1./sqrt(normal(0) * normal(0) + normal(1)*normal(1));
+               segment_normals.append(py::make_tuple(normal(0),normal(1)));
+               leftdom.append(self.GetSpline(i).leftdom);
+               rightdom.append(self.GetSpline(i).rightdom);
+             }
+           py::list bcnames;
+           for (int i = 1; i<max_bcnr + 1; i++)
+             bcnames.append(self.GetBCName(i));
+           segment_data["midpoints"] = segment_points;
+           segment_data["normals"] = segment_normals;
+           segment_data["leftdom"] = leftdom;
+           segment_data["rightdom"] = rightdom;
+           data["segment_data"] = segment_data;
+           data["vertices"] = vertices;
+           data["domains"] = domains;
+           data["min"] = min_val;
+           data["max"] = max_val;
+           data["bcnames"] = bcnames;
+           return data;
+         })
 	.def("PointData", FunctionPointer([](SplineGeometry2d &self)
 	  {
 		  py::list xpoints, ypoints, pointindex;
diff --git a/libsrc/occ/python_occ.cpp b/libsrc/occ/python_occ.cpp
index b565f7c0..8662ad10 100644
--- a/libsrc/occ/python_occ.cpp
+++ b/libsrc/occ/python_occ.cpp
@@ -30,6 +30,82 @@ DLL_HEADER void ExportNgOCC(py::module &m)
            self.HealGeometry();
            self.BuildFMap();
          },py::arg("tolerance")=1e-3, py::arg("fixsmalledges")=true, py::arg("fixspotstripfaces")=true, py::arg("sewfaces")=true, py::arg("makesolids")=true, py::arg("splitpartitions")=false,R"raw_string(Heal the OCCGeometry.)raw_string",py::call_guard<py::gil_scoped_release>())
+    .def("_visualizationData", [] (shared_ptr<OCCGeometry> occ_geo)
+         {
+           std::vector<float> vertices;
+           std::vector<int> trigs;
+           std::vector<float> normals;
+           std::vector<float> min = {std::numeric_limits<float>::max(),
+                               std::numeric_limits<float>::max(),
+                               std::numeric_limits<float>::max()};
+           std::vector<float> max = {std::numeric_limits<float>::lowest(),
+                               std::numeric_limits<float>::lowest(),
+                               std::numeric_limits<float>::lowest()};
+           std::vector<string> surfnames;
+           auto box = occ_geo->GetBoundingBox();
+           for(int i = 0; i < 3; i++)
+             {
+               min[i] = box.PMin()[i];
+               max[i] = box.PMax()[i];
+             }
+           occ_geo->BuildVisualizationMesh(0.01);
+           gp_Pnt2d uv;
+           gp_Pnt pnt;
+           gp_Vec n;
+           gp_Pnt p[3];
+           int count = 0;
+           for (int i = 1; i <= occ_geo->fmap.Extent(); i++)
+             {
+               surfnames.push_back("occ_surface" + to_string(i));
+               auto face = TopoDS::Face(occ_geo->fmap(i));
+               auto surf = BRep_Tool::Surface(face);
+               TopLoc_Location loc;
+               BRepAdaptor_Surface sf(face, Standard_False);
+               BRepLProp_SLProps prop(sf, 1, 1e-5);
+               Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc);
+               if (triangulation.IsNull())
+                 cout << "cannot visualize face " << i << endl;
+               trigs.reserve(trigs.size() + triangulation->NbTriangles()*4);
+               vertices.reserve(vertices.size() + triangulation->NbTriangles()*3*3);
+               normals.reserve(normals.size() + triangulation->NbTriangles()*3*3);
+               for (int j = 1; j < triangulation->NbTriangles()+1; j++)
+                 {
+                   auto triangle = (triangulation->Triangles())(j);
+                   for (int k = 1; k < 4; k++)
+                     p[k-1] = (triangulation->Nodes())(triangle(k)).Transformed(loc);
+                   for (int k = 1; k < 4; k++)
+                     {
+                       vertices.insert(vertices.end(),{float(p[k-1].X()), float(p[k-1].Y()), float(p[k-1].Z())});
+                       trigs.insert(trigs.end(),{count, count+1, count+2,i});
+                       count += 3;
+                       uv = (triangulation->UVNodes())(triangle(k));
+                       prop.SetParameters(uv.X(), uv.Y());
+                       if (prop.IsNormalDefined())
+                         n = prop.Normal();
+                       else
+                         {
+                           gp_Vec a(p[0], p[1]);
+                           gp_Vec b(p[0], p[2]);
+                           n = b^a;
+                         }
+                       if (face.Orientation() == TopAbs_REVERSED) n*= -1;
+                       normals.insert(normals.end(),{float(n.X()), float(n.Y()), float(n.Z())});
+                     }
+                 }
+             }
+            py::gil_scoped_acquire ac;
+            py::dict res;
+            py::list snames;
+            for(auto name : surfnames)
+              snames.append(py::cast(name));
+            res["vertices"] = MoveToNumpy(vertices);
+            res["triangles"] = MoveToNumpy(trigs);
+            res["normals"] = MoveToNumpy(normals);
+            res["surfnames"] = snames;
+            res["min"] = MoveToNumpy(min);
+            res["max"] = MoveToNumpy(max);
+            return res;
+         }, py::call_guard<py::gil_scoped_release>())
       ;
     m.def("LoadOCCGeometry",FunctionPointer([] (const string & filename)
            {
diff --git a/libsrc/stlgeom/python_stl.cpp b/libsrc/stlgeom/python_stl.cpp
index 8cef239d..66101d0f 100644
--- a/libsrc/stlgeom/python_stl.cpp
+++ b/libsrc/stlgeom/python_stl.cpp
@@ -20,6 +20,55 @@ DLL_HEADER void ExportSTL(py::module & m)
 {
   py::class_<STLGeometry,shared_ptr<STLGeometry>, NetgenGeometry> (m,"STLGeometry")
     .def(py::init<>())
+    .def("_visualizationData", [](shared_ptr<STLGeometry> stl_geo)
+         {
+           std::vector<float> vertices;
+           std::vector<int> trigs;
+           std::vector<float> normals;
+           std::vector<float> min = {std::numeric_limits<float>::max(),
+                               std::numeric_limits<float>::max(),
+                               std::numeric_limits<float>::max()};
+           std::vector<float> max = {std::numeric_limits<float>::lowest(),
+                               std::numeric_limits<float>::lowest(),
+                               std::numeric_limits<float>::lowest()};
+           std::vector<string> surfnames;
+
+           surfnames.push_back("stl");
+           vertices.reserve(stl_geo->GetNT()*3*3);
+           trigs.reserve(stl_geo->GetNT()*4);
+           normals.reserve(stl_geo->GetNT()*3*3);
+           size_t ii = 0;
+           for(int i = 0; i < stl_geo->GetNT(); i++)
+             {
+               auto& trig = stl_geo->GetTriangle(i+1);
+               for(int k = 0; k < 3; k++)
+                 {
+                   trigs.push_back(ii++);
+                   auto& pnt = stl_geo->GetPoint(trig[k]);
+                   for (int l = 0; l < 3; l++)
+                     {
+                       float val = pnt[l];
+                       vertices.push_back(val);
+                       min[l] = min2(min[l], val);
+                       max[l] = max2(max[l], val);
+                       normals.push_back(trig.Normal()[l]);
+                     }
+                 }
+               trigs.push_back(0);
+             }
+            py::gil_scoped_acquire ac;
+            py::dict res;
+            py::list snames;
+            for(auto name : surfnames)
+              snames.append(py::cast(name));
+            res["vertices"] = MoveToNumpy(vertices);
+            res["triangles"] = MoveToNumpy(trigs);
+            res["normals"] = MoveToNumpy(normals);
+            res["surfnames"] = snames;
+            res["min"] = MoveToNumpy(min);
+            res["max"] = MoveToNumpy(max);
+            return res;
+         }, py::call_guard<py::gil_scoped_release>())
     ;
   m.def("LoadSTLGeometry", FunctionPointer([] (const string & filename)
 					   {