From ea7e980c8d25f02f5ee597a505d8488abbf14f5b Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 23 Feb 2022 12:22:19 +0100 Subject: [PATCH] [occ] ZRefinement --- libsrc/general/hashtabl.hpp | 37 +++++++- libsrc/general/table.hpp | 2 + libsrc/meshing/meshclass.cpp | 162 +++++++++++++++++++++++++++++++++ libsrc/meshing/meshclass.hpp | 3 +- libsrc/meshing/python_mesh.cpp | 2 + 5 files changed, 201 insertions(+), 5 deletions(-) diff --git a/libsrc/general/hashtabl.hpp b/libsrc/general/hashtabl.hpp index 5c5fe4f8..7841d02b 100644 --- a/libsrc/general/hashtabl.hpp +++ b/libsrc/general/hashtabl.hpp @@ -149,7 +149,14 @@ public: int pos = Position (bnr, ahash); return cont.Get (bnr, pos); } - + + T & Get (const INDEX_2 & ahash) + { + int bnr = HashValue (ahash); + int pos = Position (bnr, ahash); + return cont.Get (bnr, pos); + } + /// bool Used (const INDEX_2 & ahash) const { @@ -214,9 +221,14 @@ public: int BagNr() const { return bagnr; } int Pos() const { return pos; } - void operator++ (int) + Iterator operator++ (int) + { + Iterator it(ht, bagnr, pos); + ++(*this); + return it; + } + Iterator& operator++() { - // cout << "begin Operator ++: bagnr = " << bagnr << " - pos = " << pos << endl; pos++; while (bagnr < ht.GetNBags() && pos == ht.GetBagSize(bagnr+1)) @@ -224,7 +236,12 @@ public: pos = 0; bagnr++; } - // cout << "end Operator ++: bagnr = " << bagnr << " - pos = " << pos << endl; + return *this; + } + + std::pair operator*() + { + return std::make_pair(ht.hash[bagnr][pos], ht.cont[bagnr][pos]); } bool operator != (int i) const @@ -246,6 +263,18 @@ public: return GetNBags(); } + Iterator begin () const + { + Iterator it(*this, 0, -1); + it++; + return it; + } + + int end() const + { + return GetNBags(); + } + void GetData (const Iterator & it, INDEX_2 & ahash, T & acont) const { diff --git a/libsrc/general/table.hpp b/libsrc/general/table.hpp index 832fc751..e07f4d8a 100644 --- a/libsrc/general/table.hpp +++ b/libsrc/general/table.hpp @@ -201,6 +201,8 @@ public: inline const T & Get (int i, int nr) const { return ((T*)data.Get(i).col)[nr-1]; } + inline T & Get (int i, int nr) + { return ((T*)data.Get(i).col)[nr-1]; } /** Returns pointer to the first element in row i. */ inline const T * GetLine (int i) const diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 945b9ea8..ed0a11a0 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -6213,7 +6213,169 @@ namespace netgen */ } + void Mesh :: ZRefine(const string& name, const Array& slices) + { + auto nr = GetIdentifications().GetNr(name); + auto& identpts = GetIdentifications().GetIdentifiedPoints(); + UpdateTopology(); + + std::map, + Array> inserted_points; + BitArray mapped_points(GetNV()+1); + mapped_points = false; + + // Add new points + for(auto [p1p2, idnr] : identpts) + { + if(idnr != nr) + continue; + auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}]; + auto p1 = Point(p1p2.I1()); + auto p2 = Point(p1p2.I2()); + ipts.Append(p1p2.I1()); + mapped_points.SetBit(p1p2.I1()); + for(auto slice : slices) + { + auto np = p1 + slice * (p2-p1); + auto npi = AddPoint(np); + ipts.Append(npi); + } + ipts.Append(p1p2.I2()); + } + + // Split segments + for(auto si : Range(segments)) + { + auto& seg = segments[si]; + auto p1 = seg[0]; + auto p2 = seg[1]; + + auto c1 = inserted_points.count({p1, p2}); + auto c2 = inserted_points.count({p2, p1}); + + if(c1 == 0 && c2 == 0) + continue; + + if(c2) + Swap(p1,p2); + + const auto& ipts = inserted_points[{p1,p2}]; + if(c2) + seg[1] = ipts[ipts.Size()-2]; + else + seg[1] = ipts[1]; + for(auto i : Range(size_t(1), ipts.Size()-1)) + { + Segment snew = seg; + if(c2) + { + seg[0] = ipts[ipts.Size()-1-i]; + seg[1] = ipts[ipts.Size()-2-i]; + } + else + { + snew[0] = ipts[i]; + snew[1] = ipts[i+1]; + } + AddSegment(snew); + } + } + + BitArray sel_done(surfelements.Size()); + sel_done = false; + + // Split surface elements + auto p2sel = CreatePoint2SurfaceElementTable(); + for(const auto& [pair, inserted] : inserted_points) + { + for(auto si : p2sel[pair.first]) + { + if(sel_done[si]) + continue; + sel_done.SetBit(si); + auto sel = surfelements[si]; + map> mapped_points; + int nmapped = 0; + for(auto i : Range(sel.GetNP())) + { + auto p1 = sel[i]; + auto p2 = sel[(i+1)%sel.GetNP()]; + auto c1 = inserted_points.count({p1, p2}); + auto c2 = inserted_points.count({p2, p1}); + if(c1 == 0 && c2 == 0) + continue; + if(c2) + Swap(p1, p2); + auto& ipts = inserted_points[{p1, p2}]; + auto& a1 = mapped_points[p1]; + auto& a2 = mapped_points[p2]; + a1 = ipts.Range(0, ipts.Size()-1); + a2 = ipts.Range(1, ipts.Size()); + nmapped = ipts.Size()-1; + } + for(auto i : Range(nmapped)) + { + Element2d nsel = sel; + for(auto& pi : nsel.PNums()) + if(mapped_points.count(pi)) + pi = mapped_points[pi][i]; + AddSurfaceElement(nsel); + } + if(nmapped) + surfelements[si].Delete(); + } + } + + // Split volume elements + BitArray vol_done(volelements.Size()); + vol_done = false; + auto p2el = CreatePoint2ElementTable(); // mapped_points); + for(const auto& [pair, inserted] : inserted_points) + { + for(auto ei : p2el[pair.first]) + { + if(vol_done[ei]) + continue; + vol_done.SetBit(ei); + auto el = volelements[ei]; + map> mapped_points; + int nmapped = 0; + NgArray eledges; + topology.GetElementEdges(ei+1, eledges); + for(auto edgei : eledges) + { + int p1, p2; + topology.GetEdgeVertices(edgei, p1, p2); + auto c1 = inserted_points.count({p1, p2}); + auto c2 = inserted_points.count({p2, p1}); + if(c1 == 0 && c2 == 0) + continue; + if(c2) + Swap(p1, p2); + auto& ipts = inserted_points[{p1, p2}]; + auto& a1 = mapped_points[p1]; + auto& a2 = mapped_points[p2]; + a1 = ipts.Range(0, ipts.Size()-1); + a2 = ipts.Range(1, ipts.Size()); + nmapped = ipts.Size()-1; + } + + for(auto i : Range(nmapped)) + { + Element nel = el; + for(auto& pi : nel.PNums()) + if(mapped_points.count(pi)) + pi = mapped_points[pi][i]; + AddVolumeElement(nel); + } + if(nmapped) + volelements[ei].Delete(); + } + } + + Compress(); + } void Mesh :: RebuildSurfaceElementLists () { diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 8a2a8166..f16ba622 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -498,7 +498,8 @@ namespace netgen /// Refines mesh and projects points to true surface // void Refine (int levels, const CSGeometry * geom); - + + void ZRefine(const string& name, const Array& slices); bool BoundaryEdge (PointIndex pi1, PointIndex pi2) const { diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index eb71a381..f335565e 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1127,6 +1127,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) self.UpdateTopology(); }),py::call_guard()) + .def("ZRefine", &Mesh::ZRefine) + .def ("SecondOrder", FunctionPointer ([](Mesh & self) {