From 00edc92c000f875110f6d82961367f623ea46249 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sat, 28 Dec 2024 13:01:20 +0100 Subject: [PATCH] improve3 with consistent PointIndex --- libsrc/meshing/improve3.cpp | 50 +++++++++++++++++------------------ libsrc/meshing/improve3.hpp | 8 +++--- libsrc/meshing/meshclass.cpp | 40 ++++++++++++++-------------- libsrc/meshing/meshtype.hpp | 14 +++++++--- libsrc/meshing/smoothing3.cpp | 2 +- 5 files changed, 61 insertions(+), 53 deletions(-) diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index c081e6ba..6b285591 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -100,7 +100,7 @@ static double SplitElementBadness (const Mesh::T_POINTS & points, const MeshingP { double badness = 0; auto np = old.GetNP(); - PointIndex dummy{-1}; + PointIndex dummy{PointIndex::INVALID}; if(np == 4) { // Split tet into two tets @@ -451,7 +451,7 @@ void MeshOptimize3d :: CombineImprove () -double MeshOptimize3d :: SplitImproveEdge (Table & elementsonnode, NgArray &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only) +double MeshOptimize3d :: SplitImproveEdge (Table & elementsonnode, NgArray> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only) { double d_badness = 0.0; // int cnt = 0; @@ -512,11 +512,11 @@ double MeshOptimize3d :: SplitImproveEdge (Table & elem for (int l = 0; l < 4; l++) if (el[l] == pi1 || el[l] == pi2) { - INDEX_3 i3; + PointIndices<3> i3; Element2d face(TRIG); el.GetFace (l+1, face); - for (int kk = 1; kk <= 3; kk++) - i3.I(kk) = face.PNum(kk); + for (int kk = 0; kk < 3; kk++) + i3[kk] = face[kk]; locfaces.Append (i3); } } @@ -536,7 +536,7 @@ double MeshOptimize3d :: SplitImproveEdge (Table & elem { int pok = pf.Func (px) < 1e10; if (!pok) - pok = FindInnerPoint (mesh.Points(), locfaces, pnew); + pok = FindInnerPoint (mesh.Points(), locfaces, pnew); if(pok) { @@ -647,7 +647,7 @@ void MeshOptimize3d :: SplitImprove () tsearch.Start(); ParallelForRange(Range(edges), [&] (auto myrange) { - NgArray locfaces; + NgArray> locfaces; for(auto i : myrange) { @@ -671,7 +671,7 @@ void MeshOptimize3d :: SplitImprove () // Apply actual optimizations topt.Start(); int cnt = 0; - NgArray locfaces; + NgArray> locfaces; for(auto [d_badness, ei] : edges_with_improvement) { auto [p0,p1] = edges[ei]; @@ -871,7 +871,7 @@ double MeshOptimize3d :: SwapImproveEdge ( if ((goal == OPT_CONFORM) && NotTooBad(bad1, bad2)) { - INDEX_3 face(pi3, pi4, pi5); + PointIndices<3> face(pi3, pi4, pi5); face.Sort(); if (faces.Used(face)) { @@ -1237,9 +1237,9 @@ double MeshOptimize3d :: SwapImproveEdge ( for (int k = l+1; k <= nsuround + l - 2; k++) { - INDEX_3 hi3(suroundpts[l], - suroundpts[k % nsuround], - suroundpts[(k+1) % nsuround]); + PointIndices<3> hi3(suroundpts[l], + suroundpts[k % nsuround], + suroundpts[(k+1) % nsuround]); hi3.Sort(); if (faces.Used(hi3)) { @@ -1336,7 +1336,7 @@ void MeshOptimize3d :: SwapImprove (const BitArray * working_elements) // int ne = mesh.GetNE(); mesh.BuildBoundaryEdges(false); - BitArray free_points(mesh.GetNP()+PointIndex::BASE); + TBitArray free_points(mesh.GetNP()); free_points.Clear(); ParallelForRange(mesh.VolumeElements().Range(), [&] (auto myrange) @@ -1372,7 +1372,7 @@ void MeshOptimize3d :: SwapImprove (const BitArray * working_elements) for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & hel = mesh.OpenElement(i); - INDEX_3 face(hel[0], hel[1], hel[2]); + PointIndices<3> face(hel[0], hel[1], hel[2]); face.Sort(); faces.Set (face, i); } @@ -1439,7 +1439,7 @@ void MeshOptimize3d :: SwapImprove (const BitArray * working_elements) { Element2d sel; el.GetFace(i, sel); - INDEX_3 face(sel[0], sel[1], sel[2]); + PointIndices<3> face(sel[0], sel[1], sel[2]); face.Sort(); if(faces.Used(face)) open_els[faces.Get(face)-1].Delete(); @@ -1502,9 +1502,9 @@ void MeshOptimize3d :: SwapImproveSurface ( // contains at least all elements at node - TABLE elementsonnode(np); - TABLE surfaceelementsonnode(np); - TABLE surfaceindicesonnode(np); + DynamicTable elementsonnode(np); + DynamicTable surfaceelementsonnode(np); + DynamicTable surfaceindicesonnode(np); NgArray hasbothpoints; NgArray hasbothpointsother; @@ -1607,14 +1607,14 @@ void MeshOptimize3d :: SwapImproveSurface ( othermattype = -1; - INDEX_2 i2 (pi1, pi2); + PointIndices<2> i2 (pi1, pi2); i2.Sort(); if (edgeused.Used(i2)) continue; edgeused.Set (i2, 1); if(periodic) { - i2.I1() = pi1other; - i2.I2() = pi2other; + i2[0] = pi1other; + i2[1] = pi2other; i2.Sort(); edgeused.Set(i2,1); } @@ -1770,7 +1770,7 @@ void MeshOptimize3d :: SwapImproveSurface ( //(*testout) << "sel1 " << sel1 << " sel2 " << sel2 << " el " << mesh[sel1] << " resp. " << mesh[sel2] << endl; - PointIndex sp1(0), sp2(0); + PointIndex sp1(PointIndex::INVALID), sp2(PointIndex::INVALID); PointIndex sp1other, sp2other; for(int l=0; l & elementsonnode, - TABLE & belementsonnode, bool check_only ) + Table & elementsonnode, + DynamicTable & belementsonnode, bool check_only ) { PointIndex pi1, pi2, pi3, pi4, pi5; Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); @@ -2509,7 +2509,7 @@ void MeshOptimize3d :: SwapImprove2 () int nse = mesh.GetNSE(); // contains at least all elements at node - TABLE belementsonnode(np); + DynamicTable belementsonnode(np); PrintMessage (3, "SwapImprove2 "); (*testout) << "\n" << "Start SwapImprove2" << "\n"; diff --git a/libsrc/meshing/improve3.hpp b/libsrc/meshing/improve3.hpp index 1f02df22..bfbf89c0 100644 --- a/libsrc/meshing/improve3.hpp +++ b/libsrc/meshing/improve3.hpp @@ -35,7 +35,7 @@ public: void CombineImprove (); void SplitImprove (); - double SplitImproveEdge (Table & elementsonnode, NgArray &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false); + double SplitImproveEdge (Table & elementsonnode, NgArray> &locfaces, double badmax, PointIndex pi1, PointIndex pi2, PointIndex ptmp, bool check_only=false); void SplitImprove2 (); double SplitImprove2Element (ElementIndex ei, const Table & elements_of_point, bool check_only); @@ -46,7 +46,7 @@ public: void SwapImproveSurface (const BitArray * working_elements = NULL, const NgArray< idmap_type* > * idmaps = NULL); void SwapImprove2 (); - double SwapImprove2 (ElementIndex eli1, int face, Table & elementsonnode, TABLE & belementsonnode, bool check_only=false ); + double SwapImprove2 (ElementIndex eli1, int face, Table & elementsonnode, DynamicTable & belementsonnode, bool check_only=false ); void ImproveMesh() { mesh.ImproveMesh(mp, goal); } @@ -108,12 +108,12 @@ public: class PointFunction1 : public MinFunction { Mesh::T_POINTS & points; - const NgArray & faces; + const NgArray> & faces; const MeshingParameters & mp; double h; public: PointFunction1 (Mesh::T_POINTS & apoints, - const NgArray & afaces, + const NgArray> & afaces, const MeshingParameters & amp, double ah); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 02dccf0e..0ecaf4fa 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -1052,7 +1052,7 @@ namespace netgen // for (PointIndex pi = points.Begin(); pi < points.End(); pi++) for (PointIndex pi : points.Range()) if ((*this)[pi].Singularity()>=1.) - outfile << int(pi) << "\t" << (*this)[pi].Singularity() << endl; + outfile << pi << "\t" << (*this)[pi].Singularity() << endl; } cnt_sing = 0; @@ -1341,7 +1341,7 @@ namespace netgen el.SetNP(nep); el.SetCurved (nep != 4); for (int j = 0; j < nep; j++) - infile >> (int&)(el[j]); + infile >> el[j]; if (inverttets) el.Invert(); @@ -2113,7 +2113,7 @@ namespace netgen for (int j = 0; j < nep; j++) { - infile >> (int&)(el[j]); + infile >> el[j]; el[j] = el[j]+oldnp; } @@ -2182,7 +2182,7 @@ namespace netgen for (ElementIndex ei = 0; ei < volelements.Size(); ei++) { for (int j = 0; j < 4; j++) - if ( (*this)[ei][j] <= PointIndex::BASE-1) + if ( !(*this)[ei][j].IsValid()) { (*testout) << "El " << ei << " has 0 nodes: "; for (int k = 0; k < 4; k++) @@ -2223,9 +2223,9 @@ namespace netgen if (sel.GetNP() <= 4) for (int j = 0; j < sel.GetNP(); j++) { - INDEX_2 i2; - i2.I1() = sel.PNumMod(j+1); - i2.I2() = sel.PNumMod(j+2); + PointIndices<2> i2; + i2[0] = sel.PNumMod(j+1); + i2[1] = sel.PNumMod(j+2); i2.Sort(); boundaryedges->Set (i2, 1); } @@ -2233,9 +2233,9 @@ namespace netgen { for (int j = 0; j < 3; j++) { - INDEX_2 i2; - i2.I1() = sel[j]; - i2.I2() = sel[(j+1)%3]; + PointIndices<2> i2; + i2[0] = sel[j]; + i2[1] = sel[(j+1)%3]; i2.Sort(); boundaryedges->Set (i2, 1); } @@ -2298,7 +2298,7 @@ namespace netgen static Timer tn2se("Mesh::CalcSurfacesOfNode - surf on node"); static Timer tht("Mesh::CalcSurfacesOfNode - surfelementht"); // surfacesonnode.SetSize (GetNP()); - TABLE surfacesonnode(GetNP()); + DynamicTable surfacesonnode(GetNP()); // delete boundaryedges; // boundaryedges = NULL; @@ -2379,10 +2379,10 @@ namespace netgen const Element2d & sel = surfelements[sei]; if (sel.IsDeleted()) continue; - INDEX_3 i3; - i3.I1() = sel.PNum(1); - i3.I2() = sel.PNum(2); - i3.I3() = sel.PNum(3); + PointIndices<3> i3; + i3[0] = sel.PNum(1); + i3[1] = sel.PNum(2); + i3[2] = sel.PNum(3); i3.Sort(); surfelementht -> Set (i3, sei); // war das wichtig ??? sel.GetIndex()); } @@ -2479,7 +2479,7 @@ namespace netgen for (int i = 0; i < GetNSeg(); i++) { const Segment & seg = segments[i]; - INDEX_2 i2(seg[0], seg[1]); + PointIndices<2> i2(seg[0], seg[1]); i2.Sort(); //boundaryedges -> Set (i2, 2); @@ -2817,7 +2817,7 @@ namespace netgen hel.NormalizeNumbering(); if (hel.PNum(1) == pi) { - INDEX_3 i3(hel[0], hel[1], hel[2]); + PointIndices<3> i3(hel[0], hel[1], hel[2]); tval i2; i2.index = GetFaceDescriptor(ind).DomainIn(); i2.p4 = (hel.GetNP() == 3) @@ -2833,7 +2833,7 @@ namespace netgen hel.NormalizeNumbering(); if (hel.PNum(1) == pi) { - INDEX_3 i3(hel[0], hel[1], hel[2]); + PointIndices<3> i3(hel[0], hel[1], hel[2]); tval i2; i2.index = GetFaceDescriptor(ind).DomainOut(); i2.p4 = (hel.GetNP() == 3) @@ -2860,7 +2860,7 @@ namespace netgen if (hel[0] == pi) { - INDEX_3 i3(hel[0], hel[1], hel[2]); + PointIndices<3> i3(hel[0], hel[1], hel[2]); if (faceht.Used (i3)) { @@ -2889,7 +2889,7 @@ namespace netgen { hel.Invert(); hel.NormalizeNumbering(); - INDEX_3 i3(hel[0], hel[1], hel[2]); + PointIndices<3> i3(hel[0], hel[1], hel[2]); tval i2; i2.index = el.GetIndex(); diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 2ed0e481..7034f6ce 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -220,11 +220,10 @@ namespace netgen #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); @@ -237,6 +236,7 @@ namespace netgen friend bool operator== (PointIndex a, PointIndex b); friend bool operator!= (PointIndex a, PointIndex b); */ + public: constexpr PointIndex (t_invalid inv) : i(long(PointIndex::BASE)-1) { ; } // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; } @@ -260,6 +260,7 @@ namespace netgen void DoArchive (Archive & ar) { ar & i; } }; + /* inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); } inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); } @@ -367,12 +368,19 @@ namespace netgen { Sort(); } }; + constexpr inline size_t HashValue2 (const PointIndex & ind, size_t mask) + { return (ind-IndexBASE()) & mask; } } namespace ngcore { template constexpr inline T InvalidHash(); + + + template <> + constexpr inline netgen::PointIndex InvalidHash () + { return netgen::PointIndex::INVALID; } template <> constexpr inline netgen::PointIndices<2> InvalidHash> () diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 311b8315..8d3ef551 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -91,7 +91,7 @@ namespace netgen } PointFunction1 :: PointFunction1 (Mesh::T_POINTS & apoints, - const NgArray & afaces, + const NgArray> & afaces, const MeshingParameters & amp, double ah) : points(apoints), faces(afaces), mp(amp)