diff --git a/libsrc/general/hashtabl.cpp b/libsrc/general/hashtabl.cpp index fd2fbfbd..9856431e 100644 --- a/libsrc/general/hashtabl.cpp +++ b/libsrc/general/hashtabl.cpp @@ -213,19 +213,24 @@ namespace netgen BASE_INDEX_2_CLOSED_HASHTABLE :: - BASE_INDEX_2_CLOSED_HASHTABLE (int size) - : hash(size) + BASE_INDEX_2_CLOSED_HASHTABLE (size_t size) + : hash(RoundUp2(size)) { + size = hash.Size(); + mask = size-1; // hash.SetName ("i2-hashtable, hash"); invalid = -1; - for (int i = 1; i <= size; i++) - hash.Elem(i).I1() = invalid; + for (size_t i = 0; i < size; i++) + hash[i].I1() = invalid; } void BASE_INDEX_2_CLOSED_HASHTABLE :: BaseSetSize (int size) { + size = RoundUp2 (size); + mask = size-1; + hash.SetSize(size); for (int i = 1; i <= size; i++) hash.Elem(i).I1() = invalid; @@ -245,25 +250,28 @@ namespace netgen } } - int BASE_INDEX_2_CLOSED_HASHTABLE :: + bool BASE_INDEX_2_CLOSED_HASHTABLE :: PositionCreate2 (const INDEX_2 & ind, int & apos) { int i = HashValue(ind); int startpos = i; while (1) { + /* i++; if (i > hash.Size()) i = 1; - if (hash.Get(i) == ind) + */ + i = (i+1) % hash.Size(); + if (hash[i] == ind) { apos = i; - return 0; + return false; } - if (hash.Get(i).I1() == invalid) + if (hash[i].I1() == invalid) { - hash.Elem(i) = ind; + hash[i] = ind; apos = i; - return 1; + return true; } if (i == startpos) throw NgException ("Try to set new element in full closed hashtable"); diff --git a/libsrc/general/hashtabl.hpp b/libsrc/general/hashtabl.hpp index b856bbc6..72263cba 100644 --- a/libsrc/general/hashtabl.hpp +++ b/libsrc/general/hashtabl.hpp @@ -628,7 +628,12 @@ public: - + inline size_t RoundUp2 (size_t i) + { + size_t res = 1; + while (res < i) res *= 2; // hope it will never be too large + return res; + } /// Closed Hashing HT @@ -640,47 +645,52 @@ protected: Array hash; /// int invalid; + size_t mask; public: /// - BASE_INDEX_2_CLOSED_HASHTABLE (int size); + BASE_INDEX_2_CLOSED_HASHTABLE (size_t size); int Size() const { return hash.Size(); } - int UsedPos (int pos) const { return ! (hash.Get(pos).I1() == invalid); } + bool UsedPos0 (int pos) const { return ! (hash[pos].I1() == invalid); } int UsedElements () const; /// int HashValue (const INDEX_2 & ind) const { - return (ind.I1() + 71 * ind.I2()) % hash.Size() + 1; + // return (ind.I1() + 71 * ind.I2()) % hash.Size() + 1; + return (ind.I1() + 71 * ind.I2()) & mask; } - int Position (const INDEX_2 & ind) const + int Position0 (const INDEX_2 & ind) const { int i = HashValue(ind); while (1) { - if (hash.Get(i) == ind) return i; - if (hash.Get(i).I1() == invalid) return 0; + if (hash[i] == ind) return i; + if (hash[i].I1() == invalid) return -1; + i = (i+1) & mask; + /* i++; if (i > hash.Size()) i = 1; + */ } } // returns 1, if new postion is created - int PositionCreate (const INDEX_2 & ind, int & apos) + bool PositionCreate0 (const INDEX_2 & ind, int & apos) { int i = HashValue (ind); - if (hash.Get(i) == ind) + if (hash[i] == ind) { apos = i; - return 0; + return false; } - if (hash.Get(i).I1() == invalid) + if (hash[i].I1() == invalid) { - hash.Elem(i) = ind; + hash[i] = ind; apos = i; - return 1; + return true; } return PositionCreate2 (ind, apos); } @@ -689,7 +699,7 @@ protected: /// int Position2 (const INDEX_2 & ind) const; - int PositionCreate2 (const INDEX_2 & ind, int & apos); + bool PositionCreate2 (const INDEX_2 & ind, int & apos); void BaseSetSize (int asize); }; @@ -711,15 +721,15 @@ public: /// inline bool Used (const INDEX_2 & ahash) const; /// - inline void SetData (int pos, const INDEX_2 & ahash, const T & acont); + inline void SetData0 (int pos, const INDEX_2 & ahash, const T & acont); /// - inline void GetData (int pos, INDEX_2 & ahash, T & acont) const; + inline void GetData0 (int pos, INDEX_2 & ahash, T & acont) const; /// - inline void SetData (int pos, const T & acont); + inline void SetData0 (int pos, const T & acont); /// - inline void GetData (int pos, T & acont) const; + inline void GetData0 (int pos, T & acont) const; /// - const T & GetData (int pos) { return cont.Get(pos); } + const T & GetData0 (int pos) { return cont[pos]; } /// inline void SetSize (int size); /// @@ -755,13 +765,6 @@ inline ostream & operator<< (ostream & ost, const INDEX_2_CLOSED_HASHTABLE & - inline size_t RoundUp2 (size_t i) - { - size_t res = 1; - while (res < i) res *= 2; // hope it will never be too large - return res; - } - class BASE_INDEX_3_CLOSED_HASHTABLE { protected: @@ -1192,7 +1195,7 @@ inline void INDEX_HASHTABLE :: PrintMemInfo (ostream & ost) const template inline INDEX_2_CLOSED_HASHTABLE :: INDEX_2_CLOSED_HASHTABLE (int size) - : BASE_INDEX_2_CLOSED_HASHTABLE(size), cont(size) + : BASE_INDEX_2_CLOSED_HASHTABLE(size), cont(RoundUp2(size)) { // cont.SetName ("i2-hashtable, contents"); } @@ -1202,55 +1205,55 @@ inline void INDEX_2_CLOSED_HASHTABLE :: Set (const INDEX_2 & ahash, const T & acont) { int pos; - PositionCreate (ahash, pos); - hash.Elem(pos) = ahash; - cont.Elem(pos) = acont; + PositionCreate0 (ahash, pos); + hash[pos] = ahash; + cont[pos] = acont; } template inline const T & INDEX_2_CLOSED_HASHTABLE :: Get (const INDEX_2 & ahash) const { - int pos = Position (ahash); - return cont.Get(pos); + int pos = Position0 (ahash); + return cont[pos]; } template inline bool INDEX_2_CLOSED_HASHTABLE :: Used (const INDEX_2 & ahash) const { - int pos = Position (ahash); - return (pos != 0); + int pos = Position0 (ahash); + return (pos != -1); } template inline void INDEX_2_CLOSED_HASHTABLE :: -SetData (int pos, const INDEX_2 & ahash, const T & acont) +SetData0 (int pos, const INDEX_2 & ahash, const T & acont) { - hash.Elem(pos) = ahash; - cont.Elem(pos) = acont; + hash[pos] = ahash; + cont[pos] = acont; } template inline void INDEX_2_CLOSED_HASHTABLE :: -GetData (int pos, INDEX_2 & ahash, T & acont) const +GetData0 (int pos, INDEX_2 & ahash, T & acont) const { - ahash = hash.Get(pos); - acont = cont.Get(pos); + ahash = hash[pos]; + acont = cont[pos]; } template inline void INDEX_2_CLOSED_HASHTABLE :: -SetData (int pos, const T & acont) +SetData0 (int pos, const T & acont) { - cont.Elem(pos) = acont; + cont[pos] = acont; } template inline void INDEX_2_CLOSED_HASHTABLE :: -GetData (int pos, T & acont) const +GetData0 (int pos, T & acont) const { - acont = cont.Get(pos); + acont = cont[pos]; } diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 1662417b..d4941b97 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -1723,27 +1723,28 @@ namespace netgen - int MarkHangingTris (T_MTRIS & mtris, + bool MarkHangingTris (T_MTRIS & mtris, const INDEX_2_CLOSED_HASHTABLE & cutedges) { - int hanging = 0; - for (int i = 1; i <= mtris.Size(); i++) + bool hanging = false; + // for (int i = 1; i <= mtris.Size(); i++) + for (auto & tri : mtris) { - if (mtris.Get(i).marked) + if (tri.marked) { - hanging = 1; + hanging = true; continue; } for (int j = 0; j < 2; j++) for (int k = j+1; k < 3; k++) { - INDEX_2 edge(mtris.Get(i).pnums[j], - mtris.Get(i).pnums[k]); + INDEX_2 edge(tri.pnums[j], + tri.pnums[k]); edge.Sort(); if (cutedges.Used (edge)) { - mtris.Elem(i).marked = 1; - hanging = 1; + tri.marked = 1; + hanging = true; } } } @@ -2620,14 +2621,22 @@ namespace netgen static int timer = NgProfiler::CreateTimer ("Bisect"); static int timer1 = NgProfiler::CreateTimer ("Bisect 1"); + static int timer1a = NgProfiler::CreateTimer ("Bisect 1a"); + static int timer1b = NgProfiler::CreateTimer ("Bisect 1b"); static int timer2 = NgProfiler::CreateTimer ("Bisect 2"); + static int timer2a = NgProfiler::CreateTimer ("Bisect 2a"); + static int timer2b = NgProfiler::CreateTimer ("Bisect 2b"); + static int timer2c = NgProfiler::CreateTimer ("Bisect 2c"); static int timer3 = NgProfiler::CreateTimer ("Bisect 3"); + static int timer3a = NgProfiler::CreateTimer ("Bisect 3a"); + static int timer3b = NgProfiler::CreateTimer ("Bisect 3b"); static int timer_bisecttet = NgProfiler::CreateTimer ("Bisect tets"); static int timer_bisecttrig = NgProfiler::CreateTimer ("Bisect trigs"); + static int timer_bisectsegms = NgProfiler::CreateTimer ("Bisect segms"); NgProfiler::RegionTimer reg1 (timer); NgProfiler::StartTimer (timer1); - + NgProfiler::StartTimer (timer1a); static int localizetimer = NgProfiler::CreateTimer("localize edgepoints"); NgProfiler::RegionTimer * loct = new NgProfiler::RegionTimer(localizetimer); LocalizeEdgePoints(mesh); @@ -2756,7 +2765,8 @@ namespace netgen INDEX_2_CLOSED_HASHTABLE cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size())); bool noprojection = false; - + NgProfiler::StopTimer (timer1a); + for (int l = 1; l <= 1; l++) { int marked = 0; @@ -3169,7 +3179,6 @@ namespace netgen //string yn; //cin >> yn; - NgProfiler::StopTimer (timer1); (*testout) << "refine volume elements" << endl; do @@ -3369,7 +3378,10 @@ namespace netgen mtris.Append (newtri2); mesh.mlparentsurfaceelement.Append (i); } - + + NgProfiler::StopTimer (timer_bisecttrig); + + int nquad = mquads.Size(); for (int i = 1; i <= nquad; i++) if (mquads.Elem(i).marked) @@ -3461,15 +3473,15 @@ namespace netgen mquads.Elem(i) = newquad1; mquads.Append (newquad2); } - - NgProfiler::StopTimer (timer_bisecttrig); - + + NgProfiler::StartTimer (timer1b); hangingsurf = MarkHangingTris (mtris, cutedges) + MarkHangingQuads (mquads, cutedges); hangingedge = 0; - + NgProfiler::StopTimer (timer1b); + NgProfiler::StartTimer (timer_bisectsegms); int nseg = mesh.GetNSeg (); for (int i = 1; i <= nseg; i++) { @@ -3506,6 +3518,8 @@ namespace netgen mesh.AddSegment (nseg2); } } + + NgProfiler::StopTimer (timer_bisectsegms); } while (hangingvol || hangingsurf || hangingedge); @@ -3534,6 +3548,8 @@ namespace netgen PrintMessage (4, mesh.GetNP(), " points"); } + NgProfiler::StopTimer (timer1); + // (*testout) << "mtets = " << mtets << endl; @@ -3579,6 +3595,8 @@ namespace netgen } NgProfiler::StartTimer (timer2); + + NgProfiler::StartTimer (timer2a); mtets.SetAllocSize (mtets.Size()); mprisms.SetAllocSize (mprisms.Size()); @@ -3654,15 +3672,18 @@ namespace netgen } mesh.ClearSurfaceElements(); - for (int i = 1; i <= mtris.Size(); i++) + mesh.SurfaceElements().SetAllocSize (mtris.Size()+mquads.Size()); + NgProfiler::StopTimer (timer2a); + NgProfiler::StartTimer (timer2b); + for (auto & trig : mtris) { Element2d el(TRIG); - el.SetIndex (mtris.Get(i).surfid); - el.SetOrder (mtris.Get(i).order); - for (int j = 1; j <= 3; j++) + el.SetIndex (trig.surfid); + el.SetOrder (trig.order); + for (int j = 0; j < 3; j++) { - el.PNum(j) = mtris.Get(i).pnums[j-1]; - el.GeomInfoPi(j) = mtris.Get(i).pgeominfo[j-1]; + el[j] = trig.pnums[j]; + el.GeomInfoPi(j+1) = trig.pgeominfo[j]; } mesh.AddSurfaceElement (el); } @@ -3675,7 +3696,7 @@ namespace netgen Swap (el.PNum(3), el.PNum(4)); mesh.AddSurfaceElement (el); } - + NgProfiler::StopTimer (timer2b); // write multilevel hierarchy to mesh: @@ -3705,12 +3726,12 @@ namespace netgen BitArray isnewpoint(np); isnewpoint.Clear(); - for (int i = 1; i <= cutedges.Size(); i++) - if (cutedges.UsedPos(i)) + for (int i = 0; i < cutedges.Size(); i++) + if (cutedges.UsedPos0(i)) { INDEX_2 edge; PointIndex newpi; - cutedges.GetData (i, edge, newpi); + cutedges.GetData0 (i, edge, newpi); isnewpoint.Set(newpi); mesh.mlbetweennodes.Elem(newpi) = edge; } @@ -3811,12 +3832,12 @@ namespace netgen } */ - for (int j = 1; j <= cutedges.Size(); j++) - if (cutedges.UsedPos(j)) + for (int j = 0; j < cutedges.Size(); j++) + if (cutedges.UsedPos0(j)) { INDEX_2 i2; PointIndex newpi; - cutedges.GetData (j, i2, newpi); + cutedges.GetData0 (j, i2, newpi); INDEX_2 oi2(identmap.Get(i2.I1()), identmap.Get(i2.I2())); oi2.Sort(); @@ -3960,9 +3981,10 @@ namespace netgen NgProfiler::StopTimer (timer2); NgProfiler::StartTimer (timer3); - - mesh.UpdateTopology(opt.task_manager); + NgProfiler::StartTimer (timer3a); + mesh.UpdateTopology(opt.task_manager); + NgProfiler::StopTimer (timer3a); @@ -3976,8 +3998,9 @@ namespace netgen ofst.close(); } - + NgProfiler::StartTimer (timer3b); mesh.CalcSurfacesOfNode(); + NgProfiler::StopTimer (timer3b); PrintMessage (1, "Bisection done"); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 396bb167..17fd6ba4 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -116,7 +116,7 @@ namespace netgen surfelements.SetSize(0); volelements.SetSize(0); lockedpoints.SetSize(0); - surfacesonnode.SetSize(0); + // surfacesonnode.SetSize(0); delete boundaryedges; boundaryedges = NULL; @@ -1593,12 +1593,14 @@ namespace netgen void Mesh :: CalcSurfacesOfNode () { - surfacesonnode.SetSize (GetNP()); + // surfacesonnode.SetSize (GetNP()); + TABLE surfacesonnode(GetNP()); delete boundaryedges; boundaryedges = NULL; delete surfelementht; + surfelementht = nullptr; delete segmentht; /* @@ -1606,9 +1608,11 @@ namespace netgen segmentht = new INDEX_2_HASHTABLE (GetNSeg() + 1); */ - surfelementht = new INDEX_3_CLOSED_HASHTABLE (3*GetNSE() + 1); + if (dimension == 3) + surfelementht = new INDEX_3_CLOSED_HASHTABLE (3*GetNSE() + 1); segmentht = new INDEX_2_CLOSED_HASHTABLE (3*GetNSeg() + 1); + if (dimension == 3) for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { const Element2d & sel = surfelements[sei]; @@ -1619,6 +1623,9 @@ namespace netgen for (int j = 0; j < sel.GetNP(); j++) { PointIndex pi = sel[j]; + if (!surfacesonnode[pi].Contains(si)) + surfacesonnode.Add (pi, si); + /* bool found = 0; for (int k = 0; k < surfacesonnode[pi].Size(); k++) if (surfacesonnode[pi][k] == si) @@ -1629,6 +1636,7 @@ namespace netgen if (!found) surfacesonnode.Add (pi, si); + */ } } /* @@ -1647,6 +1655,8 @@ namespace netgen surfelementht -> AllocateElements(); */ + + if (dimension==3) for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { const Element2d & sel = surfelements[sei]; @@ -3665,11 +3675,11 @@ namespace netgen { bool sege = false, be = false; - int pos = boundaryedges -> Position(INDEX_2::Sort(el[i], el[j])); - if (pos) + int pos = boundaryedges -> Position0(INDEX_2::Sort(el[i], el[j])); + if (pos != -1) { be = true; - if (boundaryedges -> GetData(pos) == 2) + if (boundaryedges -> GetData0(pos) == 2) sege = true; } @@ -5545,6 +5555,7 @@ namespace netgen if (el[j] > numvertices) numvertices = el[j]; } + /* for (i = 1; i <= nse; i++) { const Element2d & el = SurfaceElement(i); @@ -5553,6 +5564,10 @@ namespace netgen if (el.PNum(j) > numvertices) numvertices = el.PNum(j); } + */ + for (auto & el : SurfaceElements()) + for (PointIndex v : el.Vertices()) + if (v > numvertices) numvertices = v; numvertices += 1- PointIndex::BASE; } @@ -5917,8 +5932,8 @@ namespace netgen << sizeof (Element) << " = " << GetNE() * sizeof(Element) << endl; - ost << "surfs on node:"; - surfacesonnode.PrintMemInfo (cout); + // ost << "surfs on node:"; + // surfacesonnode.PrintMemInfo (cout); ost << "boundaryedges: "; if (boundaryedges) diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index cf5e7de2..6ec7627a 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -43,7 +43,7 @@ namespace netgen /// surface indices at boundary nodes - TABLE surfacesonnode; + // TABLE surfacesonnode; /// boundary edges (1..normal bedge, 2..segment) INDEX_2_CLOSED_HASHTABLE * boundaryedges; /// diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 92acaaf7..1471c541 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -422,6 +422,8 @@ namespace netgen { return FlatArray (np, &pnum[0]); } FlatArray PNums () { return FlatArray (np, &pnum[0]); } + auto Vertices() const + { return FlatArray (GetNV(), &pnum[0]); } /// PointIndex & PNum (int i) { return pnum[i-1]; }