From ae6b23fffc149551b8e5a343a0ac4aa1e1ff8f4e Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Thu, 12 Aug 2021 11:20:07 +0200 Subject: [PATCH] add functionality to draw occ shapes with webgui --- libsrc/occ/python_occ_shapes.cpp | 97 ++++++++++++ python/CMakeLists.txt | 1 + python/webgui.py | 247 +++++++++++++++++++++++++++++++ 3 files changed, 345 insertions(+) create mode 100644 python/webgui.py diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index abd08a04..ea661ccb 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -56,6 +56,35 @@ using namespace netgen; +void ExtractFaceData( const TopoDS_Face & face, int index, std::vector * p, Box<3> & box ) +{ + Handle(Geom_Surface) surf = BRep_Tool::Surface (face); + + TopLoc_Location loc; + Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation (face, loc); + + bool flip = TopAbs_REVERSED == face.Orientation(); + + int ntriangles = triangulation -> NbTriangles(); + for (int j = 1; j <= ntriangles; j++) + { + Poly_Triangle triangle = (triangulation -> Triangles())(j); + std::array,3> pts; + for (int k = 0; k < 3; k++) + pts[k] = occ2ng( (triangulation -> Nodes())(triangle(k+1)).Transformed(loc) ); + + if(flip) + Swap(pts[1], pts[2]); + + for (int k = 0; k < 3; k++) + { + box.Add(pts[k]); + for (int d = 0; d < 3; d++) + p[k].push_back( pts[k][d] ); + p[k].push_back( index ); + } + } +} py::object CastShape(const TopoDS_Shape & s) { @@ -543,6 +572,74 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) // return MoveToNumpyArray(triangles); return triangles; }) + .def("_webgui_data", [](const TopoDS_Shape & shape) + { + BRepTools::Clean (shape); + double deflection = 0.01; + BRepMesh_IncrementalMesh (shape, deflection, true); + // triangulation = BRep_Tool::Triangulation (face, loc); + + std::vector p[3]; + py::list names, colors; + + int index = 0; + + Box<3> box(Box<3>::EMPTY_BOX); + for (TopExp_Explorer e(shape, TopAbs_FACE); e.More(); e.Next()) + { + TopoDS_Face face = TopoDS::Face(e.Current()); + // Handle(TopoDS_Face) face = e.Current(); + ExtractFaceData(face, index, p, box); + auto & props = OCCGeometry::global_shape_properties[face.TShape()]; + if(props.col) + { + auto & c = *props.col; + colors.append(py::make_tuple(c[0], c[1], c[2])); + } + else + colors.append(py::make_tuple(0.0, 1.0, 0.0)); + if(props.name) + { + names.append(*props.name); + } + else + names.append(""); + index++; + } + + auto center = box.Center(); + + py::list mesh_center; + mesh_center.append(center[0]); + mesh_center.append(center[1]); + mesh_center.append(center[2]); + py::dict data; + data["ngsolve_version"] = "Netgen x.x"; // TODO + data["mesh_dim"] = 3; // TODO + data["mesh_center"] = mesh_center; + data["mesh_radius"] = box.Diam()/2; + data["order2d"] = 1; + data["order3d"] = 0; + data["draw_vol"] = false; + data["draw_surf"] = true; + data["funcdim"] = 0; + data["show_wireframe"] = false; + data["show_mesh"] = true; + data["edges"] = py::list{}; + data["Bezier_points"] = py::list{}; + py::list points; + points.append(p[0]); + points.append(p[1]); + points.append(p[2]); + data["Bezier_trig_points"] = points; + data["funcmin"] = 0; + data["funcmax"] = 1; + data["mesh_regions_2d"] = index; + data["autoscale"] = false; + data["colors"] = colors; + data["names"] = names; + return data; + }) ; py::class_ (m, "TopoDS_Vertex") diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 47910190..c3dc982b 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -3,6 +3,7 @@ configure_file(__init__.py ${CMAKE_CURRENT_BINARY_DIR}/__init__.py @ONLY) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/__init__.py meshing.py csg.py geom2d.py stl.py gui.py NgOCC.py occ.py read_gmsh.py + webgui.py DESTINATION ${NG_INSTALL_DIR_PYTHON}/${NG_INSTALL_SUFFIX} COMPONENT netgen ) diff --git a/python/webgui.py b/python/webgui.py new file mode 100644 index 00000000..48f65811 --- /dev/null +++ b/python/webgui.py @@ -0,0 +1,247 @@ +import math +import numpy as np +from time import time +import os + +from webgui_jupyter_widgets import BaseWebGuiScene, encodeData, WebGuiDocuWidget +import webgui_jupyter_widgets.widget as wg + +class WebGLScene(BaseWebGuiScene): + def __init__(self, mesh, clipping, on_init): + from IPython.display import display, Javascript + import threading + self.mesh = mesh + self.clipping = clipping + self.on_init = on_init + + def GetData(self, set_minmax=True): + import json + # d = BuildRenderData(self.mesh, self.cf, self.order, draw_surf=self.draw_surf, draw_vol=self.draw_vol, deformation=self.deformation, region=self.region) + d = self.mesh._webgui_data() + bp = d['Bezier_trig_points'] + for i in range(len(bp)): + bp[i] = encodeData(np.array(bp[i])) + + + if self.clipping is not None: + d['clipping'] = True + if isinstance(self.clipping, dict): + allowed_args = ("x", "y", "z", "dist", "function", "pnt", "vec") + if "vec" in self.clipping: + vec = self.clipping["vec"] + self.clipping["x"] = vec[0] + self.clipping["y"] = vec[1] + self.clipping["z"] = vec[2] + if "pnt" in self.clipping: + d['mesh_center'] = list(self.clipping["pnt"]) + for name, val in self.clipping.items(): + if not (name in allowed_args): + raise Exception('Only {} allowed as arguments for clipping!'.format(", ".join(allowed_args))) + d['clipping_' + name] = val + + if self.on_init: + d['on_init'] = self.on_init + + return d + + +bezier_trig_trafos = { } # cache trafos for different orders + +def BuildRenderData(mesh, func, order=2, draw_surf=True, draw_vol=True, deformation=None, region=True): + d = {} + d['ngsolve_version'] = "Netgen" + d['mesh_dim'] = 3 # mesh.dim TODO + + d['order2d'] = 1 + d['order3d'] = 0 + + d['draw_vol'] = False + d['draw_surf'] = True + d['funcdim'] = 1 + + func2 = None + if func and func.is_complex: + d['is_complex'] = True + func1 = func[0].real + func2 = ngs.CoefficientFunction( (func[0].imag, 0.0) ) + elif func and func.dim>1: + func1 = func[0] + func2 = ngs.CoefficientFunction( tuple(func[i] if i we are just drawing a mesh, eval mesh element index instead + mats = mesh.GetMaterials() + bnds = mesh.GetBoundaries() + nmats = len(mesh.GetMaterials()) + nbnds = len(mesh.GetBoundaries()) + n = max(nmats, nbnds) + func1 = ngs.CoefficientFunction(list(range(n))) + n_regions = [0, 0, nmats, nbnds] + d['mesh_regions_2d'] = n_regions[mesh.dim] + d['mesh_regions_3d'] = nmats if mesh.dim==3 else 0 + d['funcdim'] = 0 + func1 = ngs.CoefficientFunction( (ngs.x, ngs.y, ngs.z, func1 ) ) + func0 = ngs.CoefficientFunction( (ngs.x, ngs.y, ngs.z, 0.0 ) ) + if deformation is not None: + func1 += ngs.CoefficientFunction((deformation, 0.0)) + func0 += ngs.CoefficientFunction((deformation, 0.0)) + + d['show_wireframe'] = False + d['show_mesh'] = True + if order2d>0: + og = order2d + d['show_wireframe'] = True + d['show_mesh'] = True + timer2.Start() + + timer3Bvals.Start() + + # transform point-values to Bernsteinbasis + def Binomial(n,i): return math.factorial(n) / math.factorial(i) / math.factorial(n-i) + def Bernstein(x, i, n): return Binomial(n,i) * x**i*(1-x)**(n-i) + Bvals = ngs.Matrix(og+1,og+1) + for i in range(og+1): + for j in range(og+1): + Bvals[i,j] = Bernstein(i/og, j, og) + iBvals = Bvals.I + timer3Bvals.Stop() + # print (Bvals) + # print (iBvals) + + + Bezier_points = [] + + ipts = [(i/og,0) for i in range(og+1)] + [(0, i/og) for i in range(og+1)] + [(i/og,1.0-i/og) for i in range(og+1)] + ir_trig = ngs.IntegrationRule(ipts, [0,]*len(ipts)) + ipts = [(i/og,0) for i in range(og+1)] + [(0, i/og) for i in range(og+1)] + [(i/og,1.0) for i in range(og+1)] + [(1.0, i/og) for i in range(og+1)] + ir_quad = ngs.IntegrationRule(ipts, [0,]*len(ipts)) + + vb = [ngs.VOL, ngs.BND][mesh.dim-2] + if region and region.VB() == vb: + vb = region + cf = func1 if draw_surf else func0 + timer2map.Start() + pts = mesh.MapToAllElements({ngs.ET.TRIG: ir_trig, ngs.ET.QUAD: ir_quad}, vb) + timer2map.Stop() + pmat = cf(pts) + + timermult.Start() + pmat = pmat.reshape(-1, og+1, 4) + if False: + BezierPnts = np.tensordot(iBvals.NumPy(), pmat, axes=(1,1)) + else: + BezierPnts = np.zeros( (og+1, pmat.shape[0], 4) ) + for i in range(4): + ngsmat = ngs.Matrix(pmat[:,:,i].transpose()) + BezierPnts[:,:,i] = iBvals * ngsmat + timermult.Stop() + + timer2list.Start() + for i in range(og+1): + Bezier_points.append(encodeData(BezierPnts[i])) + timer2list.Stop() + + d['Bezier_points'] = Bezier_points + + ipts = [(i/og,0) for i in range(og+1)] + ir_seg = ngs.IntegrationRule(ipts, [0,]*len(ipts)) + vb = [ngs.VOL, ngs.BND, ngs.BBND][mesh.dim-1] + if region and region.VB() == vb: + vb = region + pts = mesh.MapToAllElements(ir_seg, vb) + pmat = func0(pts) + pmat = pmat.reshape(-1, og+1, 4) + edge_data = np.tensordot(iBvals.NumPy(), pmat, axes=(1,1)) + edges = [] + for i in range(og+1): + edges.append(encodeData(edge_data[i])) + d['edges'] = edges + + ndtrig = int((og+1)*(og+2)/2) + + if og in bezier_trig_trafos.keys(): + iBvals_trig = bezier_trig_trafos[og] + else: + def BernsteinTrig(x, y, i, j, n): + return math.factorial(n)/math.factorial(i)/math.factorial(j)/math.factorial(n-i-j) \ + * x**i*y**j*(1-x-y)**(n-i-j) + Bvals = ngs.Matrix(ndtrig, ndtrig) + ii = 0 + for ix in range(og+1): + for iy in range(og+1-ix): + jj = 0 + for jx in range(og+1): + for jy in range(og+1-jx): + Bvals[ii,jj] = BernsteinTrig(ix/og, iy/og, jx, jy, og) + jj += 1 + ii += 1 + iBvals_trig = Bvals.I + bezier_trig_trafos[og] = iBvals_trig + + + # Bezier_points = [ [] for i in range(ndtrig) ] + Bezier_points = [] + + ipts = [(i/og,j/og) for j in range(og+1) for i in range(og+1-j)] + ir_trig = ngs.IntegrationRule(ipts, [0,]*len(ipts)) + ipts = ([(i/og,j/og) for j in range(og+1) for i in range(og+1-j)] + + [(1-i/og,1-j/og) for j in range(og+1) for i in range(og+1-j)]) + ir_quad = ngs.IntegrationRule(ipts, [0,]*len(ipts)) + + vb = [ngs.VOL, ngs.BND][mesh.dim-2] + if region and region.VB() == vb: + vb = region + pts = mesh.MapToAllElements({ngs.ET.TRIG: ir_trig, ngs.ET.QUAD: ir_quad}, vb) + + pmat = ngs.CoefficientFunction( func1 if draw_surf else func0 ) (pts) + + + funcmin = np.min(pmat[:,3]) + funcmax = np.max(pmat[:,3]) + pmin = np.min(pmat[:,0:3], axis=0) + pmax = np.max(pmat[:,0:3], axis=0) + + mesh_center = 0.5*(pmin+pmax) + mesh_radius = np.linalg.norm(pmax-pmin)/2 + + pmat = pmat.reshape(-1, len(ir_trig), 4) + + BezierPnts = np.tensordot(iBvals_trig.NumPy(), pmat, axes=(1,1)) + + for i in range(ndtrig): + Bezier_points.append(encodeData(BezierPnts[i])) + + + d['Bezier_trig_points'] = Bezier_points + d['mesh_center'] = list(mesh_center) + d['mesh_radius'] = mesh_radius + + + + if func: + d['funcmin'] = funcmin + d['funcmax'] = funcmax + return d + +def Draw(shape, clipping=None, js_code=None, filename=""): + # todo: also handle occ geometry, list of shapes, etc. + + scene = WebGLScene(shape, clipping=clipping, on_init=js_code) + + if wg._IN_IPYTHON: + if wg._IN_GOOGLE_COLAB: + from IPython.display import display, HTML + html = scene.GenerateHTML() + display(HTML(html)) + else: + # render scene using widgets.DOMWidget + scene.Draw() + return scene + else: + if filename: + scene.GenerateHTML(filename=filename) + return scene +