diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp index 52e7a20e..0d4bcb00 100644 --- a/libsrc/geom2d/geometry2d.hpp +++ b/libsrc/geom2d/geometry2d.hpp @@ -175,6 +175,7 @@ namespace netgen void CopyEdgeMesh (int from, int to, Mesh & mesh2d, Point3dTree & searchtree); + size_t GetNDomains() const { return materials.Size(); } void GetMaterial (int domnr, char* & material ); void SetMaterial (int domnr, const string & material); diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index a16c13f7..f2b20127 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -15,6 +15,22 @@ namespace netgen DLL_HEADER void ExportGeom2d(py::module &m) { + py::class_> + (m, "Spline", "Spline of a SplineGeometry object") + .def_property("leftdom", [] (SplineSegExt& self) { return self.leftdom; }, + [](SplineSegExt& self, int dom) { self.leftdom = dom; }) + .def_property("rightdom", [] (SplineSegExt& self) { return self.rightdom; }, + [](SplineSegExt& self, int dom) { self.rightdom = dom; }) + .def_property_readonly("bc", [] (SplineSegExt& self) { return self.bc; }) + .def("GetNormal", [](SplineSegExt& self, double t) + { + auto tang = self.GetTangent(t).Normalize(); + return Vec<2>(tang[1], -tang[0]); + }) + .def("StartPoint", [](SplineSegExt& self) { return Point<2>(self.StartPI()); }) + .def("EndPoint", [](SplineSegExt& self) { return Point<2>(self.EndPI()); }) + ; + py::class_> (m, "SplineGeometry", "a 2d boundary representation geometry model by lines and splines", @@ -121,38 +137,21 @@ DLL_HEADER void ExportGeom2d(py::module &m) seg->copyfrom = -1; self.AppendSegment(seg); }), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = py::int_(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)); - // SplineSegExt * seg = new SplineSegExt(*l); - // seg->leftdom = 1;// leftdomain; - // seg->rightdom = 0;// rightdomain; - // seg->hmax = 1e99; - // seg->reffak = 1; - // seg->copyfrom = -1; - - // self.AppendSegment(seg); - // }))//, (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)); - // SplineSegExt * seg = new SplineSegExt(*seg3); - // seg->leftdom = 1;// leftdomain; - // seg->rightdom = 0;// rightdomain; - // seg->hmax = 1e99; - // seg->reffak = 1; - // seg->copyfrom = -1; - // self.AppendSegment(seg); - // }))//, (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) .def("SetDomainMaxH", &SplineGeometry2d::SetDomainMaxh) + .def("GetBCName", [](SplineGeometry2d& self, size_t index) { return self.GetBCName(index); }) + .def("GetNDomains", [](SplineGeometry2d& self) { return self.GetNDomains(); }) + .def("GetNSplines", [](SplineGeometry2d& self) { return self.splines.Size(); }) + .def("GetSpline", [](SplineGeometry2d& self, size_t index) + { return shared_ptr(&self.GetSpline(index), NOOP_Deleter); }, + py::return_value_policy::reference_internal) + .def("GetNPoints", [](SplineGeometry2d& self) { return self.GetNP(); }) + .def("GetPoint", [](SplineGeometry2d& self, size_t index) { return Point<2>(self.GetPoint(index)); }) - .def("PlotData", FunctionPointer([](SplineGeometry2d &self) { Box<2> box(self.GetBoundingBox()); diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index b0c2a5c7..436835ef 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -80,6 +80,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(py::self-py::self) .def(py::self+Vec<2>()) .def(py::self-Vec<2>()) + .def("__getitem__", [](Point<2>& self, int index) { return self[index]; }) ; py::class_> (m, "Point3d") @@ -88,6 +89,7 @@ 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]; }) ; m.def ("Pnt", FunctionPointer diff --git a/python/geom2d.py b/python/geom2d.py index 88856964..34b37870 100644 --- a/python/geom2d.py +++ b/python/geom2d.py @@ -39,11 +39,107 @@ def MakeCircle (geo, c, r, **args): geo.Append( ["spline3", pts[p1], pts[p2], pts[p3]], **args) + +def CreatePML(geo, pml_size, tol=1e-12): + """Create a pml layer around the geometry. This function works only on convex geometries and +the highest existing domain number must be named by using the function geo.SetMaterial(domnr, name). +Points in the geometry are assumed to be the same if (pt1 - pt2).Norm() < tol. +Returned is a dict with information to create the pml layer: + normals: A dict from the names of the linear pml domains to the normal vectors pointing inside the pml.""" + + def Start(spline): + if spline.rightdom == 0: + return spline.StartPoint() + return spline.EndPoint() + def End(spline): + if spline.rightdom == 0: + return spline.EndPoint() + return spline.StartPoint() + + splines = [] + for i in range(geo.GetNSplines()): + splines.append(geo.GetSpline(i)) + border = [] + is_closed = False + current_endpoint = None + while not is_closed: + for spline in splines: + if spline.leftdom == 0 or spline.rightdom == 0: + if current_endpoint is not None: + if (Start(spline)-current_endpoint).Norm() < tol: + border.append(spline) + current_endpoint = End(spline) + if (current_endpoint - startpoint).Norm() < tol: + is_closed = True + break + else: + startpoint = Start(spline) + current_endpoint = End(spline) + border.append(spline) + break + else: + raise Exception("Couldn't find closed spline around domain") + endpointindex_map = [] + for spline in border: + pnt = End(spline) + for i in range(geo.GetNPoints()): + if (pnt - geo.GetPoint(i)).Norm() < tol: + endpointindex_map.append(i) + break + else: + raise Exception("Couldn't find endpoint of spline in geometry") + start_ndoms = ndoms = geo.GetNDomains() + 1 + new_spline_domains = [] + normals = {} + for i, spline in enumerate(border): + if i == 0: + global_start = Start(spline) + pml_size * spline.GetNormal(0) + global_start_pnt = current_start = geo.AppendPoint(global_start[0], global_start[1]) + next_spline = border[(i+1)%len(border)] + new_end = End(spline) + pml_size * spline.GetNormal(1) + spline_name = geo.GetBCName(spline.bc) + if (new_end - global_start).Norm() < tol: + new_spline_domains.append(ndoms) + geo.Append(["line", current_start, global_start_pnt], bc="outer_" + spline_name, leftdomain = ndoms) + geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms) + geo.SetMaterial(ndoms, "pml_" + spline_name) + normals["pml_" + spline_name] = spline.GetNormal(0) + ndoms += 1 + break + end = geo.AppendPoint(new_end[0], new_end[1]) + new_spline_domains.append(ndoms) + geo.Append(["line", current_start, end], bc="outer_" + spline_name, leftdomain = ndoms) + geo.Append(["line", end, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1) + geo.SetMaterial(ndoms, "pml_" + spline_name) + normals["pml_" + spline_name] = spline.GetNormal(0) + ndoms += 1 + new_start = Start(next_spline) + pml_size * next_spline.GetNormal(0) + if (new_start - global_start).Norm() < tol: + geo.Append(["line", end, global_start_pnt], bc="outer", leftdomain = ndoms) + geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms) + geo.SetMaterial(ndoms, "pml_corner") + ndoms += 1 + break + if (new_end - new_start).Norm() < tol: + current_start = end + else: + current_start = geo.AppendPoint(new_start[0], new_start[1]) + geo.Append(["line", end, current_start], bc="outer", leftdomain = ndoms) + geo.Append(["line", current_start, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1) + geo.SetMaterial(ndoms, "pml_corner") + ndoms += 1 + for spline, domnr in zip(border, new_spline_domains): + if spline.leftdom == 0: + spline.leftdom = domnr + else: + spline.rightdom = domnr + return {"normals" : normals} + SplineGeometry.AddCircle = lambda geo, c, r, **args : MakeCircle(geo, c, r, **args) SplineGeometry.AddRectangle = lambda geo, p1, p2, **args : MakeRectangle(geo, p1, p2, **args) SplineGeometry.AddSegment = lambda *args, **kwargs : SplineGeometry.Append(*args, **kwargs) SplineGeometry.AddPoint = lambda *args, **kwargs : SplineGeometry.AppendPoint(*args, **kwargs) - +SplineGeometry.CreatePML = CreatePML __all__ = ['SplineGeometry', 'unit_square']