From 356e78b809c8536d153644e23c2c60a06dc7487d Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 25 Jul 2023 12:11:00 +0200 Subject: [PATCH 1/3] Fix Point3d Python operators --- libsrc/meshing/python_mesh.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 7f103936..0ef1332a 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -183,6 +183,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(py::self+Vec<2>()) .def(py::self-Vec<2>()) .def("__getitem__", [](Point<2>& self, int index) { return self[index]; }) + .def("__len__", [](Point<2>& /*unused*/) { return 2; }) ; py::implicitly_convertible>(); @@ -198,7 +199,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(py::self-py::self) .def(py::self+Vec<3>()) .def(py::self-Vec<3>()) - .def("__getitem__", [](Point<2>& self, int index) { return self[index]; }) + .def("__getitem__", [](Point<3>& self, int index) { return self[index]; }) + .def("__len__", [](Point<3>& /*unused*/) { return 3; }) ; py::implicitly_convertible>(); From 9ae05ab71260c392981c172e89f8caad7e72bb6d Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 25 Jul 2023 12:11:13 +0200 Subject: [PATCH 2/3] add mesh.bounding_box in Python --- libsrc/meshing/python_mesh.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 0ef1332a..eea4c3be 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -753,6 +753,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def_property_readonly("_timestamp", &Mesh::GetTimeStamp) .def_property_readonly("ne", [](Mesh& m) { return m.GetNE(); }) + .def_property_readonly("bounding_box", [](Mesh& m) { + Point3d pmin, pmax; + m.GetBox(pmin, pmax); + return py::make_tuple( Point<3>(pmin),Point<3>(pmax)); + }) .def("Partition", [](shared_ptr self, int numproc) { self->ParallelMetis(numproc); }, py::arg("numproc")) From 11a898442826753fd206efaa5c01b09419cd615a Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 25 Jul 2023 12:30:24 +0200 Subject: [PATCH 3/3] Webgui Draw for netgen.mesh --- python/webgui.py | 183 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 182 insertions(+), 1 deletion(-) diff --git a/python/webgui.py b/python/webgui.py index 7131e234..d357ac2c 100644 --- a/python/webgui.py +++ b/python/webgui.py @@ -9,6 +9,8 @@ import webgui_jupyter_widgets.widget as wg from packaging.version import parse +import netgen.meshing as ng + if parse(webgui_jupyter_widgets.__version__) >= parse("0.2.18"): _default_width = None _default_height = None @@ -28,6 +30,176 @@ def register_draw_type(*types): return inner +_bernstein_cache = {} + + +def GetIBernsteinBasis(etype, order): + if (etype, order) in _bernstein_cache: + return _bernstein_cache[(etype, order)] + bvals = None + + if etype == "segment": + + 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 = np.zeros( + (order + 1, order + 1), dtype=float + ) # .Matrix(order+1,order+1) + for i in range(order + 1): + for j in range(order + 1): + bvals[i, j] = Bernstein(i / order, j, order) + + if etype == "trig": + + 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) + ) + + og = order + ndtrig = int((og + 1) * (og + 2) / 2) + bvals = np.zeros((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 + + if bvals is None: + raise RuntimeError(f"Unkown element type {etype}") + + ibvals = _bernstein_cache[(etype, order)] = np.linalg.inv(bvals) + return ibvals + + +def GetWireframePoints(etype, order): + n = order + if etype == "trig": + return np.array( + [(i / n, 0) for i in range(n + 1)] + + [(0, i / n) for i in range(n + 1)] + + [(i / n, 1.0 - i / n) for i in range(n + 1)] + ) + if etype == "quad": + return np.array( + [(i / n, 0) for i in range(n + 1)] + + [(0, i / n) for i in range(n + 1)] + + [(i / n, 1.0) for i in range(n + 1)] + + [(1.0, i / n) for i in range(n + 1)] + ) + + raise RuntimeError(f"Unknown element type {etype}") + + +def GetElementPoints(etype, order): + n = order + if etype == "trig": + return np.array( + [(i / n, j / n) for j in range(n + 1) for i in range(n + 1 - j)] + ) + if etype == "quad": + return np.array( + [(i / n, j / n) for j in range(n + 1) for i in range(n + 1 - j)] + + [(1 - i / n, 1 - j / n) for j in range(n + 1) for i in range(n + 1 - j)] + ) + + raise RuntimeError(f"Unknown element type {etype}") + + +def MapBernstein(pnts, etype, order): + """ + Maps function values at equidistant control points to the Bernstein basis function. + Parameters: + pnts (numpy.ndarray): The input control points with shape (number_of_elements, points_per_element, function_dimension) + point_per_element must be a multiple of the basis size + etype (str): Element type (currently ignored and trig assumed) + order (int): Polynomial order + + Returns: + numpy.ndarray: The mapped points with the shape (points_per_element, number_of_elements, function_dimension) + """ + ibvals = GetIBernsteinBasis(etype, order) + # for wireframe or subdivided elements, we have multiple point sets per element + # so do a reshape to simulate more elements with correct number of control points per element instead + if pnts.shape[1] != ibvals.shape[0]: + pnts = pnts.reshape((-1, ibvals.shape[0], pnts.shape[2])) + + points = np.zeros(pnts.shape, dtype=np.float32).transpose(1, 0, 2) + for i in range(points.shape[2]): + points[:, :, i] = np.tensordot(ibvals, pnts[:, :, i], axes=(1, 1)) + return points + + +@register_draw_type(ng.Mesh) +def GetData(mesh, args, kwargs): + d = {} + d["gui_settings"] = kwargs["settings"] + d["mesh_dim"] = mesh.dim + + pmin, pmax = mesh.bounding_box + diag = pmax - pmin + pmid = pmin + 0.5 * diag + d["mesh_center"] = [pmid[i] for i in range(3)] + d["mesh_radius"] = diag.Norm() + + d["funcdim"] = 0 + d["show_mesh"] = True + d["draw_surf"] = True + d["funcmin"] = 0.0 + d["funcmax"] = 1.0 + + # Generate surface element data + # webgui code assumes 4 scalar fields (x,y,z, mesh_index) + # TODO: other element types than trigs + order = kwargs["order"] + refpts = GetElementPoints("trig", order) + pnts = np.ndarray((len(mesh.Elements2D()), refpts.shape[0], 4)) + mesh.CalcElementMapping(refpts, pnts) + + # set mesh_index + for i, el in enumerate(mesh.Elements2D()): + pnts[i, :, 3] = el.index - 1 + fds = mesh.FaceDescriptors() + d["colors"] = [fd.color for fd in fds] + d["mesh_regions_2d"] = len(fds) + d["names"] = [fd.bcname for fd in fds] + + d["Bezier_trig_points"] = MapBernstein(pnts, "trig", order) + d["order2d"] = order + + # Generate wireframe data + refpts = GetWireframePoints("trig", order) + pnts = np.ndarray((len(mesh.Elements2D()), refpts.shape[0], 4)) + mesh.CalcElementMapping(refpts, pnts) + d["Bezier_points"] = MapBernstein(pnts, "segment", order) + d["show_wireframe"] = True + + # TODO: Generate edge data + d["edges"] = [] + + # encode data as b64 + for name in ["Bezier_trig_points", "edges", "Bezier_points"]: + pnew = [] + for plist in d[name]: + pnew.append(encodeData(np.array(plist, dtype=np.float32))) + d[name] = pnew + return d + + class WebGLScene(BaseWebGuiScene): def __init__(self, obj, args=[], kwargs={}): self.obj = obj @@ -91,7 +263,7 @@ class WebGLScene(BaseWebGuiScene): d["funcmin"] = kwargs["min"] if "max" in kwargs: d["funcmax"] = kwargs["max"] - d['autoscale'] = kwargs['autoscale'] + d["autoscale"] = kwargs["autoscale"] if "vectors" in kwargs: d["vectors"] = True @@ -122,6 +294,14 @@ class WebGLScene(BaseWebGuiScene): d["fullscreen"] = kwargs["fullscreen"] if "gui_settings" not in d: d["gui_settings"] = self.kwargs["settings"] + + d["objects"] = [] + for obj in kwargs["objects"]: + if isinstance(obj, dict): + d["objects"].append(obj) + else: + d["objects"].append(obj._GetWebguiData()) + return d @@ -146,6 +326,7 @@ bezier_trig_trafos = {} # cache trafos for different orders # scene.GenerateHTML(filename=filename) # return scene + def _get_draw_default_args(): return dict( name="function",