prepare LocalH tree before blockfill sequentially

libsrc/meshing/meshfunc.cpp

@@ -227,6 +227,26 @@ namespace netgen
    }
   }
 
+  void PrepareForBlockFillLocalH(MeshingData & md)
+  {
+    static Timer t("PrepareForBlockFillLocalH"); RegionTimer rt(t);
+    md.meshing = make_unique<Meshing3>(nullptr);
+
+    auto & mesh = md.mesh;
+
+    mesh.CalcSurfacesOfNode();
+    mesh.FindOpenElements(md.domain);
+
+    for (PointIndex pi : mesh.Points().Range())
+       md.meshing->AddPoint (mesh[pi], pi);
+
+    for (int i = 1; i <= mesh.GetNOpenElements(); i++)
+       md.meshing->AddBoundaryElement (mesh.OpenElement(i));
+
+    if (mesh.HasLocalHFunction())
+        md.meshing->PrepareBlockFillLocalH(mesh, md.mp);
+  }
+
 
   void MeshDomain( MeshingData & md)
   {
@@ -237,15 +257,7 @@ namespace netgen
     if(mp.only3D_domain_nr && mp.only3D_domain_nr != domain)
         return;
 
-    mesh.CalcSurfacesOfNode();
-    mesh.FindOpenElements(domain);
 
-    md.meshing = make_unique<Meshing3>(nullptr);
-    for (PointIndex pi : mesh.Points().Range())
-       md.meshing->AddPoint (mesh[pi], pi);
-
-    for (int i = 1; i <= mesh.GetNOpenElements(); i++)
-       md.meshing->AddBoundaryElement (mesh.OpenElement(i));
 
 
    if (mp.delaunay && mesh.GetNOpenElements())
@@ -261,9 +273,18 @@ namespace netgen
          mesh.GetNE(), " elements");
    }
 
-   // mesh.Save("before_findopenels.vol");
+   Box<3> domain_bbox( Box<3>::EMPTY_BOX ); 
-   // mesh.CalcSurfacesOfNode();
+   
-   // mesh.FindOpenElements(domain);
+   for (auto & sel : mesh.SurfaceElements())
+     {
+       if (sel.IsDeleted() ) continue;
+
+       for (auto pi : sel.PNums())
+         domain_bbox.Add (mesh[pi]);
+     }
+   domain_bbox.Increase (0.01 * domain_bbox.Diam());
+
+   mesh.FindOpenElements(domain);
 
    int cntsteps = 0;
    int meshed;
@@ -273,7 +294,7 @@ namespace netgen
          if (multithread.terminate)
            break;
          
-         mesh.FindOpenElements();
+         mesh.FindOpenElements(domain);
          PrintMessage (5, mesh.GetNOpenElements(), " open faces");
          GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol");
          cntsteps++;
@@ -289,15 +310,14 @@ namespace netgen
          Array<PointIndex, PointIndex> glob2loc(mesh.GetNP());
          glob2loc = PointIndex::INVALID;
 
+         for (PointIndex pi : mesh.Points().Range())
+           if (domain_bbox.IsIn (mesh[pi]))
+               glob2loc[pi] = meshing.AddPoint (mesh[pi], pi);
+
          for (auto sel : mesh.OpenElements() )
          {
            for(auto & pi : sel.PNums())
-             {
-               if(!glob2loc[pi].IsValid())
-                 glob2loc[pi] = meshing.AddPoint (mesh[pi], pi);
-
                pi = glob2loc[pi];
-             }
            meshing.AddBoundaryElement (sel);
          }
 
@@ -312,7 +332,7 @@ namespace netgen
          
 
          mesh.CalcSurfacesOfNode();
-         mesh.FindOpenElements();
+         mesh.FindOpenElements(domain);
 
          // teterrpow = 2;
          if (mesh.GetNOpenElements() != 0)
@@ -358,7 +378,7 @@ namespace netgen
      {
         PrintMessage (3, "Check subdomain ", domain, " / ", mesh.GetNDomains());
 
-        mesh.FindOpenElements();
+        mesh.FindOpenElements(domain);
 
         bool res = (mesh.CheckConsistentBoundary() != 0);
         if (res)
@@ -710,14 +730,12 @@ namespace netgen
          CloseOpenQuads( md[i] );
        });
 
+     for(auto & md_ : md)
+         PrepareForBlockFillLocalH(md_);
+
      ParallelFor( md.Range(), [&](int i)
        {
-         // try {
+         MeshDomain(md[i]);
-           MeshDomain(md[i]);
-           // }
-         // catch( const NgException & e ) {
-         //   md[i].mesh.Save("error_"+ToString(i)+".vol");
-         // }
        });
 
      MergeMeshes(mesh3d, md);

libsrc/meshing/meshing3.cpp

@@ -180,6 +180,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
   // static int meshing3_timer_d = NgProfiler::CreateTimer ("Meshing3::GenerateMesh d");
   // NgProfiler::RegionTimer reg (meshing3_timer);
 
+  cout << "start tet meshing with " << adfront->GetNP() << " points and " << adfront->GetNF() << " faces " << endl;
 
   NgArray<Point3d, PointIndex::BASE> locpoints;      // local points
   NgArray<MiniElement2d> locfaces;                   // local faces
@@ -1094,11 +1095,10 @@ static int TestSameSide (const Point3d & p1, const Point3d & p2)
 */
 
 
-
+void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh, 
-void Meshing3 :: BlockFillLocalH (Mesh & mesh, 
 				  const MeshingParameters & mp)
 {
-  static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t);
+  static Timer t("Mesing3::PrepareBlockFillLocalH"); RegionTimer reg(t);
   
   double filldist = mp.filldist;
   
@@ -1107,10 +1107,90 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
   PrintMessage (3, "blockfill local h");
 
 
-  NgArray<Point<3> > npoints;
-  
   adfront -> CreateTrees();
 
+  double maxh = 0;
+
+  for (int i = 1; i <= adfront->GetNF(); i++)
+    {
+      const MiniElement2d & el = adfront->GetFace(i);
+      for (int j = 1; j <= 3; j++)
+	{
+	  const Point3d & p1 = adfront->GetPoint (el.PNumMod(j));
+	  const Point3d & p2 = adfront->GetPoint (el.PNumMod(j+1));
+
+	  double hi = Dist (p1, p2);
+	  if (hi > maxh) maxh = hi;
+	}
+    }
+
+
+  if (mp.maxh < maxh) maxh = mp.maxh;
+
+  // auto loch_ptr = mesh.LocalHFunction().Copy();
+  // auto & loch = *loch_ptr;
+  auto & loch = mesh.LocalHFunction();
+
+  bool changed;
+  static Timer t1("loop1");
+  t1.Start();
+  do 
+    {
+      loch.ClearFlags();
+
+      static Timer tbox("adfront-bbox");
+      tbox.Start();
+      for (int i = 1; i <= adfront->GetNF(); i++)
+	{
+	  const MiniElement2d & el = adfront->GetFace(i);
+	  
+	  Box<3> bbox (adfront->GetPoint (el[0]));
+	  bbox.Add (adfront->GetPoint (el[1]));
+	  bbox.Add (adfront->GetPoint (el[2]));
+
+
+	  double filld = filldist * bbox.Diam();
+	  bbox.Increase (filld);
+      
+      	  loch.CutBoundary (bbox); // .PMin(), bbox.PMax());
+	}
+      tbox.Stop();
+
+      //      locadfront = adfront;
+      loch.FindInnerBoxes (adfront, NULL);
+
+      npoints.SetSize(0);
+      loch.GetInnerPoints (npoints);
+
+      changed = false;
+      for (int i = 1; i <= npoints.Size(); i++)
+	{
+	  if (loch.GetH(npoints.Get(i)) > 1.5 * maxh)
+	    {
+	      loch.SetH (npoints.Get(i), maxh);
+	      changed = true;
+	    }
+	}
+    }
+  while (changed);
+  t1.Stop();
+
+
+
+}
+
+void Meshing3 :: BlockFillLocalH (Mesh & mesh, 
+				  const MeshingParameters & mp)
+{
+  static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t);
+  // PrepareBlockFillLocalH(mesh, mp);
+  
+  double filldist = mp.filldist;
+  
+  // (*testout) << "blockfill local h" << endl;
+  // (*testout) << "rel filldist = " << filldist << endl;
+  PrintMessage (3, "blockfill local h");
+
   Box<3> bbox ( Box<3>::EMPTY_BOX );
   double maxh = 0;
 
@@ -1144,47 +1224,6 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
 
   if (mp.maxh < maxh) maxh = mp.maxh;
 
-  auto loch_ptr = mesh.LocalHFunction().Copy();
-  auto & loch = *loch_ptr;
-
-  bool changed;
-  do 
-    {
-      loch.ClearFlags();
-
-      for (int i = 1; i <= adfront->GetNF(); i++)
-	{
-	  const MiniElement2d & el = adfront->GetFace(i);
-	  
-	  Box<3> bbox (adfront->GetPoint (el[0]));
-	  bbox.Add (adfront->GetPoint (el[1]));
-	  bbox.Add (adfront->GetPoint (el[2]));
-
-
-	  double filld = filldist * bbox.Diam();
-	  bbox.Increase (filld);
-      
-      	  loch.CutBoundary (bbox); // .PMin(), bbox.PMax());
-	}
-
-      //      locadfront = adfront;
-      loch.FindInnerBoxes (adfront, NULL);
-
-      npoints.SetSize(0);
-      loch.GetInnerPoints (npoints);
-
-      changed = false;
-      for (int i = 1; i <= npoints.Size(); i++)
-	{
-	  if (loch.GetH(npoints.Get(i)) > 1.5 * maxh)
-	    {
-	      loch.SetH (npoints.Get(i), maxh);
-	      changed = true;
-	    }
-	}
-    }
-  while (changed);
-
   if (debugparam.slowchecks)
     (*testout) << "Blockfill with points: " << endl;
   for (int i = 1; i <= npoints.Size(); i++)
@@ -1211,6 +1250,8 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
 
   // find outer points
   
+  static Timer tloch2("build loch2");
+  tloch2.Start();
   loch2.ClearFlags();
 
   for (int i = 1; i <= adfront->GetNF(); i++)
@@ -1248,6 +1289,7 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
       // loch2.CutBoundary (pmin, pmax);
       loch2.CutBoundary (Box<3> (pmin, pmax)); // pmin, pmax);
     }
+  tloch2.Stop();
 
   // locadfront = adfront;
   loch2.FindInnerBoxes (adfront, NULL);

libsrc/meshing/meshing3.hpp

@@ -28,6 +28,7 @@ class Meshing3
   NgArray<char*> problems;
   /// tolerance criterion
   double tolfak;
+  NgArray<Point<3> > npoints;
 public:
   /// 
   Meshing3 (const string & rulefilename); 
@@ -63,6 +64,7 @@ public:
   ///
   void BlockFill (Mesh & mesh, double gh);
   ///
+  void PrepareBlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
   void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp);
 
   /// uses points of adfront, and puts new elements into mesh