diff --git a/libsrc/core/table.hpp b/libsrc/core/table.hpp index 134ed8a9..4476a331 100644 --- a/libsrc/core/table.hpp +++ b/libsrc/core/table.hpp @@ -329,8 +329,8 @@ namespace ngcore case 1: { size_t oldval = nd; - while (blocknr+1>nd) { - nd.compare_exchange_weak (oldval, blocknr+1); + while (blocknr-IndexBASE()+1>nd) { + nd.compare_exchange_weak (oldval, blocknr-IndexBASE()+1); oldval = nd; } break; diff --git a/libsrc/general/template.hpp b/libsrc/general/template.hpp index 6f1571f1..f8fccddf 100644 --- a/libsrc/general/template.hpp +++ b/libsrc/general/template.hpp @@ -115,14 +115,19 @@ class INDEX_2 public: /// INDEX_2 () { } + INDEX_2 (const INDEX_2&) = default; + INDEX_2 (INDEX_2&&) = default; + + INDEX_2 & operator= (const INDEX_2&) = default; + INDEX_2 & operator= (INDEX_2&&) = default; /// constexpr INDEX_2 (INDEX ai1, INDEX ai2) : i{ai1, ai2} { } // { i[0] = ai1; i[1] = ai2; } /// - constexpr INDEX_2 (const INDEX_2 & in2) - : i{in2.i[0], in2.i[1]} { } + // constexpr INDEX_2 (const INDEX_2 & in2) + // : i{in2.i[0], in2.i[1]} { } // { i[0] = in2.i[0]; i[1] = in2.i[1]; } @@ -206,12 +211,14 @@ public: /// INDEX_3 () { } /// - INDEX_3 (INDEX ai1, INDEX ai2, INDEX ai3) - { i[0] = ai1; i[1] = ai2; i[2] = ai3; } + constexpr INDEX_3 (INDEX ai1, INDEX ai2, INDEX ai3) + : i{ai1, ai2, ai3} { } + // { i[0] = ai1; i[1] = ai2; i[2] = ai3; } /// - INDEX_3 (const INDEX_3 & in2) - { i[0] = in2.i[0]; i[1] = in2.i[1]; i[2] = in2.i[2]; } + constexpr INDEX_3 (const INDEX_3 & in2) + : i{in2.i[0], in2.i[1], in2.i[2]} { } + // { i[0] = in2.i[0]; i[1] = in2.i[1]; i[2] = in2.i[2]; } static INDEX_3 Sort (INDEX_3 i3) @@ -479,5 +486,12 @@ namespace netgen { return HashValue2(IVec<2,netgen::INDEX>(ind[0], ind[1]), mask); } + + constexpr inline size_t HashValue2 (const netgen::INDEX_3 & ind, size_t mask) + { + return HashValue2(IVec<3,netgen::INDEX>(ind[0], ind[1], ind[2]), mask); + } + + } #endif diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index d7e29ceb..e55dd657 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -698,7 +698,7 @@ namespace netgen // int i, j; int cnt = 0; - int found; + bool found; double len2, maxlen2; INDEX_2 ep; @@ -707,7 +707,7 @@ namespace netgen do { - found = 0; + found = false; maxlen2 = 1e30; // for (int i = 1; i <= mesh.GetNE(); i++) @@ -775,7 +775,7 @@ namespace netgen { maxlen2 = len2; ep = i2; - found = 1; + found = true; } } } @@ -833,8 +833,8 @@ namespace netgen e1.Sort(); e2.Sort(); - int used1 = edgenumber.Used (e1); - int used2 = edgenumber.Used (e2); + bool used1 = edgenumber.Used (e1); + bool used2 = edgenumber.Used (e2); if (used1 && !used2) { @@ -2906,7 +2906,8 @@ namespace netgen // INDEX_2_HASHTABLE cutedges(10 + 5 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); // INDEX_2_CLOSED_HASHTABLE cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); - ClosedHashTable cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); + // ClosedHashTable cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); + ClosedHashTable, PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); bool noprojection = false; NgProfiler::StopTimer (timer1a); diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index aa9a0536..3595f8f6 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -851,16 +851,17 @@ namespace netgen tets_with_3_bnd_points.SetSize(cnt); static Timer t1("Build face table"); t1.Start(); - ngcore::ClosedHashTable< ngcore::IVec<3>, int > face_table( 4*cnt + 3 ); - for(auto ei : tets_with_3_bnd_points) - for(auto j : Range(4)) + // ngcore::ClosedHashTable< ngcore::IVec<3>, int > face_table( 4*cnt + 3 ); + ngcore::ClosedHashTable< PointIndices<3>, int > face_table( 4*cnt + 3 ); + for (auto ei : tets_with_3_bnd_points) + for (auto j : Range(4)) { - PointIndices<3> i3_ = tempels[ei].GetFace (j); - ngcore::IVec<3> i3 = {i3_[0], i3_[1], i3_[2]}; - if(bnd_points[i3[0]] && bnd_points[i3[1]] && bnd_points[i3[2]]) + PointIndices<3> i3 = tempels[ei].GetFace (j); + // ngcore::IVec<3> i3 = {i3_[0], i3_[1], i3_[2]}; + if(bnd_points[i3[0]] && bnd_points[i3[1]] && bnd_points[i3[2]]) { - i3.Sort(); - face_table.Set( i3, true ); + i3.Sort(); + face_table.Set( i3, true ); } } t1.Stop(); @@ -871,7 +872,8 @@ namespace netgen for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & tri = mesh.OpenElement(i); - ngcore::IVec<3,PointIndex> i3(tri[0], tri[1], tri[2]); + // ngcore::IVec<3,PointIndex> i3(tri[0], tri[1], tri[2]); + PointIndices<3> i3(tri[0], tri[1], tri[2]); i3.Sort(); if(!face_table.Used(i3)) openels.Append(i); @@ -1016,7 +1018,7 @@ namespace netgen for (int j = 0; j < 4; j++) { pp[j] = &mesh.Point(el[j]); - tetpi[j] = el[j]; + tetpi[j] = el[j]-IndexBASE()+1; } Point3d tetpmin(*pp[0]); @@ -1045,7 +1047,7 @@ namespace netgen for (int k = 1; k <= 3; k++) { tripp[k-1] = &mesh.Point (tri.PNum(k)); - tripi[k-1] = tri.PNum(k); + tripi[k-1] = tri.PNum(k)-IndexBASE()+1; } if (IntersectTetTriangle (&pp[0], &tripp[0], tetpi, tripi)) @@ -1133,13 +1135,10 @@ namespace netgen for (auto & el : tempels) for (int j = 0; j < 4; j++) el.NB(j) = 0; - - TABLE elsonpoint(mesh.GetNP()); + /* - for (int i = 0; i < tempels.Size(); i++) - { - const DelaunayTet & el = tempels[i]; - */ + TABLE elsonpoint(mesh.GetNP()); + for (const DelaunayTet & el : tempels) { PointIndices<4> i4(el[0], el[1], el[2], el[3]); @@ -1158,7 +1157,25 @@ namespace netgen elsonpoint.Add (i4.I1(), i+1); elsonpoint.Add (i4.I2(), i+1); } + */ + TableCreator creator(mesh.GetNP()); + while (!creator.Done()) + { + for (int i = 0; i < tempels.Size(); i++) + { + const DelaunayTet & el = tempels[i]; + PointIndices<4> i4(el[0], el[1], el[2], el[3]); + i4.Sort(); + creator.Add (i4[0], i+1); + creator.Add (i4[1], i+1); + } + creator++; + } + auto elsonpoint = creator.MoveTable(); + + + // cout << "elsonpoint mem: "; // elsonpoint.PrintMemInfo(cout); diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp index e8ae2dbc..97b7a02c 100644 --- a/libsrc/meshing/improve2.cpp +++ b/libsrc/meshing/improve2.cpp @@ -713,7 +713,7 @@ namespace netgen if (mesh.IsSegment (pi1, pi2)) continue; - INDEX_2 ii2 (pi1, pi2); + PointIndices<2> ii2 (pi1, pi2); ii2.Sort(); if (els_on_edge.Used (ii2)) { @@ -739,7 +739,7 @@ namespace netgen if (mesh.LegalTrig(sel)) continue; // find longest edge - INDEX_2 edge; + PointIndices<2> edge; double edge_len = 0; PointIndex pi1, pi2, pi3, pi4; PointGeomInfo gi1, gi2, gi3, gi4; diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 4e036170..84e70af5 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -452,8 +452,10 @@ namespace netgen lock.Lock(); timestamp = NextTimeStamp(); - int maxn = max2 (s[0], s[1]); - maxn += 1-PointIndex::BASE; + // int maxn = max2 (s[0], s[1]); + // maxn += 1-PointIndex::BASE; + int maxn = max2 (s[0]-IndexBASE()+1, + s[1]-IndexBASE()+1); /* if (maxn > ptyps.Size()) @@ -540,12 +542,19 @@ namespace netgen void Mesh :: SetSurfaceElement (SurfaceElementIndex sei, const Element2d & el) { + /* int maxn = el[0]; for (int i = 1; i < el.GetNP(); i++) if (el[i] > maxn) maxn = el[i]; maxn += 1-PointIndex::BASE; + */ + PointIndex maxpi = el[0]; + for (int i = 1; i < el.GetNP(); i++) + if (el[i] > maxpi) maxpi = el[i]; + int maxn = maxpi-IndexBASE()+1; + if (maxn <= points.Size()) { for (int i = 0; i < el.GetNP(); i++) @@ -843,9 +852,12 @@ namespace netgen outfile.setf (ios::fixed, ios::floatfield); outfile.setf (ios::showpoint); + /* PointIndex pi; for (pi = PointIndex::BASE; pi < GetNP()+PointIndex::BASE; pi++) + */ + for (PointIndex pi : (*this).Points().Range()) { outfile.width(22); outfile << (*this)[pi](0)/scale << " "; @@ -1743,7 +1755,7 @@ namespace netgen } maxglob = comm.AllReduce (maxglob, NG_MPI_MAX); - int numglob = maxglob+1-PointIndex::BASE; + int numglob = maxglob+1-IndexBASE(); if (comm.Rank() > 0) { comm.Send (globnum, 0, 200); diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index f9078e4b..5ba984ba 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -299,7 +299,7 @@ namespace netgen - Element2d :: Element2d (int pi1, int pi2, int pi3) + Element2d :: Element2d (PointIndex pi1, PointIndex pi2, PointIndex pi3) { pnum[0] = pi1; pnum[1] = pi2; @@ -322,7 +322,7 @@ namespace netgen is_curved = false; } - Element2d :: Element2d (int pi1, int pi2, int pi3, int pi4) + Element2d :: Element2d (PointIndex pi1, PointIndex pi2, PointIndex pi3, PointIndex pi4) { pnum[0] = pi1; pnum[1] = pi2; diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 017432a5..9c0d1415 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -297,8 +297,13 @@ namespace netgen { public: PointIndices () = default; - PointIndices (INDEX_2 i2) : INDEX_2(i2) { ; } - PointIndices (PointIndex i1, PointIndex i2) : INDEX_2(i1,i2) { ; } + PointIndices (const PointIndices&) = default; + PointIndices (PointIndices&&) = default; + PointIndices & operator= (const PointIndices&) = default; + PointIndices & operator= (PointIndices&&) = default; + + constexpr PointIndices (INDEX_2 i2) : INDEX_2(i2) { ; } + constexpr PointIndices (PointIndex i1, PointIndex i2) : INDEX_2(i1,i2) { ; } PointIndex operator[] (int i) const { return PointIndex(INDEX_2::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast(INDEX_2::operator[](i)); } using INDEX_2::Sort; @@ -310,8 +315,12 @@ namespace netgen { public: PointIndices () = default; - PointIndices (INDEX_3 i3) : INDEX_3(i3) { ; } - PointIndices (PointIndex i1, PointIndex i2, PointIndex i3) : INDEX_3(i1,i2,i3) { ; } + PointIndices (const PointIndices&) = default; + PointIndices (PointIndices&&) = default; + PointIndices & operator= (const PointIndices&) = default; + PointIndices & operator= (PointIndices&&) = default; + constexpr PointIndices (INDEX_3 i3) : INDEX_3(i3) { ; } + constexpr PointIndices (PointIndex i1, PointIndex i2, PointIndex i3) : INDEX_3(i1,i2,i3) { ; } PointIndex operator[] (int i) const { return PointIndex(INDEX_3::operator[](i)); } PointIndex & operator[] (int i) { return reinterpret_cast(INDEX_3::operator[](i)); } using INDEX_3::Sort; @@ -336,6 +345,20 @@ namespace netgen } +namespace ngcore +{ + template constexpr inline T InvalidHash(); + + template <> + constexpr inline netgen::PointIndices<2> InvalidHash> () + { return netgen::PointIndices<2>{netgen::PointIndex::INVALID, netgen::PointIndex::INVALID}; } + + template <> + constexpr inline netgen::PointIndices<3> InvalidHash> () + { return netgen::PointIndices<3>{netgen::PointIndex::INVALID, netgen::PointIndex::INVALID, netgen::PointIndex::INVALID}; } +} + + namespace std { // structured binding support @@ -584,9 +607,9 @@ namespace netgen /// DLL_HEADER Element2d (ELEMENT_TYPE type); /// - DLL_HEADER Element2d (int pi1, int pi2, int pi3); + DLL_HEADER Element2d (PointIndex pi1, PointIndex pi2, PointIndex pi3); /// - DLL_HEADER Element2d (int pi1, int pi2, int pi3, int pi4); + DLL_HEADER Element2d (PointIndex pi1, PointIndex pi2, PointIndex pi3, PointIndex pi4); /// ELEMENT_TYPE GetType () const { return typ; } ///