From 73bcb1bd297ab4a39f06326d92a85e027d68b4ba Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Mon, 23 Dec 2024 19:24:48 +0100 Subject: [PATCH] PointIndex in bisect --- libsrc/core/bitarray.hpp | 18 ++++ libsrc/include/nginterface_v2_impl.hpp | 9 +- libsrc/interface/nginterface.cpp | 4 +- libsrc/meshing/bisect.cpp | 128 +++++++++++++------------ libsrc/meshing/meshclass.hpp | 2 +- libsrc/meshing/meshtype.hpp | 3 +- libsrc/meshing/parallelmesh.cpp | 2 +- libsrc/meshing/paralleltop.cpp | 2 +- libsrc/meshing/validate.cpp | 4 +- libsrc/meshing/validate.hpp | 4 +- 10 files changed, 99 insertions(+), 77 deletions(-) diff --git a/libsrc/core/bitarray.hpp b/libsrc/core/bitarray.hpp index 91259e91..a354952b 100644 --- a/libsrc/core/bitarray.hpp +++ b/libsrc/core/bitarray.hpp @@ -210,6 +210,24 @@ private: NGCORE_API std::ostream & operator<<(std::ostream & s, const BitArray & ba); + + + template + class TBitArray : public BitArray + { + public: + using BitArray::BitArray; + + void SetBit (IndexType i) { BitArray::SetBit(i-IndexBASE()); } + void Clear () { BitArray::Clear(); } + void Clear (IndexType i) { BitArray::Clear(i-IndexBASE()); } + void SetBitAtomic (IndexType i) { BitArray::SetBitAtomic(i-IndexBASE()); } + bool Test (IndexType i) const { return BitArray::Test(i-IndexBASE()); } + bool operator[] (IndexType i) const { return Test(i); } + + }; + } // namespace ngcore + #endif // NETGEN_CORE_BITARRAY diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index facd6c81..e3783ea8 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -349,12 +349,9 @@ NGX_INLINE DLL_HEADER Ng_Buffer Ngx_Mesh :: GetPeriodicVertices(int idnr NGX_INLINE void Ngx_Mesh :: GetParentNodes (int ni, int * parents) const { - ni++; - if (ni <= mesh->mlbetweennodes.Size()) - { - parents[0] = mesh->mlbetweennodes.Get(ni).I1()-1; - parents[1] = mesh->mlbetweennodes.Get(ni).I2()-1; - } + if (ni < mesh->mlbetweennodes.Size()) + for (int j = 0; j < 2; j++) + parents[j] = mesh->mlbetweennodes[IndexBASE()+ni][j] - IndexBASE(); else parents[0] = parents[1] = -1; } diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index 083418df..ea66ef41 100644 --- a/libsrc/interface/nginterface.cpp +++ b/libsrc/interface/nginterface.cpp @@ -1736,8 +1736,8 @@ void Ng_GetParentNodes (int ni, int * parents) { if (ni <= mesh->mlbetweennodes.Size()) { - parents[0] = mesh->mlbetweennodes.Get(ni).I1(); - parents[1] = mesh->mlbetweennodes.Get(ni).I2(); + parents[0] = mesh->mlbetweennodes[ni].I1(); + parents[1] = mesh->mlbetweennodes[ni].I2(); } else parents[0] = parents[1] = 0; diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index b0347e23..70de2010 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -3230,7 +3230,7 @@ namespace netgen if (opt.refine_hp) { PrintMessage(3,"refine hp"); - NgBitArray singv(np); + TBitArray singv(np); singv.Clear(); if (mesh.GetDimension() == 3) @@ -3238,8 +3238,8 @@ namespace netgen for (int i = 1; i <= mesh.GetNSeg(); i++) { const Segment & seg = mesh.LineSegment(i); - singv.Set (seg[0]); - singv.Set (seg[1]); + singv.SetBit (seg[0]); + singv.SetBit (seg[1]); } /* for ( i=1; i<= mesh.GetNSE(); i++) @@ -3253,18 +3253,18 @@ namespace netgen else { // vertices with 2 different bnds - NgArray bndind(np); + Array bndind(np); bndind = 0; for (int i = 1; i <= mesh.GetNSeg(); i++) { const Segment & seg = mesh.LineSegment(i); for (int j = 0; j < 2; j++) { - int pi = (j == 0) ? seg[0] : seg[1]; - if (bndind.Elem(pi) == 0) - bndind.Elem(pi) = seg.edgenr; - else if (bndind.Elem(pi) != seg.edgenr) - singv.Set (pi); + PointIndex pi = (j == 0) ? seg[0] : seg[1]; + if (bndind[pi] == 0) + bndind[pi] = seg.edgenr; + else if (bndind[pi] != seg.edgenr) + singv.SetBit (pi); } } } @@ -3339,8 +3339,8 @@ namespace netgen { MarkedTet oldtet = mtets[i]; - INDEX_2 edge(oldtet.pnums[oldtet.tetedge1], - oldtet.pnums[oldtet.tetedge2]); + PointIndices<2> edge(oldtet.pnums[oldtet.tetedge1], + oldtet.pnums[oldtet.tetedge2]); edge.Sort(); PointIndex newp; @@ -3350,8 +3350,7 @@ namespace netgen } else { - Point<3> npt = Center (mesh.Point (edge.I1()), - mesh.Point (edge.I2())); + Point<3> npt = Center (mesh[edge[0]], mesh[edge[1]]); newp = mesh.AddPoint (npt); cutedges.Set (edge, newp); } @@ -3378,11 +3377,11 @@ namespace netgen if (pi1 == oldprism.markededge) pi1++; int pi2 = 3-pi1-oldprism.markededge; - - INDEX_2 edge1(oldprism.pnums[pi1], - oldprism.pnums[pi2]); - INDEX_2 edge2(oldprism.pnums[pi1+3], - oldprism.pnums[pi2+3]); + + PointIndices<2> edge1(oldprism.pnums[pi1], + oldprism.pnums[pi2]); + PointIndices<2> edge2(oldprism.pnums[pi1+3], + oldprism.pnums[pi2+3]); edge1.Sort(); edge2.Sort(); @@ -3390,8 +3389,7 @@ namespace netgen newp1 = cutedges.Get(edge1); else { - Point<3> npt = Center (mesh.Point (edge1.I1()), - mesh.Point (edge1.I2())); + Point<3> npt = Center (mesh[edge1[0]], mesh[edge1[1]]); newp1 = mesh.AddPoint (npt); cutedges.Set (edge1, newp1); } @@ -3399,8 +3397,7 @@ namespace netgen newp2 = cutedges.Get(edge2); else { - Point<3> npt = Center (mesh.Point (edge2.I1()), - mesh.Point (edge2.I2())); + Point<3> npt = Center (mesh[edge2[0]], mesh[edge2[1]]); newp2 = mesh.AddPoint (npt); cutedges.Set (edge2, newp2); } @@ -3422,18 +3419,22 @@ namespace netgen oldid = mids.Get(i); - NgArray edges; - edges.Append(INDEX_2(oldid.pnums[oldid.markededge], - oldid.pnums[(oldid.markededge+1)%oldid.np])); - edges.Append(INDEX_2(oldid.pnums[oldid.markededge + oldid.np], - oldid.pnums[(oldid.markededge+1)%oldid.np + oldid.np])); - + NgArray> edges; + edges.Append( { + oldid.pnums[oldid.markededge], + oldid.pnums[(oldid.markededge+1)%oldid.np] } ); + edges.Append( { + oldid.pnums[oldid.markededge + oldid.np], + oldid.pnums[(oldid.markededge+1)%oldid.np + oldid.np] } ); + if(oldid.np == 4) { - edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np], - oldid.pnums[(oldid.markededge+3)%oldid.np])); - edges.Append(INDEX_2(oldid.pnums[(oldid.markededge+2)%oldid.np + oldid.np], - oldid.pnums[(oldid.markededge+3)%oldid.np + oldid.np])); + edges.Append( { + oldid.pnums[(oldid.markededge+2)%oldid.np], + oldid.pnums[(oldid.markededge+3)%oldid.np]} ); + edges.Append( { + oldid.pnums[(oldid.markededge+2)%oldid.np + oldid.np], + oldid.pnums[(oldid.markededge+3)%oldid.np + oldid.np] } ); } for (int j = 0; j < edges.Size(); j++) { @@ -3479,7 +3480,7 @@ namespace netgen MarkedTri oldtri = mtris[i]; PointIndex oldpi1 = oldtri.pnums[(oldtri.markededge+1)%3]; PointIndex oldpi2 = oldtri.pnums[(oldtri.markededge+2)%3]; - INDEX_2 edge(oldpi1, oldpi2); + PointIndices<2> edge(oldpi1, oldpi2); edge.Sort(); int si = mesh.GetFaceDescriptor (oldtri.surfid).SurfNr(); @@ -3529,21 +3530,21 @@ namespace netgen INDEX_2 edge2(oldquad.pnums[2], oldquad.pnums[3]); */ - INDEX_2 edge1, edge2; + PointIndices<2> edge1, edge2; PointGeomInfo pgi11, pgi12, pgi21, pgi22; if (oldquad.markededge==0 || oldquad.markededge==2) { - edge1.I1()=oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; - edge1.I2()=oldquad.pnums[1]; pgi12=oldquad.pgeominfo[1]; - edge2.I1()=oldquad.pnums[2]; pgi21=oldquad.pgeominfo[2]; - edge2.I2()=oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; + edge1[0] = oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; + edge1[1] = oldquad.pnums[1]; pgi12=oldquad.pgeominfo[1]; + edge2[0] = oldquad.pnums[2]; pgi21=oldquad.pgeominfo[2]; + edge2[1] = oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; } else // 3 || 1 { - edge1.I1()=oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; - edge1.I2()=oldquad.pnums[2]; pgi12=oldquad.pgeominfo[2]; - edge2.I1()=oldquad.pnums[1]; pgi21=oldquad.pgeominfo[1]; - edge2.I2()=oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; + edge1[0] = oldquad.pnums[0]; pgi11=oldquad.pgeominfo[0]; + edge1[1] = oldquad.pnums[2]; pgi12=oldquad.pgeominfo[2]; + edge2[0] = oldquad.pnums[1]; pgi21=oldquad.pgeominfo[1]; + edge2[1] = oldquad.pnums[3]; pgi22=oldquad.pgeominfo[3]; } edge1.Sort(); @@ -3555,8 +3556,8 @@ namespace netgen } else { - Point<3> np1 = Center (mesh.Point (edge1.I1()), - mesh.Point (edge1.I2())); + Point<3> np1 = Center (mesh.Point (edge1[0]), + mesh.Point (edge1[1])); newp1 = mesh.AddPoint (np1); cutedges.Set (edge1, newp1); } @@ -3567,8 +3568,8 @@ namespace netgen } else { - Point<3> np2 = Center (mesh.Point (edge2.I1()), - mesh.Point (edge2.I2())); + Point<3> np2 = Center (mesh.Point (edge2[0]), + mesh.Point (edge2[1])); newp2 = mesh.AddPoint (np2); cutedges.Set (edge2, newp2); } @@ -3606,7 +3607,7 @@ namespace netgen for (int i = 1; i <= nseg; i++) { Segment & seg = mesh.LineSegment (i); - INDEX_2 edge(seg[0], seg[1]); + PointIndices<2> edge(seg[0], seg[1]); edge.Sort(); if (cutedges.Used (edge)) { @@ -3614,7 +3615,7 @@ namespace netgen Segment nseg1 = seg; Segment nseg2 = seg; - int newpi = cutedges.Get(edge); + PointIndex newpi = cutedges.Get(edge); nseg1[1] = newpi; nseg2[0] = newpi; @@ -3670,7 +3671,7 @@ namespace netgen if (opt.refine_hp) { // - NgArray v_order (mesh.GetNP()); + Array v_order (mesh.GetNP()); v_order = 0; if (mesh.GetDimension() == 3) { @@ -3680,12 +3681,12 @@ namespace netgen for (int i = 0; i < mtets.Size(); i++) for (int j = 0; j < 4; j++) - if (int(mtets[i].order) > v_order.Elem(mtets[i].pnums[j])) - v_order.Elem(mtets[i].pnums[j]) = mtets[i].order; + if (int(mtets[i].order) > v_order[mtets[i].pnums[j]]) + v_order[mtets[i].pnums[j]] = mtets[i].order; for (int i = 0; i < mtets.Size(); i++) for (int j = 0; j < 4; j++) - if (int(mtets[i].order) < v_order.Elem(mtets[i].pnums[j])-1) - mtets[i].order = v_order.Elem(mtets[i].pnums[j])-1; + if (int(mtets[i].order) < v_order[mtets[i].pnums[j]]-1) + mtets[i].order = v_order[mtets[i].pnums[j]]-1; } else { @@ -3697,13 +3698,13 @@ namespace netgen for (int i = 0; i < mtris.Size(); i++) for (int j = 0; j < 3; j++) - if (int(mtris[i].order) > v_order.Elem(mtris[i].pnums[j])) - v_order.Elem(mtris[i].pnums[j]) = mtris[i].order; + if (int(mtris[i].order) > v_order[mtris[i].pnums[j]]) + v_order[mtris[i].pnums[j]] = mtris[i].order; for (int i = 0; i < mtris.Size(); i++) { for (int j = 0; j < 3; j++) - if (int(mtris[i].order) < v_order.Elem(mtris[i].pnums[j])-1) - mtris[i].order = v_order.Elem(mtris[i].pnums[j])-1; + if (int(mtris[i].order) < v_order[mtris[i].pnums[j]]-1) + mtris[i].order = v_order[mtris[i].pnums[j]]-1; } } } @@ -3866,11 +3867,15 @@ namespace netgen if (mesh.level_nv.Size() <= 1) { PrintMessage(4,"RESETTING mlbetweennodes"); + /* for (int i = 1; i <= np; i++) { mesh.mlbetweennodes.Elem(i).I1() = 0; mesh.mlbetweennodes.Elem(i).I2() = 0; } + */ + for (auto i : mesh.mlbetweennodes.Range()) + mesh.mlbetweennodes[i] = { PointIndex::INVALID, PointIndex::INVALID }; } mesh.level_nv.Append (np); @@ -3888,7 +3893,7 @@ namespace netgen */ - NgBitArray isnewpoint(np); + TBitArray isnewpoint(np); isnewpoint.Clear(); { @@ -3899,8 +3904,8 @@ namespace netgen INDEX_2 edge; PointIndex newpi; cutedges.GetData0 (i, edge, newpi); - isnewpoint.Set(newpi); - mesh.mlbetweennodes.Elem(newpi) = edge; + isnewpoint.SetBit(newpi); + mesh.mlbetweennodes[newpi] = edge; } } /* @@ -4034,7 +4039,8 @@ namespace netgen if(noprojection) { do_repair = false; - for(int ii=1; ii<=mesh.GetNP(); ii++) + // for(int ii=1; ii<=mesh.GetNP(); ii++) + for (auto ii : mesh.Points().Range()) { if(isnewpoint.Test(ii) && mesh.mlbetweennodes[ii][0].IsValid()) { diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 26d71127..87bd6f99 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -222,7 +222,7 @@ namespace netgen // number of vertices on each refinement level: NgArray level_nv; /// refinement hierarchy - NgArray,PointIndex::BASE> mlbetweennodes; + Array,PointIndex> mlbetweennodes; /// parent element of volume element NgArray mlparentelement; /// parent element of surface element diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 02355270..939b1e7b 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -217,11 +217,11 @@ namespace netgen // throw Exception("illegal PointIndex, use PointIndex::INVALID instead"); #endif } + /* friend constexpr netgen::PointIndex ngcore::IndexBASE (); friend istream & operator>> (istream &, PointIndex &); friend ostream & operator<< (ostream &, const PointIndex &); template friend class PointIndices; - /* friend PointIndex operator+ (PointIndex, int); friend PointIndex operator+ (PointIndex, size_t); friend PointIndex operator+ (int, PointIndex); @@ -234,6 +234,7 @@ namespace netgen friend bool operator== (PointIndex a, PointIndex b); friend bool operator!= (PointIndex a, PointIndex b); */ + public: constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; } // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index 3e3ac61d..a8b5b63b 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -955,7 +955,7 @@ namespace netgen self.openelements = NgArray(0); self.opensegments = NgArray(0); self.numvertices = 0; - self.mlbetweennodes = NgArray,PointIndex::BASE> (0); + self.mlbetweennodes = Array,PointIndex> (0); self.mlparentelement = NgArray(0); self.mlparentsurfaceelement = NgArray(0); self.curvedelems = make_unique (self); diff --git a/libsrc/meshing/paralleltop.cpp b/libsrc/meshing/paralleltop.cpp index a2107726..7462f5ad 100644 --- a/libsrc/meshing/paralleltop.cpp +++ b/libsrc/meshing/paralleltop.cpp @@ -184,7 +184,7 @@ namespace netgen if (mesh.mlbetweennodes.Size() == mesh.Points().Size()) { - NgArray,PointIndex::BASE> hml { mesh.mlbetweennodes }; + Array,PointIndex> hml { mesh.mlbetweennodes }; for (PointIndex pi : Range(mesh.Points())) mesh.mlbetweennodes[inv_index[pi]] = hml[pi]; } diff --git a/libsrc/meshing/validate.cpp b/libsrc/meshing/validate.cpp index 31b4abfe..781e6cd0 100644 --- a/libsrc/meshing/validate.cpp +++ b/libsrc/meshing/validate.cpp @@ -6,7 +6,7 @@ namespace netgen { void GetPureBadness(Mesh & mesh, NgArray & pure_badness, - const NgBitArray & isnewpoint) + const TBitArray & isnewpoint) { //const int ne = mesh.GetNE(); const int np = mesh.GetNP(); @@ -152,7 +152,7 @@ namespace netgen void RepairBisection(Mesh & mesh, NgArray & bad_elements, - const NgBitArray & isnewpoint, const Refinement & refinement, + const TBitArray & isnewpoint, const Refinement & refinement, const NgArray & pure_badness, double max_worsening, const bool uselocalworsening, const NgArray< idmap_type* > & idmaps) diff --git a/libsrc/meshing/validate.hpp b/libsrc/meshing/validate.hpp index ae61a133..3defaf02 100644 --- a/libsrc/meshing/validate.hpp +++ b/libsrc/meshing/validate.hpp @@ -5,13 +5,13 @@ namespace netgen { void GetPureBadness(Mesh & mesh, NgArray & pure_badness, - const NgBitArray & isnewpoint); + const TBitArray & isnewpoint); double Validate(const Mesh & mesh, NgArray & bad_elements, const NgArray & pure_badness, double max_worsening, const bool uselocalworsening, NgArray * quality_loss = NULL); void RepairBisection(Mesh & mesh, NgArray & bad_elements, - const NgBitArray & isnewpoint, const Refinement & refinement, + const TBitArray & isnewpoint, const Refinement & refinement, const NgArray & pure_badness, double max_worsening, const bool uselocalworsening, const NgArray< idmap_type* > & idmaps);