hashtables roundup to power of 2, optimize bisect

libsrc/general/hashtabl.cpp

@@ -213,19 +213,24 @@ namespace netgen
 
 
   BASE_INDEX_2_CLOSED_HASHTABLE ::
-  BASE_INDEX_2_CLOSED_HASHTABLE (int size)
+  BASE_INDEX_2_CLOSED_HASHTABLE (size_t size)
-    : hash(size)
+    : hash(RoundUp2(size))
   {
+    size = hash.Size();
+    mask = size-1;
     // hash.SetName ("i2-hashtable, hash");
 
     invalid = -1;
-    for (int i = 1; i <= size; i++)
+    for (size_t i = 0; i < size; i++)
-      hash.Elem(i).I1() = invalid;
+      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");

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<INDEX_2> 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[i] == ind) return i;
-	if (hash.Get(i).I1() == invalid) return 0;
+	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<T> &
 
 
 
-  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<T> :: PrintMemInfo (ostream & ost) const
 template<class T>
 inline INDEX_2_CLOSED_HASHTABLE<T> :: 
 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<T> ::
 Set (const INDEX_2 & ahash, const T & acont)
 {
   int pos;
-  PositionCreate (ahash, pos);
+  PositionCreate0 (ahash, pos);
-  hash.Elem(pos) = ahash;
+  hash[pos] = ahash;
-  cont.Elem(pos) = acont;
+  cont[pos] = acont;
 }
 
 template<class T>
 inline const T & INDEX_2_CLOSED_HASHTABLE<T> :: 
 Get (const INDEX_2 & ahash) const
 {
-  int pos = Position (ahash);
+  int pos = Position0 (ahash);
-  return cont.Get(pos);
+  return cont[pos];
 }
 
 template<class T>
 inline bool INDEX_2_CLOSED_HASHTABLE<T> :: 
 Used (const INDEX_2 & ahash) const
 {
-  int pos = Position (ahash);
+  int pos = Position0 (ahash);
-  return (pos != 0);
+  return (pos != -1);
 }
 
 template<class T>
 inline void INDEX_2_CLOSED_HASHTABLE<T> :: 
-SetData (int pos, const INDEX_2 & ahash, const T & acont)
+SetData0 (int pos, const INDEX_2 & ahash, const T & acont)
 {
-  hash.Elem(pos) = ahash;
+  hash[pos] = ahash;
-  cont.Elem(pos) = acont;
+  cont[pos] = acont;
 }
   
 template<class T>
 inline void INDEX_2_CLOSED_HASHTABLE<T> :: 
-GetData (int pos, INDEX_2 & ahash, T & acont) const
+GetData0 (int pos, INDEX_2 & ahash, T & acont) const
 {
-  ahash = hash.Get(pos);
+  ahash = hash[pos];
-  acont = cont.Get(pos);
+  acont = cont[pos];
 }
 
 template<class T>
 inline void INDEX_2_CLOSED_HASHTABLE<T> :: 
-SetData (int pos, const T & acont)
+SetData0 (int pos, const T & acont)
 {
-  cont.Elem(pos) = acont;
+  cont[pos] = acont;
 }
   
 template<class T>
 inline void INDEX_2_CLOSED_HASHTABLE<T> :: 
-GetData (int pos, T & acont) const
+GetData0 (int pos, T & acont) const
 {
-  acont = cont.Get(pos);
+  acont = cont[pos];
 }
 
 

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<PointIndex> & cutedges)
   {
-    int hanging = 0;
+    bool hanging = false;
-    for (int i = 1; i <= mtris.Size(); i++)
+    // 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],
+	      INDEX_2 edge(tri.pnums[j],
-			   mtris.Get(i).pnums[k]);
+			   tri.pnums[k]);
 	      edge.Sort();
 	      if (cutedges.Used (edge))
 		{
-		  mtris.Elem(i).marked = 1;
+		  tri.marked = 1;
-		  hanging = 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,6 +2765,7 @@ namespace netgen
     INDEX_2_CLOSED_HASHTABLE<PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
 
     bool noprojection = false;
+    NgProfiler::StopTimer (timer1a);
     
     for (int l = 1; l <= 1; l++)
       {
@@ -3169,7 +3179,6 @@ namespace netgen
 	//string yn;
 	//cin >> yn;
         
-        NgProfiler::StopTimer (timer1);
     
         (*testout) << "refine volume elements" << endl;
 	do
@@ -3370,6 +3379,9 @@ namespace netgen
 		  mesh.mlparentsurfaceelement.Append (i);
 		}
 
+            NgProfiler::StopTimer (timer_bisecttrig);
+            
+            
 	    int nquad = mquads.Size();
 	    for (int i = 1; i <= nquad; i++)
 	      if (mquads.Elem(i).marked)
@@ -3462,14 +3474,14 @@ namespace netgen
 		  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;
 
@@ -3580,6 +3596,8 @@ namespace netgen
 
     NgProfiler::StartTimer (timer2);
 
+    NgProfiler::StartTimer (timer2a);    
+    
     mtets.SetAllocSize (mtets.Size());
     mprisms.SetAllocSize (mprisms.Size());
     mids.SetAllocSize (mids.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.SetIndex (trig.surfid);
-	el.SetOrder (mtris.Get(i).order);
+	el.SetOrder (trig.order);
-	for (int j = 1; j <= 3; j++)
+	for (int j = 0; j < 3; j++)
 	  {
-	    el.PNum(j) = mtris.Get(i).pnums[j-1];
+	    el[j] = trig.pnums[j];
-	    el.GeomInfoPi(j) = mtris.Get(i).pgeominfo[j-1];
+	    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++)
+    for (int i = 0; i < cutedges.Size(); i++)
-      if (cutedges.UsedPos(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++)
+	for (int j = 0; j < cutedges.Size(); j++)
-	  if (cutedges.UsedPos(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();
@@ -3961,8 +3982,9 @@ namespace netgen
     NgProfiler::StopTimer (timer2);
     NgProfiler::StartTimer (timer3);
 
+    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");
 

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<int,PointIndex::BASE> surfacesonnode(GetNP());
 
     delete boundaryedges;
     boundaryedges = NULL;
 
     delete surfelementht;
+    surfelementht = nullptr;
     delete segmentht;
 
     /*
@@ -1606,9 +1608,11 @@ namespace netgen
       segmentht = new INDEX_2_HASHTABLE<int> (GetNSeg() + 1);
     */
 
+    if (dimension == 3)
       surfelementht = new INDEX_3_CLOSED_HASHTABLE<int> (3*GetNSE() + 1);
     segmentht = new INDEX_2_CLOSED_HASHTABLE<int> (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]));
+          int pos = boundaryedges -> Position0(INDEX_2::Sort(el[i], el[j]));
-          if (pos)
+          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:";
+    // ost << "surfs on node:";
-    surfacesonnode.PrintMemInfo (cout);
+    // surfacesonnode.PrintMemInfo (cout);
 
     ost << "boundaryedges: ";
     if (boundaryedges)

libsrc/meshing/meshclass.hpp

@@ -43,7 +43,7 @@ namespace netgen
 
 
     /// surface indices at boundary nodes
-    TABLE<int,PointIndex::BASE> surfacesonnode;
+    // TABLE<int,PointIndex::BASE> surfacesonnode;
     /// boundary edges  (1..normal bedge, 2..segment)
     INDEX_2_CLOSED_HASHTABLE<int> * boundaryedges;
     ///

libsrc/meshing/meshtype.hpp

@@ -422,6 +422,8 @@ namespace netgen
     { return FlatArray<const PointIndex> (np, &pnum[0]); }
     FlatArray<PointIndex> PNums ()
     { return FlatArray<PointIndex> (np, &pnum[0]); }
+    auto Vertices() const
+    { return FlatArray<const PointIndex> (GetNV(), &pnum[0]); }
     
     ///
     PointIndex & PNum (int i) { return pnum[i-1]; }