Delaunay2d also for small sub-domains

libsrc/geom2d/genmesh2d.cpp

@@ -687,7 +687,7 @@ namespace netgen
         
         t_points.Stop();
 
-        if(mp.delaunay2d && cnt>100)
+        if(mp.delaunay2d && cnt>1)
           meshing.Delaunay(*mesh, domnr, mp);
         else
         {
@@ -707,6 +707,8 @@ namespace netgen
           mesh->SetMaterial (domnr, material);
       }
 
+    mesh->Compress();
+
     mp.quad = hquad;
 
 

libsrc/meshing/adfront2.hpp

@@ -206,7 +206,7 @@ public:
 
   const FrontLine & GetLine (int nr) { return lines[nr]; }
   const FrontPoint2 & GetPoint (int nr) { return points[nr]; }
-
+  const auto & GetLines () const { return lines; }
 
   ///
   int SelectBaseLine (Point<3> & p1, Point<3> & p2, 

libsrc/meshing/delaunay2d.cpp

@@ -454,8 +454,14 @@ namespace netgen
     bool changed;
     do 
       {
-	mesh.LocalHFunction().ClearFlags();
-	
+        static Timer tcf("clear flags");
+        tcf.Start();
+        // mesh.LocalHFunction().ClearFlags();
+        mesh.LocalHFunction().ClearRootFlags();
+	tcf.Stop();
+        
+        static Timer tcut("tcut");
+        tcut.Start();
 	for (int i = 0; i < adfront.GetNFL(); i++)
 	  {
 	    const FrontLine & line = adfront.GetLine(i);
@@ -469,7 +475,7 @@ namespace netgen
 	    
 	    mesh.LocalHFunction().CutBoundary (bbox); 
 	  }
-	
+	tcut.Stop();
 
 	mesh.LocalHFunction().FindInnerBoxes (&adfront, NULL);
 	
@@ -506,11 +512,11 @@ namespace netgen
   
     {
       int i = 0;
-      while (npoints.Size() % prims[i] == 0) i++;
+      if (npoints.Size())
+        while (npoints.Size() % prims[i] == 0) i++;
       prim = prims[i];
     }
 
-    
     for (int i = 0; i < npoints.Size(); i++)
       {
         size_t hi = (size_t(prim) * size_t(i)) % npoints.Size();
@@ -570,7 +576,8 @@ namespace netgen
       // outer points : smooth mesh-grading
     npoints.SetSize(0);
     loch2.GetOuterPoints (npoints);
-    
+
+    /*
     for (int i = 1; i <= npoints.Size(); i++)
       {
 	if (meshbox.IsIn (npoints.Get(i)))
@@ -579,7 +586,14 @@ namespace netgen
 	    adfront.AddPoint (npoints.Get(i), gpnum);
 	  }
       }  
+    */
 
+    for (const Point<3> p : npoints)
+      if (meshbox.IsIn(p))
+        {
+          PointIndex gpnum = mesh.AddPoint (p);
+          adfront.AddPoint (p, gpnum);
+        }
     timer4.Stop();
   }
 
@@ -589,6 +603,9 @@ namespace netgen
   void Meshing2 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp)
   {
     static Timer timer("Meshing2::Delaunay");
+    static Timer t1("Meshing2::Delaunay1");
+    static Timer t2("Meshing2::Delaunay2");
+    static Timer t3("Meshing2::Delaunay3");
     static Timer timer_addpoints("add points");
     RegionTimer reg (timer);
 
@@ -600,6 +617,7 @@ namespace netgen
 
     auto last_point_blockfill = mesh.Points().Range().Next();
 
+    t1.Start();
     // Bounding box for starting trig in delaunay
     Box<2> bbox (Box<2>::EMPTY_BOX);
 
@@ -615,10 +633,14 @@ namespace netgen
 
     for (int i = 0; i < mesh.LockedPoints().Size(); i++)
       bbox.Add (P2(mesh.Point (mesh.LockedPoints()[i])));
+    t1.Stop();
 
+    t2.Start();
     Array<PointIndex> old_points;
     BitArray add_point(mesh.Points().Size()+1);
+    Array<PointIndex> addpoints;
     add_point.Clear();
+    /*
     for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
     {
       const auto & s = mesh[si];
@@ -628,7 +650,28 @@ namespace netgen
         add_point.SetBit(s[1]);
       }
     }
+    */
+    /*
+    for (int i = 0; i < adfront.GetNFL(); i++)
+      {
+	const FrontLine & line = adfront.GetLine(i);
+        for (int j = 0; j < 2; j++)
+          add_point.SetBit (adfront.GetGlobalIndex (line.L()[j]))adfront.GetGlobalIndex (line.L()[j]));
+      }
+    */
+    for (const auto & line : adfront.GetLines())
+      for (int j = 0; j < 2; j++)
+        {
+          PointIndex pnum = adfront.GetGlobalIndex (line.L()[j]);
+          if (!add_point.Test(pnum))
+            addpoints.Append(pnum);
+          add_point.SetBit (pnum);
+        }
+      
+    
+    t2.Stop();
 
+    t3.Start();
     Mesh tempmesh;
     tempmesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
     tempmesh.AddFaceDescriptor (FaceDescriptor (2, 1, 0, 0));
@@ -637,8 +680,11 @@ namespace netgen
     Array<PointIndex, PointIndex> compress;
     Array<PointIndex, PointIndex> icompress(mesh.Points().Size());
 
+    /*
     for(auto pi : mesh.Points().Range())
       if(add_point.Test(pi))
+    */
+    for (PointIndex pi : addpoints)
       {
         icompress[pi] = tempmesh.AddPoint(mesh[pi]);
         compress.Append(pi);
@@ -649,7 +695,7 @@ namespace netgen
         icompress[pi] = tempmesh.AddPoint(mesh[pi]);
         compress.Append(pi);
       }
-
+    t3.Stop();
     // DelaunayMesh adds surrounding trig (don't add the last 3 points to delaunay AGAIN!
     auto tempmesh_points = tempmesh.Points().Range();
 
@@ -676,7 +722,10 @@ namespace netgen
 
     timer_addpoints.Stop();
 
+    static Timer taddseg("addseg");
+    taddseg.Start();
 
+    /*
     for (auto seg : mesh.LineSegments())
     {
       if ( seg.domin == domainnr || seg.domout == domainnr )
@@ -690,7 +739,19 @@ namespace netgen
         tempmesh.AddSegment(seg);
       }
     }
-
+    */
+    for (const auto & line : adfront.GetLines())
+      {
+        Segment seg;
+        for (int j = 0; j < 2; j++)
+          seg[j] = icompress [adfront.GetGlobalIndex (line.L()[j])];
+        seg.domin = domainnr;
+        seg.domout = 0;
+        tempmesh.AddSegment(seg);
+      }
+           
+    taddseg.Stop();
+    
     for (auto & trig : dmesh.GetElements())
     {
       if (trig[0] < 0) continue;
@@ -883,7 +944,7 @@ namespace netgen
       }
     }
 
-   mesh.Compress();
+    // mesh.Compress();  // don't compress whole mesh after every sub-domain
   }
 
 }

libsrc/meshing/localh.cpp

@@ -381,7 +381,12 @@ namespace netgen
           return;
       }
 
-    box->flags.cutboundary = 1;
+    if (!box->flags.cutboundary)
+      for (int i = 0; i < 8; i++)
+        if (box->childs[i])
+          box->childs[i]->flags.cutboundary = false;
+    
+    box->flags.cutboundary = true;
     for (int i = 0; i < 8; i++)
       if (box->childs[i])
 	CutBoundaryRec (pmin, pmax, box->childs[i]);
@@ -547,13 +552,19 @@ namespace netgen
   void LocalH :: FindInnerBoxes (AdFront2 * adfront,
 				 int (*testinner)(const Point<2> & p1))
   {
-    static int timer = NgProfiler::CreateTimer ("LocalH::FindInnerBoxes 2d");
-    NgProfiler::RegionTimer reg (timer);
+    static Timer t("LocalH::FindInnerBoxes 2d"); RegionTimer reg (t);
+    static Timer trec("LocalH::FindInnerBoxes 2d - rec");
+    static Timer tinit("LocalH::FindInnerBoxes 2d - init");     
 
+    /*
+    tinit.Start();
     for (int i = 0; i < boxes.Size(); i++)
       boxes[i] -> flags.isinner = 0;
-
+    tinit.Stop();
+    */
+    
     root->flags.isinner = 0;
+    root->flags.cutboundary = true;
 
     Point<2> rpmid(root->xmid[0], root->xmid[1]);   // , root->xmid[2]);
     Vec<2> rv(root->h2, root->h2);
@@ -582,7 +593,8 @@ namespace netgen
 	faceboxes[i].Set (Point<2> (p1(0), p1(1)));
 	faceboxes[i].Add (Point<2> (p2(0), p2(1)));
       }
-  
+
+    RegionTimer regrc(trec);
     for (int i = 0; i < 8; i++)
       FindInnerBoxesRec2 (root->childs[i], adfront, faceboxes, faceinds); // , nf);
   }
@@ -602,11 +614,13 @@ namespace netgen
       {
 	box->flags.isinner = father->flags.isinner;
 	box->flags.pinner = father->flags.pinner;
+        box->flags.cutboundary = false;
       }
     else
       {        
 	if (father->flags.isinner)
           {
+            cout << "how is this possible ???" << endl;
             box->flags.pinner = 1;
           }
 	else
@@ -672,9 +686,11 @@ namespace netgen
               }
           }
       }
-    
-    for (int i = 0; i < 8; i++)
-      FindInnerBoxesRec2 (box->childs[i], adfront, faceboxes, faceinds.Range(0,iused));
+
+
+    if (box->flags.isinner || box->flags.cutboundary)
+      for (int i = 0; i < 8; i++)
+        FindInnerBoxesRec2 (box->childs[i], adfront, faceboxes, faceinds.Range(0,iused));
   }
 
 
@@ -720,6 +736,13 @@ namespace netgen
 	ClearFlagsRec (box->childs[i]);
   }
 
+  void LocalH :: ClearRootFlags ()
+  {
+    root->flags.cutboundary = false;
+    root->flags.isinner = false;
+  }
+
+  
   void LocalH :: ClearFlagsRec (GradingBox * box)
   {
     box->flags.cutboundary = 0;
@@ -746,13 +769,17 @@ namespace netgen
       }
   }
 
-  void LocalH :: GetInnerPoints (NgArray<Point<3> > & points)
+  void LocalH :: GetInnerPoints (NgArray<Point<3> > & points) const
   {
+    static Timer t("GetInnerPoints"); RegionTimer reg(t);
     if (dimension == 2)
       {
+        GetInnerPointsRec (root, points);
+        /*
         for (int i = 0; i < boxes.Size(); i++)
           if (boxes[i] -> flags.isinner && boxes[i] -> HasChilds())
             points.Append ( boxes[i] -> PMid() );
+        */
       }
     else
       {
@@ -763,7 +790,18 @@ namespace netgen
           
   }
 
+  void LocalH :: GetInnerPointsRec (const GradingBox * box, NgArray<Point<3> > & points) const
+  {
+    if (box -> flags.isinner && box -> HasChilds())
+      points.Append ( box -> PMid() );
 
+    if (box->flags.isinner || box->flags.cutboundary)
+      for (int i = 0; i < 8; i++)
+        if (box->childs[i])
+          GetInnerPointsRec (box->childs[i], points);
+  }
+
+  
   void LocalH :: GetOuterPoints (NgArray<Point<3> > & points)
   {
     for (int i = 0; i < boxes.Size(); i++)

libsrc/meshing/localh.hpp

@@ -121,11 +121,14 @@ namespace netgen
     void ClearFlags ()
     { ClearFlagsRec(root); }
 
+    void ClearRootFlags ();
+
     /// widen refinement zone
     void WidenRefinement ();
 
     /// get points in inner elements
-    void GetInnerPoints (NgArray<Point<3> > & points);
+    void GetInnerPoints (NgArray<Point<3> > & points) const;
+    void GetInnerPointsRec (const GradingBox * box, NgArray<Point<3> > & points) const;
 
     /// get points in outer closure
     void GetOuterPoints (NgArray<Point<3> > & points);