tuning of Delaunay2d: FindInnerPoints, use of edgeHT

libsrc/meshing/adfront2.cpp

@@ -503,7 +503,7 @@ namespace netgen
   }
 
   bool AdFront2 :: SameSide (const Point<2> & lp1, const Point<2> & lp2, 
-                             const NgArray<int> * testfaces) const
+                             const FlatArray<int> * testfaces) const
   {
     int cnt = 0;
 

libsrc/meshing/adfront2.hpp

@@ -262,7 +262,7 @@ public:
   bool Inside (const Point<2> & p) const;
 
   bool SameSide (const Point<2> & lp1, const Point<2> & lp2, 
-                 const NgArray<int> * /* testfaces */ = NULL) const;
+                 const FlatArray<int> * /* testfaces */ = NULL) const;
   /*
   {
     return Inside (lp1) == Inside (lp2);

libsrc/meshing/delaunay2d.cpp

@@ -89,12 +89,18 @@ namespace netgen
         Swap(p0,p1);
 
       INT<2> hash = {p0,p1};
+      /*
       if(!edge_to_trig.Used(hash))
         return -1;
 
       auto i2 = edge_to_trig.Get({p0,p1});
 
       return i2[0] == eli ? i2[1] : i2[0];
+      */
+      auto pos = edge_to_trig.Position(hash);
+      if (pos == -1) return -1;
+      auto i2 = edge_to_trig.GetData(pos);
+      return i2[0] == eli ? i2[1] : i2[0];
     }
 
     void SetNeighbour( int eli, int edge )
@@ -105,6 +111,22 @@ namespace netgen
         Swap(p0,p1);
 
       INT<2> hash = {p0,p1};
+      auto pos = edge_to_trig.Position(hash);
+      if (pos == -1)
+        edge_to_trig[hash] = {eli, -1};
+      else
+        {
+          auto i2 = edge_to_trig.GetData(pos);
+          if(i2[0]==-1)
+            i2[0] = eli;
+          else
+            {
+              if(i2[1]==-1)
+                i2[1] = eli;
+            }
+          edge_to_trig.SetData (pos, i2);
+        }
+      /*
       if(!edge_to_trig.Used(hash))
         edge_to_trig[hash] = {eli, -1};
       else
@@ -122,6 +144,7 @@ namespace netgen
 
         edge_to_trig[hash] = i2;
       }
+      */
     }
 
     void UnsetNeighbours( int eli )
@@ -477,18 +500,31 @@ namespace netgen
         *testout << "npoints = " << endl << npoints << endl;
       }
 
-    for (int i = 1; i <= npoints.Size(); i++)
+
+    int prims[] = { 211, 223, 227, 229, 233, 239, 241, 251, 257, 263 }; 
+    int prim;
+  
     {
-	if (meshbox.IsIn (npoints.Get(i)))
+      int i = 0;
+      while (npoints.Size() % prims[i] == 0) i++;
+      prim = prims[i];
+    }
+
+    
+    for (int i = 0; i < npoints.Size(); i++)
       {
-	    PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
+        size_t hi = (size_t(prim) * size_t(i)) % npoints.Size();
-	    adfront.AddPoint (npoints.Get(i), gpnum);
+        
+	if (meshbox.IsIn (npoints[hi]))
+	  {
+	    PointIndex gpnum = mesh.AddPoint (npoints[hi]);
+	    adfront.AddPoint (npoints[hi], gpnum);
 	    
 	    if (debugparam.slowchecks)
 	      {
-		(*testout) << npoints.Get(i) << endl;
+		(*testout) << npoints[hi] << endl;
 
-		Point<2> p2d (npoints.Get(i)(0), npoints.Get(i)(1));
+		Point<2> p2d (npoints[hi](0), npoints[hi](1));
 		if (!adfront.Inside(p2d))
 		  {
 		    cout << "add outside point" << endl;

libsrc/meshing/localh.cpp

@@ -569,66 +569,35 @@ namespace netgen
 
 
     int nf = adfront->GetNFL();
-    NgArray<int> faceinds(nf);
+    Array<int> faceinds(nf);
-    NgArray<Box<3> > faceboxes(nf);
+    Array<Box<2>> faceboxes(nf);
 
     for (int i = 0; i < nf; i++)
       {
 	faceinds[i] = i;
-	// adfront->GetFaceBoundingBox(i, faceboxes.Elem(i));
-
 	const FrontLine & line = adfront->GetLine(i);
-	faceboxes[i].Set (adfront->GetPoint (line.L().I1()));
+        Point<3> p1 = adfront->GetPoint (line.L().I1());
-	faceboxes[i].Add (adfront->GetPoint (line.L().I2()));
+        Point<3> p2 = adfront->GetPoint (line.L().I2());
+        
+	faceboxes[i].Set (Point<2> (p1(0), p1(1)));
+	faceboxes[i].Add (Point<2> (p2(0), p2(1)));
       }
   
     for (int i = 0; i < 8; i++)
-      FindInnerBoxesRec2 (root->childs[i], adfront, faceboxes, faceinds, nf);
+      FindInnerBoxesRec2 (root->childs[i], adfront, faceboxes, faceinds); // , nf);
   }
 
 
   void LocalH :: 
   FindInnerBoxesRec2 (GradingBox * box,
 		      class AdFront2 * adfront, 
-		      NgArray<Box<3> > & faceboxes,
+		      FlatArray<Box<2>> faceboxes,
-		      NgArray<int> & faceinds, int nfinbox)
+		      FlatArray<int> faceinds) // , int nfinbox)
   {
     if (!box) return;
 
     GradingBox * father = box -> father;    
     
-    Point3d c(box->xmid[0], box->xmid[1], 0); // box->xmid[2]);
-    Vec3d v(box->h2, box->h2, box->h2);
-    Box3d boxc(c-v, c+v);
-
-    Point3d fc(father->xmid[0], father->xmid[1], 0); // father->xmid[2]);
-    Vec3d fv(father->h2, father->h2, father->h2);
-    Box3d fboxc(fc-fv, fc+fv);
-    Box3d boxcfc(c,fc);
-
-    NgArrayMem<int, 100> faceused;
-    NgArrayMem<int, 100> faceused2;
-    NgArrayMem<int, 100> facenotused;
-
-    for (int j = 0; j < nfinbox; j++)
-      {
-	//      adfront->GetFaceBoundingBox (faceinds.Get(j), facebox);
-	const Box3d & facebox = faceboxes[faceinds[j]];
-  
-	if (boxc.Intersect (facebox))
-	  faceused.Append(faceinds[j]);
-	else
-	  facenotused.Append(faceinds[j]);
-
-	if (boxcfc.Intersect (facebox))
-	  faceused2.Append (faceinds[j]);
-      }
-  
-    for (int j = 0; j < faceused.Size(); j++)
-      faceinds[j] = faceused[j];
-    for (int j = 0; j < facenotused.Size(); j++)
-      faceinds[j+faceused.Size()] = facenotused[j];
-  
     if (!father->flags.cutboundary)
       {
 	box->flags.isinner = father->flags.isinner;
@@ -636,17 +605,36 @@ namespace netgen
       }
     else
       {        
-	Point3d cf(father->xmid[0], father->xmid[1], father->xmid[2]);
-      
 	if (father->flags.isinner)
           {
             box->flags.pinner = 1;
           }
 	else
 	  {
-	    Point<2> c2d (c.X(), c.Y());
+            Point<2> c(box->xmid[0], box->xmid[1]); 
-	    Point<2> cf2d (cf.X(), cf.Y());
+            Point<2> fc(father->xmid[0], father->xmid[1]); 
-            bool sameside = adfront->SameSide (c2d, cf2d, &faceused2);
+            Box<2> boxcfc(c,fc);
+
+            // reorder: put faces cutting boxcfc first:
+            int iused = 0;
+            int inotused = faceinds.Size()-1;
+            while (iused <= inotused)
+              {
+                while ( (iused <= inotused) && boxcfc.Intersect (faceboxes[faceinds[iused]]))
+                  iused++;
+                while ( (iused <= inotused) && !boxcfc.Intersect (faceboxes[faceinds[inotused]]))
+                  inotused--;
+                if (iused < inotused)
+                  {
+                    Swap (faceinds[iused], faceinds[inotused]);
+                    iused++;
+                    inotused--;
+                  }
+              }
+
+            // bool sameside = adfront->SameSide (c2d, cf2d, &faceused2);
+            auto sub = faceinds.Range(0, iused);
+            bool sameside = adfront->SameSide (c, fc, &sub);
             if (sameside)
 	      box->flags.pinner = father->flags.pinner;
 	    else
@@ -659,11 +647,34 @@ namespace netgen
 	  box->flags.isinner = box->flags.pinner;
       }
 
-    // cout << "faceused: " << faceused.Size() << ", " << faceused2.Size() << ", " << facenotused.Size() << endl;
 
-    int nf = faceused.Size();
+    
+    int iused = 0;
+    if (faceinds.Size())
+      {
+        Point<2> c(box->xmid[0], box->xmid[1]); // box->xmid[2]);
+        Vec<2> v(box->h2, box->h2);
+        Box<2> boxc(c-v, c+v);
+        
+        // reorder again: put faces cutting boxc first:    
+        int inotused = faceinds.Size()-1;
+        while (iused <= inotused)
+          {
+            while ( (iused <= inotused) && boxc.Intersect (faceboxes[faceinds[iused]]))
+              iused++;
+            while ( (iused <= inotused) && !boxc.Intersect (faceboxes[faceinds[inotused]]))
+              inotused--;
+            if (iused < inotused)
+              {
+                Swap (faceinds[iused], faceinds[inotused]);
+                iused++;
+                inotused--;
+              }
+          }
+      }
+    
     for (int i = 0; i < 8; i++)
-      FindInnerBoxesRec2 (box->childs[i], adfront, faceboxes, faceinds, nf);
+      FindInnerBoxesRec2 (box->childs[i], adfront, faceboxes, faceinds.Range(0,iused));
   }
 
 

libsrc/meshing/localh.hpp

@@ -30,10 +30,16 @@ namespace netgen
 
     struct 
     {
+      /*
       unsigned int cutboundary:1;
       unsigned int isinner:1;
       unsigned int oldcell:1;
       unsigned int pinner:1;
+      */
+      bool cutboundary;
+      bool isinner;
+      bool oldcell;
+      bool pinner;
     } flags;
 
     ///
@@ -158,8 +164,8 @@ namespace netgen
     ///
     void FindInnerBoxesRec2 (GradingBox * box,
 			     class AdFront2 * adfront,
-			     NgArray<Box<3> > & faceboxes,
+			     FlatArray<Box<2>> faceboxes,
-			     NgArray<int> & finds, int nfinbox);
+			     FlatArray<int> finds); // , int nfinbox);