From c9cd9eea2ced433021e4f5cbeda115313ee11cd9 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 12 Jul 2018 16:35:52 +0200 Subject: [PATCH] add functions to collect visualization data to python export of geometries --- libsrc/csg/python_csg.cpp | 62 +++++++++++++++++++++++++++ libsrc/general/ngpython.hpp | 10 +++++ libsrc/geom2d/python_geom2d.cpp | 73 +++++++++++++++++++++++++++++++ libsrc/occ/python_occ.cpp | 76 +++++++++++++++++++++++++++++++++ libsrc/stlgeom/python_stl.cpp | 49 +++++++++++++++++++++ 5 files changed, 270 insertions(+) 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 csg_geo) + { + std::vector vertices; + std::vector trigs; + std::vector normals; + std::vector min = {std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + std::vector max = {std::numeric_limits::lowest(), + std::numeric_limits::lowest(), + std::numeric_limits::lowest()}; + std::vector 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()) ; m.def("GenerateMesh", FunctionPointer diff --git a/libsrc/general/ngpython.hpp b/libsrc/general/ngpython.hpp index e0980d85..fa9862b1 100644 --- a/libsrc/general/ngpython.hpp +++ b/libsrc/general/ngpython.hpp @@ -2,11 +2,21 @@ #include #include +#include namespace py = pybind11; #include #include +template +py::array MoveToNumpy(std::vector& vec) +{ + auto newvec = new std::vector(); + std::swap(*newvec, vec); + auto capsule = py::capsule(newvec, [](void *v) { delete reinterpret_cast*>(v); }); + return py::array(newvec->size(), newvec->data(), capsule); +} + namespace PYBIND11_NAMESPACE { template bool CheckCast( py::handle obj ) { diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index 5ca29046..fa82c2ed 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -200,6 +200,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> 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()) + .def("_visualizationData", [] (shared_ptr occ_geo) + { + std::vector vertices; + std::vector trigs; + std::vector normals; + std::vector min = {std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + std::vector max = {std::numeric_limits::lowest(), + std::numeric_limits::lowest(), + std::numeric_limits::lowest()}; + std::vector 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()) ; 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_, NetgenGeometry> (m,"STLGeometry") .def(py::init<>()) + .def("_visualizationData", [](shared_ptr stl_geo) + { + std::vector vertices; + std::vector trigs; + std::vector normals; + std::vector min = {std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + std::vector max = {std::numeric_limits::lowest(), + std::numeric_limits::lowest(), + std::numeric_limits::lowest()}; + std::vector 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()) ; m.def("LoadSTLGeometry", FunctionPointer([] (const string & filename) {