From 666fb2ee86b2b7f40b9553035642c8a3c5c3162e Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Aug 2022 11:11:15 +0200 Subject: [PATCH] fix boundarylayer 2d code (now single line segments, not per face) --- libsrc/meshing/boundarylayer.cpp | 178 +++++++++++-------------------- libsrc/meshing/fieldlines.cpp | 2 +- libsrc/meshing/python_mesh.cpp | 1 + libsrc/occ/occgeom.hpp | 6 ++ libsrc/occ/python_occ.cpp | 6 +- 5 files changed, 75 insertions(+), 118 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 7da8a4f2..db6efc08 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -436,7 +436,6 @@ namespace netgen } } - // seg.edgenr = topo.GetEdge(segi)+1; segments.Append(seg); } } @@ -1394,20 +1393,24 @@ namespace netgen compress = PointIndex::INVALID; PointIndex cnt = PointIndex::BASE; - for(const auto & seg : mesh.LineSegments()) - if (seg.si == domain) - for (auto pi : {seg[0], seg[1]}) - if (compress[pi]==PointIndex{PointIndex::INVALID}) - { - meshing.AddPoint(mesh[pi], pi); - compress[pi] = cnt++; - } + auto p2sel = mesh.CreatePoint2SurfaceElementTable(); + Array domain_segments; PointGeomInfo gi; gi.trignum = domain; - for(const auto & seg : mesh.LineSegments()) - if (seg.si == domain) - meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], gi, gi); + for(auto seg : mesh.LineSegments()) + { + for (auto pi : {seg[0], seg[1]}) + if (compress[pi]==PointIndex{PointIndex::INVALID}) + { + meshing.AddPoint(mesh[pi], pi); + compress[pi] = cnt++; + } + if(seg.domin == domain) + meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], gi, gi); + if(seg.domout == domain) + meshing.AddBoundaryElement (compress[seg[1]], compress[seg[0]], gi, gi); + } auto oldnf = mesh.GetNSE(); auto res = meshing.GenerateMesh (mesh, mp, mp.maxh, domain); @@ -1450,17 +1453,20 @@ namespace netgen } mesh.Compress(); + mesh.CalcSurfacesOfNode(); mesh.OrderElements(); mesh.SetNextMajorTimeStamp(); - } int GenerateBoundaryLayer2 (Mesh & mesh, int domain, const Array & thicknesses, bool should_make_new_domain, const Array & boundaries) { + mesh.GetTopology().SetBuildVertex2Element(true); + mesh.UpdateTopology(); + const auto & line_segments = mesh.LineSegments(); SegmentIndex first_new_seg = mesh.LineSegments().Range().Next(); int np = mesh.GetNP(); - int nseg = mesh.GetNSeg(); + int nseg = line_segments.Size(); int ne = mesh.GetNSE(); mesh.UpdateTopology(); @@ -1482,23 +1488,10 @@ namespace netgen auto & meshtopo = mesh.GetTopology(); - Array seg2surfel(mesh.GetNSeg()); - seg2surfel = -1; - NgArray temp_els; - for(auto si : Range(mesh.LineSegments())) - { - meshtopo.GetSegmentSurfaceElements ( si+1, temp_els ); - // NgArray surfeledges; - // meshtopo.GetSurfaceElementEdges(si+1, surfeledges); - for(auto seli : temp_els) - if(mesh[seli].GetIndex() == mesh[si].si) - seg2surfel[si] = seli; - } - Array segments; // surface index map - Array si_map(mesh.GetNFD()+1); + Array si_map(mesh.GetNFD()+2); si_map = -1; int fd_old = mesh.GetNFD(); @@ -1506,12 +1499,14 @@ namespace netgen int max_edge_nr = -1; int max_domain = -1; - for(const auto& seg : mesh.LineSegments()) + for(const auto& seg : line_segments) { if(seg.epgeominfo[0].edgenr > max_edge_nr) max_edge_nr = seg.epgeominfo[0].edgenr; - if(seg.si > max_domain) - max_domain = seg.si; + if(seg.domin > max_domain) + max_domain = seg.domin; + if(seg.domout > max_domain) + max_domain = seg.domout; } int new_domain = max_domain+1; @@ -1527,95 +1522,43 @@ namespace netgen for(auto edgenr : boundaries) active_boundaries.SetBit(edgenr); - for(auto segi : Range(mesh.LineSegments())) + for(auto segi : Range(line_segments)) { - const auto seg = mesh[segi]; - if(active_boundaries.Test(seg.epgeominfo[0].edgenr) && seg.si==domain) + const auto seg = line_segments[segi]; + if(active_boundaries.Test(seg.epgeominfo[0].edgenr) && (seg.domin==domain || seg.domout==domain)) active_segments.SetBit(segi); } - for(auto segi : Range(mesh.LineSegments())) { - const auto& seg = mesh[segi]; - auto si = seg.si; - - if(si_map[si]!=-1) - continue; - - if(!active_segments.Test(segi)) - continue; - FaceDescriptor new_fd(0, 0, 0, -1); new_fd.SetBCProperty(new_domain); int new_fd_index = mesh.AddFaceDescriptor(new_fd); - si_map[si] = new_domain; if(should_make_new_domain) - mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(si-1)); + mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(domain-1)); } - for(auto si : Range(mesh.LineSegments())) + for(auto segi : Range(line_segments)) { - if(segs_done[si]) continue; - segs_done.SetBit(si); - const auto& segi = mesh[si]; - if(si_map[segi.si] == -1) continue; - if(!active_boundaries.Test(segi.epgeominfo[0].edgenr)) + if(segs_done[segi]) continue; + segs_done.SetBit(segi); + const auto& seg = line_segments[segi]; + if(seg.domin != domain && seg.domout != domain) continue; + if(!active_boundaries.Test(seg.epgeominfo[0].edgenr)) continue; - moved_segs.Append(si); + moved_segs.Append(segi); } // calculate growth vectors (average normal vectors of adjacent segments at each point) for (auto si : moved_segs) { - auto & seg = mesh[si]; - - meshtopo.GetSegmentSurfaceElements ( si+1, temp_els ); - ArrayMem seg_domains; - - temp_els.SetSize(0); - if(seg2surfel[si]!=-1) - temp_els.Append(seg2surfel[si]); - - int n_temp_els = temp_els.Size(); - if(n_temp_els==0) - continue; - - int dom0 = mesh[temp_els[0]].GetIndex(); - int dom1 = n_temp_els==2 ? mesh[temp_els[1]].GetIndex() : 0; - - bool in_dom0 = dom0 == domain; - bool in_dom1 = dom1 == domain; - - if(!in_dom0 && !in_dom1) - continue; - - int side = in_dom0 ? 0 : 1; - - auto & sel = mesh[ temp_els[side] ]; - - int domain = sel.GetIndex(); - Vec<3> pcenter = 0.0; - for(auto i : IntRange(sel.GetNP())) - { - for(auto d : IntRange(3)) - pcenter[d] += mesh[sel[i]][d]; - } - pcenter = 1.0/sel.GetNP() * pcenter; + auto & seg = line_segments[si]; auto n = mesh[seg[1]] - mesh[seg[0]]; n = {-n[1], n[0], 0}; n.Normalize(); - Vec<3> p0{mesh[seg[0]]}; - Vec<3> p1{mesh[seg[0]]}; - - - auto v = pcenter -0.5*(p0+p1); - - const auto Dot = [](Vec<3> a, Vec<3> b) - { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }; - if(Dot(n, v)<0) - n = -1*n; + if(seg.domout == domain) + n = -n; AddDirection(growthvectors[seg[0]], n); AddDirection(growthvectors[seg[1]], n); @@ -1628,7 +1571,7 @@ namespace netgen for(auto si : moved_segs) { - auto current_seg = mesh[si]; + auto current_seg = line_segments[si]; auto current_si = si; auto first = current_seg[0]; @@ -1766,8 +1709,9 @@ namespace netgen const auto & seg0 = mesh[segi0]; const auto & seg1 = mesh[segi1]; - if(seg0.si != seg1.si) - return; + if( (seg0.domin != domain && seg0.domout != domain) || + (seg1.domin != domain && seg1.domout != domain) ) + return; if(segi0 == segi1) return; @@ -1880,7 +1824,7 @@ namespace netgen { auto p2 = [](Point<3> p) { return Point<2>{p[0], p[1]}; }; - auto seg = mesh[segi]; + auto seg = line_segments[segi]; double alpha,beta; intersect( p2(mesh[seg[0]]), p2(mesh[seg[0]]+total_thickness*growthvectors[seg[0]]), p2(mesh[seg[1]]), p2(mesh[seg[1]]+total_thickness*growthvectors[seg[1]]), alpha, beta ); @@ -1918,13 +1862,13 @@ namespace netgen // insert new elements ( and move old ones ) for(auto si : moved_segs) { - auto seg = mesh[si]; + auto seg = line_segments[si]; bool swap = false; auto & pm0 = mapto[seg[0]]; auto & pm1 = mapto[seg[1]]; - auto newindex = si_map[seg.si]; + auto newindex = si_map[domain]; Segment s = seg; s.geominfo[0] = {}; @@ -1939,10 +1883,6 @@ namespace netgen s.si = seg.si; mesh.AddSegment(s); - Swap(s[0], s[1]); - s.si = newindex; - mesh.AddSegment(s); - for ( auto i : Range(thicknesses)) { PointIndex pi0, pi1, pi2, pi3; @@ -1980,14 +1920,14 @@ namespace netgen newel[1] = pi1; newel[2] = pi2; newel[3] = pi3; - newel.SetIndex(si_map[seg.si]); + newel.SetIndex(new_domain); newel.GeomInfo() = PointGeomInfo{}; -// if(swap) -// { -// Swap(newel[0], newel[1]); -// Swap(newel[2], newel[3]); -// } + if(swap) + { + Swap(newel[0], newel[1]); + Swap(newel[2], newel[3]); + } for(auto i : Range(4)) { @@ -1998,7 +1938,10 @@ namespace netgen } // segment now adjacent to new 2d-domain! - mesh[si].si = si_map[seg.si]; + if(line_segments[si].domin == domain) + line_segments[si].domin = new_domain; + if(line_segments[si].domout == domain) + line_segments[si].domout = new_domain; } for(auto pi : Range(mapto)) @@ -2043,7 +1986,12 @@ namespace netgen } for(auto segi : moved_segs) - mesh[segi].si = domain; + { + if(mesh[segi].domin == new_domain) + mesh[segi].domin = domain; + if(mesh[segi].domout == new_domain) + mesh[segi].domout = domain; + } mesh.Compress(); mesh.CalcSurfacesOfNode(); diff --git a/libsrc/meshing/fieldlines.cpp b/libsrc/meshing/fieldlines.cpp index 0db30a6e..84cf0a65 100644 --- a/libsrc/meshing/fieldlines.cpp +++ b/libsrc/meshing/fieldlines.cpp @@ -77,7 +77,7 @@ namespace netgen { bool finished = false; - if(stepcount <= steps) + if(stepcount <= steps && stepcount>0) { t = startt + c[stepcount-1]*h; val = startval; diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 01f547f0..192d98f5 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1191,6 +1191,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard()) + .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array{}) .def ("BoundaryLayer", [](Mesh & self, variant boundary, variant thickness, string material, diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index eb598f1f..7d719f54 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -202,6 +202,12 @@ namespace netgen Mesh::GEOM_TYPE GetGeomType() const override { return Mesh::GEOM_OCC; } + void SetDimension(int dim) + { + dimension = dim; + BuildFMap(); + } + void SetOCCParameters(const OCCParameters& par) { occparam = par; } diff --git a/libsrc/occ/python_occ.cpp b/libsrc/occ/python_occ.cpp index 8247d352..b74d12ac 100644 --- a/libsrc/occ/python_occ.cpp +++ b/libsrc/occ/python_occ.cpp @@ -134,7 +134,7 @@ DLL_HEADER void ExportNgOCC(py::module &m) }), py::arg("shape"), "Create Netgen OCCGeometry from existing TopoDS_Shape") - .def(py::init([] (const string& filename) + .def(py::init([] (const string& filename, int dim) { shared_ptr geo; if(EndsWith(filename, ".step") || EndsWith(filename, ".stp")) @@ -145,9 +145,11 @@ DLL_HEADER void ExportNgOCC(py::module &m) geo.reset(LoadOCC_IGES(filename)); else throw Exception("Cannot load file " + filename + "\nValid formats are: step, stp, brep, iges"); + if(dim<3) + geo->SetDimension(dim); ng_geometry = geo; return geo; - }), py::arg("filename"), + }), py::arg("filename"), py::arg("dim")=3, "Load OCC geometry from step, brep or iges file") .def(NGSPickle()) .def("Glue", &OCCGeometry::GlueGeometry)