copy localh tree (but skip parts outside of bounding box)

libsrc/meshing/localh.cpp

@@ -130,6 +130,54 @@ namespace netgen
     return lh;
   }
 
+  unique_ptr<LocalH> LocalH :: CopyRec ( const Box<3> & bbox, GradingBox *current )
+  {
+    static Timer t("LocalH::Copy with bounding box"); RegionTimer rt(t);
+    auto lh = make_unique<LocalH>(boundingbox, grading, dimension);
+    std::map<GradingBox*, GradingBox*> mapping;
+    lh->boxes.SetAllocSize(boxes.Size());
+
+    for(auto i : boxes.Range())
+    {
+      auto & b = *boxes[i];
+      auto h = b.H2();
+      Vec<3> vh = {h,h,h};
+      Box<3> box( b.PMid() - vh, b.PMid() + vh);
+      if(!box.Intersect(bbox))
+          continue;
+      lh->boxes.Append(new GradingBox());
+      auto & bnew = *lh->boxes.Last();
+      bnew.xmid[0] = b.xmid[0];
+      bnew.xmid[1] = b.xmid[1];
+      bnew.xmid[2] = b.xmid[2];
+      bnew.h2 = b.h2;
+      bnew.hopt = b.hopt;
+      bnew.flags = b.flags;
+      mapping[&b] = &bnew;
+    }
+
+    for(auto i : boxes.Range())
+    {
+      auto & b = *boxes[i];
+      if(mapping.count(&b)==0)
+          continue;
+
+      auto & bnew = *mapping[&b];
+      for(auto k : Range(8))
+      {
+        if(b.childs[k] && mapping.count(b.childs[k]))
+          bnew.childs[k] = mapping[b.childs[k]];
+      }
+
+      if(b.father && mapping.count(b.father))
+        bnew.father = mapping[b.father];
+    }
+
+    lh->root = mapping[root];
+    return lh;
+
+  }
+
   void LocalH :: Delete ()
   {
     root->DeleteChilds();

libsrc/meshing/localh.hpp

@@ -99,6 +99,7 @@ namespace netgen
     ~LocalH();
     ///
     unique_ptr<LocalH> Copy();
+    unique_ptr<LocalH> Copy( const Box<3> & bbox ) { return CopyRec(bbox, root); }
     ///
     void Delete();
     ///
@@ -191,6 +192,8 @@ namespace netgen
     ///
     void ConvexifyRec (GradingBox * box);
 
+    unique_ptr<LocalH> CopyRec( const Box<3> & bbox, GradingBox * current );
+
     friend ostream & operator<< (ostream & ost, const LocalH & loch);
   };
 

libsrc/meshing/meshfunc.cpp

@@ -68,7 +68,7 @@ namespace netgen
           ret[i].mesh = make_unique<Mesh>();
           auto & m = *ret[i].mesh;
 
-          m.SetLocalH(mesh.GetLocalH());
+          // m.SetLocalH(mesh.GetLocalH()->Copy());
 
           ipmap[i].SetSize(num_points);
           ipmap[i] = PointIndex::INVALID;

libsrc/meshing/meshing3.cpp

@@ -1109,6 +1109,7 @@ void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
   adfront -> CreateTrees();
 
   double maxh = 0;
+  Box<3> bounding_box(Box<3>::EMPTY_BOX);
 
   for (int i = 1; i <= adfront->GetNF(); i++)
     {
@@ -1117,6 +1118,7 @@ void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
 	{
 	  const Point3d & p1 = adfront->GetPoint (el.PNumMod(j));
 	  const Point3d & p2 = adfront->GetPoint (el.PNumMod(j+1));
+          bounding_box.Add(p1);
 
 	  double hi = Dist (p1, p2);
 	  if (hi > maxh) maxh = hi;
@@ -1126,9 +1128,9 @@ void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh,
 
   if (mp.maxh < maxh) maxh = mp.maxh;
 
-  // auto loch_ptr = mesh.LocalHFunction().Copy();
+  auto loch_ptr = mesh.LocalHFunction().Copy(bounding_box);
-  // auto & loch = *loch_ptr;
+  auto & loch = *loch_ptr;
-  auto & loch = mesh.LocalHFunction();
+  // auto & loch = mesh.LocalHFunction();
 
   bool changed;
   static Timer t1("loop1");