Merge branch 'use_new_btree' into 'master'

Use new BoxTree, activate search tree in STLBoundary

See merge request jschoeberl/netgen!259

libsrc/gprim/adtree.hpp

@@ -741,89 +741,355 @@ public:
   const ADTree3 & Tree() const { return *tree; };
 };
 
+template<int dim, typename T=INDEX>
+class BoxTree
+{
+public:
+  // Number of entries per leaf
+  static constexpr int N = 100;
 
-  template <int dim, typename T = INDEX>
+  struct Node;
-  class BoxTree
+
+  struct Leaf
   {
-    T_ADTree<2*dim,T> * tree;
+    Point<2*dim> p[N];
-    Point<dim> boxpmin, boxpmax;
+    T index[N];
-  public:
+    int n_elements;
-    BoxTree (const Box<dim> & abox)
+
+    Leaf() : n_elements(0)
+    { }
+
+    void Add( ClosedHashTable<T, Leaf*> &leaf_index, const Point<2*dim> &ap, T aindex )
+      {
+        p[n_elements] = ap;
+        index[n_elements] = aindex;
+        n_elements++;
+        leaf_index[aindex] = this;
+      }
+  };
+
+  struct Node
+  {
+    union
     {
-      boxpmin = abox.PMin();
+      Node *children[2];
-      boxpmax = abox.PMax();
+      Leaf *leaf;
+    };
+    double sep;
+    int level;
+
+    Node()
+        : children{nullptr,nullptr}
+    { }
+
+    ~Node()
+     { }
+
+    Leaf *GetLeaf() const
+      {
+        return children[1] ? nullptr : leaf;
+      }
+  };
+
+private:
+  Node root;
+
+  ClosedHashTable<T, Leaf*> leaf_index;
+
+  Point<dim> global_min, global_max;
+  double tol;
+  size_t n_leaves;
+  size_t n_nodes;
+  BlockAllocator ball_nodes;
+  BlockAllocator ball_leaves;
+
+public:
+
+  BoxTree (const Point<dim> & pmin, const Point<dim> & pmax)
+      : global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf))
+    {
+      root.leaf = (Leaf*) ball_leaves.Alloc(); new (root.leaf) Leaf();
+      root.level = 0;
+      tol = 1e-7 * Dist(pmax, pmin);
+    }
+
+  BoxTree (const Box<dim> & box)
+      : BoxTree(box.PMin(), box.PMax())
+    { }
+
+  size_t GetNLeaves()
+    {
+      return n_leaves;
+    }
+
+  size_t GetNNodes()
+    {
+      return n_nodes;
+    }
+
+  template<typename TFunc>
+  void GetFirstIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
+          TFunc func=[](auto pi){return false;}) const
+    {
+      // static Timer timer("BoxTree::GetIntersecting"); RegionTimer rt(timer);
+      // static Timer timer1("BoxTree::GetIntersecting-LinearSearch");
+      ArrayMem<const Node*, 100> stack;
+      ArrayMem<int, 100> dir_stack;
+
+
       Point<2*dim> tpmin, tpmax;
-      for (int i = 0; i < dim; i++)
+
+      for (size_t i : IntRange(dim))
         {
-          tpmin(i) = tpmin(i+dim) = boxpmin(i);
+          tpmin(i) = global_min(i);
-          tpmax(i) = tpmax(i+dim) = boxpmax(i);
+          tpmax(i) = pmax(i)+tol;
-        }
-      tree = new T_ADTree<2*dim,T> (tpmin, tpmax);
-    }
 
-    BoxTree (const Point<dim> & apmin, const Point<dim> & apmax)
+          tpmin(i+dim) = pmin(i)-tol;
-    {
+          tpmax(i+dim) = global_max(i);
-      boxpmin = apmin;
-      boxpmax = apmax;
-      Point<2*dim> tpmin, tpmax;
-      for (int i = 0; i < dim; i++)
-        {
-          tpmin(i) = tpmin(i+dim) = boxpmin(i);
-          tpmax(i) = tpmax(i+dim) = boxpmax(i);
-        }
-      tree = new T_ADTree<2*dim,T> (tpmin, tpmax);
-    }
-      
-    ~BoxTree ()
-    {
-      delete tree;
-    }
-    
-    void Insert (const Point<dim> & bmin, const Point<dim> & bmax, T pi)
-    {
-      Point<2*dim> tp;
-      
-      for (size_t i = 0; i < dim; i++)
-        {
-          tp(i) = bmin(i);
-          tp(i+dim) = bmax(i);
         }
 
-      tree->Insert (tp, pi);
+      stack.SetSize(0);
+      stack.Append(&root);
+      dir_stack.SetSize(0);
+      dir_stack.Append(0);
+
+      while(stack.Size())
+        {
+          const Node *node = stack.Last();
+          stack.DeleteLast();
+
+          int dir = dir_stack.Last();
+          dir_stack.DeleteLast();
+
+          if(Leaf *leaf = node->GetLeaf())
+            {
+              //               RegionTimer rt1(timer1);
+              for (auto i : IntRange(leaf->n_elements))
+                {
+                  bool intersect = true;
+                  const auto p = leaf->p[i];
+
+                  for (int d = 0; d < dim; d++)
+                      if (p[d] > tpmax[d])
+                          intersect = false;
+                  for (int d = dim; d < 2*dim; d++)
+                      if (p[d] < tpmin[d])
+                          intersect = false;
+                  if(intersect)
+                      if(func(leaf->index[i])) return;
+                }
+            }
+          else
+            {
+              int newdir = dir+1;
+              if(newdir==2*dim) newdir = 0;
+              if (tpmin[dir] <= node->sep)
+                {
+                  stack.Append(node->children[0]);
+                  dir_stack.Append(newdir);
+                }
+              if (tpmax[dir] >= node->sep)
+                {
+                  stack.Append(node->children[1]);
+                  dir_stack.Append(newdir);
+                }
+            }
+        }
     }
 
-    void Insert (const Box<dim> & box, T pi)
+  void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
+          NgArray<T> & pis) const
+    {
+      pis.SetSize(0);
+      GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
+    }
+
+  void Insert (const Box<dim> & box, T pi)
     {
       Insert (box.PMin(), box.PMax(), pi);
     }
 
-    void DeleteElement (T pi) 
+  void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
     {
-      tree->DeleteElement(pi);
+      // static Timer timer("BoxTree::Insert"); RegionTimer rt(timer);
-    }
+      int dir = 0;
-
+      Point<2*dim> p;
-    void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax, 
+      for (auto i : IntRange(dim))
-                          NgArray<T> & pis) const
-    {
-      Point<2*dim> tpmin, tpmax;
-      double tol = Tolerance();
-      for (size_t i = 0; i < dim; i++)
         {
-          tpmin(i) = boxpmin(i);
+          p(i) = pmin[i];
-          tpmax(i) = pmax(i)+tol;
+          p(i+dim) = pmax[i];
-        
-          tpmin(i+dim) = pmin(i)-tol;
-          tpmax(i+dim) = boxpmax(i);
         }
 
-      tree->GetIntersecting (tpmin, tpmax, pis);
+      Node * node = &root;
+      Leaf * leaf = node->GetLeaf();
+
+      // search correct leaf to add point
+      while(!leaf)
+        {
+          node = p[dir] < node->sep ? node->children[0] : node->children[1];
+          dir++;
+          if(dir==2*dim) dir = 0;
+          leaf = node->GetLeaf();
+        }
+
+      // add point to leaf
+      if(leaf->n_elements < N)
+          leaf->Add(leaf_index, p,pi);
+      else // assume leaf->n_elements == N
+        {
+          // add two new nodes and one new leaf
+          int n_elements = leaf->n_elements;
+          ArrayMem<double, N> coords(n_elements);
+          ArrayMem<int, N> order(n_elements);
+
+          // separate points in two halves, first sort all coordinates in direction dir
+          for (auto i : IntRange(n_elements))
+            {
+              order[i] = i;
+              coords[i] = leaf->p[i][dir];
+            }
+
+          QuickSortI(coords, order);
+          int isplit = N/2;
+          Leaf *leaf1 = (Leaf*) ball_leaves.Alloc(); new (leaf1) Leaf();
+          Leaf *leaf2 = (Leaf*) ball_leaves.Alloc(); new (leaf2) Leaf();
+
+          for (auto i : order.Range(isplit))
+              leaf1->Add(leaf_index, leaf->p[i], leaf->index[i] );
+          for (auto i : order.Range(isplit, N))
+              leaf2->Add(leaf_index, leaf->p[i], leaf->index[i] );
+
+          Node *node1 = (Node*) ball_nodes.Alloc(); new (node1) Node();
+          node1->leaf = leaf1;
+          node1->level = node->level+1;
+
+          Node *node2 = (Node*) ball_nodes.Alloc(); new (node2) Node();
+          node2->leaf = leaf2;
+          node2->level = node->level+1;
+
+          node->children[0] = node1;
+          node->children[1] = node2;
+          node->sep = 0.5 * (leaf->p[order[isplit-1]][dir] + leaf->p[order[isplit]][dir]);
+
+          // add new point to one of the new leaves
+          if (p[dir] < node->sep)
+              leaf1->Add( leaf_index, p, pi );
+          else
+              leaf2->Add( leaf_index, p, pi );
+
+          ball_leaves.Free(leaf);
+          n_leaves++;
+          n_nodes+=2;
+        }
     }
 
+  void DeleteElement (T pi)
+    {
+      // static Timer timer("BoxTree::DeleteElement"); RegionTimer rt(timer);
+      Leaf *leaf = leaf_index[pi];
+      leaf_index.Delete(pi);
+      auto & n_elements = leaf->n_elements;
+      auto & index = leaf->index;
+      auto & p = leaf->p;
 
-    double Tolerance() const { return 1e-7 * Dist(boxpmax, boxpmin); } // single precision
+      for (auto i : IntRange(n_elements))
-    const auto & Tree() const { return *tree; };
+        {
-    auto & Tree() { return *tree; };
+          if(index[i] == pi)
-  };
+            {
+              n_elements--;
+              if(i!=n_elements)
+                {
+                  index[i] = index[n_elements];
+                  p[i] = p[n_elements];
+                }
+              return;
+            }
+        }
+    }
+};
+
+//   template <int dim, typename T = INDEX>
+//   class BoxTree
+//   {
+//     T_ADTree<2*dim,T> * tree;
+//     Point<dim> boxpmin, boxpmax;
+//   public:
+//     BoxTree (const Box<dim> & abox)
+//     {
+//       boxpmin = abox.PMin();
+//       boxpmax = abox.PMax();
+//       Point<2*dim> tpmin, tpmax;
+//       for (int i = 0; i < dim; i++)
+//         {
+//           tpmin(i) = tpmin(i+dim) = boxpmin(i);
+//           tpmax(i) = tpmax(i+dim) = boxpmax(i);
+//         }
+//       tree = new T_ADTree<2*dim,T> (tpmin, tpmax);
+//     }
+//
+//     BoxTree (const Point<dim> & apmin, const Point<dim> & apmax)
+//     {
+//       boxpmin = apmin;
+//       boxpmax = apmax;
+//       Point<2*dim> tpmin, tpmax;
+//       for (int i = 0; i < dim; i++)
+//         {
+//           tpmin(i) = tpmin(i+dim) = boxpmin(i);
+//           tpmax(i) = tpmax(i+dim) = boxpmax(i);
+//         }
+//       tree = new T_ADTree<2*dim,T> (tpmin, tpmax);
+//     }
+//
+//     ~BoxTree ()
+//     {
+//       delete tree;
+//     }
+//
+//     void Insert (const Point<dim> & bmin, const Point<dim> & bmax, T pi)
+//     {
+//       Point<2*dim> tp;
+//
+//       for (size_t i = 0; i < dim; i++)
+//         {
+//           tp(i) = bmin(i);
+//           tp(i+dim) = bmax(i);
+//         }
+//
+//       tree->Insert (tp, pi);
+//     }
+//
+//     void Insert (const Box<dim> & box, T pi)
+//     {
+//       Insert (box.PMin(), box.PMax(), pi);
+//     }
+//
+//     void DeleteElement (T pi)
+//     {
+//       tree->DeleteElement(pi);
+//     }
+//
+//     void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
+//                           NgArray<T> & pis) const
+//     {
+//       Point<2*dim> tpmin, tpmax;
+//       double tol = Tolerance();
+//       for (size_t i = 0; i < dim; i++)
+//         {
+//           tpmin(i) = boxpmin(i);
+//           tpmax(i) = pmax(i)+tol;
+//
+//           tpmin(i+dim) = pmin(i)-tol;
+//           tpmax(i+dim) = boxpmax(i);
+//         }
+//
+//       tree->GetIntersecting (tpmin, tpmax, pis);
+//     }
+//
+//
+//     double Tolerance() const { return 1e-7 * Dist(boxpmax, boxpmin); } // single precision
+//     const auto & Tree() const { return *tree; };
+//     auto & Tree() { return *tree; };
+//   };
 
 }
 

libsrc/meshing/delaunay.cpp

@@ -5,265 +5,6 @@
 namespace netgen
 {
 
-  template<int dim, typename T=INDEX>
-  class BTree
-  {
-  public:
-    // Number of entries per leaf
-    static constexpr int N = 100;
-
-    struct Node;
-
-    struct Leaf
-    {
-      Point<2*dim> p[N];
-      T index[N];
-      int n_elements;
-
-      Leaf() : n_elements(0) {}
-
-      void Add( Array<Leaf*, INDEX> &leaf_index, const Point<2*dim> &ap, T aindex )
-      {
-        p[n_elements] = ap;
-        index[n_elements] = aindex;
-        n_elements++;
-        if(leaf_index.Size()<=aindex)
-            leaf_index.SetSize(2*aindex);
-        leaf_index[aindex] = this;
-      }
-    };
-
-    struct Node
-    {
-      union
-        {
-          Node *children[2];
-          Leaf *leaf;
-        };
-      double sep;
-      int level;
-
-      Node()
-          : children{nullptr,nullptr}
-      { }
-
-      ~Node()
-      { }
-
-      Leaf *GetLeaf() const
-        {
-          return children[1] ? nullptr : leaf;
-        }
-    };
-
-  private:
-    Node root;
-    Array<Leaf*, INDEX> leaf_index;
-    Point<dim> global_min, global_max;
-    double tol;
-    size_t n_leaves;
-    size_t n_nodes;
-    BlockAllocator ball_nodes;
-    BlockAllocator ball_leaves;
-
-  public:
-
-    BTree (const Point<dim> & pmin, const Point<dim> & pmax)
-        : global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf))
-    {
-      root.leaf = (Leaf*) ball_leaves.Alloc(); new (root.leaf) Leaf();
-      root.level = 0;
-      tol = 1e-7 * Dist(pmax, pmin);
-    }
-
-    size_t GetNLeaves()
-    {
-      return n_leaves;
-    }
-
-    size_t GetNNodes()
-    {
-      return n_nodes;
-    }
-
-    template<typename TFunc>
-    void GetFirstIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
-                               TFunc func=[](auto pi){return false;}) const
-    {
-      // static Timer timer("BTree::GetIntersecting"); RegionTimer rt(timer);
-      // static Timer timer1("BTree::GetIntersecting-LinearSearch");
-      static Array<const Node*> stack(1000);
-      static Array<int> dir_stack(1000);
-
-
-      Point<2*dim> tpmin, tpmax;
-
-      for (size_t i : IntRange(dim))
-        {
-          tpmin(i) = global_min(i);
-          tpmax(i) = pmax(i)+tol;
-        
-          tpmin(i+dim) = pmin(i)-tol;
-          tpmax(i+dim) = global_max(i);
-        }
-
-      stack.SetSize(0);
-      stack.Append(&root);
-      dir_stack.SetSize(0);
-      dir_stack.Append(0);
-
-      while(stack.Size())
-        {
-          const Node *node = stack.Last();
-          stack.DeleteLast();
-
-          int dir = dir_stack.Last();
-          dir_stack.DeleteLast();
-
-          if(Leaf *leaf = node->GetLeaf())
-            {
-//               RegionTimer rt1(timer1);
-              for (auto i : IntRange(leaf->n_elements))
-                {
-                  bool intersect = true;
-                  const auto p = leaf->p[i];
-
-                  for (int d = 0; d < dim; d++)
-                      if (p[d] > tpmax[d])
-                          intersect = false;
-                  for (int d = dim; d < 2*dim; d++)
-                      if (p[d] < tpmin[d])
-                          intersect = false;
-                  if(intersect)
-                    if(func(leaf->index[i])) return;
-                }
-            }
-          else
-            {
-              int newdir = dir+1;
-              if(newdir==2*dim) newdir = 0;
-              if (tpmin[dir] <= node->sep)
-                {
-                  stack.Append(node->children[0]);
-                  dir_stack.Append(newdir);
-                }
-              if (tpmax[dir] >= node->sep)
-                {
-                  stack.Append(node->children[1]);
-                  dir_stack.Append(newdir);
-                }
-            }
-        }
-    }
-
-    void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
-                          NgArray<T> & pis) const
-    {
-      pis.SetSize(0);
-      GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
-    }
-
-    void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
-    {
-      // static Timer timer("BTree::Insert"); RegionTimer rt(timer);
-      int dir = 0;
-      Point<2*dim> p;
-      for (auto i : IntRange(dim))
-        {
-          p(i) = pmin[i];
-          p(i+dim) = pmax[i];
-        }
-
-      Node * node = &root;
-      Leaf * leaf = node->GetLeaf();
-
-      // search correct leaf to add point
-      while(!leaf)
-        {
-          node = p[dir] < node->sep ? node->children[0] : node->children[1];
-          dir++;
-          if(dir==2*dim) dir = 0;
-          leaf = node->GetLeaf();
-        }
-
-      // add point to leaf
-      if(leaf->n_elements < N)
-          leaf->Add(leaf_index, p,pi);
-      else // assume leaf->n_elements == N
-        {
-          // add two new nodes and one new leaf
-          int n_elements = leaf->n_elements;
-          ArrayMem<double, N> coords(n_elements);
-          ArrayMem<int, N> order(n_elements);
-
-          // separate points in two halves, first sort all coordinates in direction dir
-          for (auto i : IntRange(n_elements))
-            {
-              order[i] = i;
-              coords[i] = leaf->p[i][dir];
-            }
-
-          QuickSortI(coords, order);
-          int isplit = N/2;
-          Leaf *leaf1 = (Leaf*) ball_leaves.Alloc(); new (leaf1) Leaf();
-          Leaf *leaf2 = (Leaf*) ball_leaves.Alloc(); new (leaf2) Leaf();
-
-          for (auto i : order.Range(isplit))
-              leaf1->Add(leaf_index, leaf->p[i], leaf->index[i] );
-          for (auto i : order.Range(isplit, N))
-              leaf2->Add(leaf_index, leaf->p[i], leaf->index[i] );
-
-          Node *node1 = (Node*) ball_nodes.Alloc(); new (node1) Node();
-          node1->leaf = leaf1;
-          node1->level = node->level+1;
-
-          Node *node2 = (Node*) ball_nodes.Alloc(); new (node2) Node();
-          node2->leaf = leaf2;
-          node2->level = node->level+1;
-
-          node->children[0] = node1;
-          node->children[1] = node2;
-          node->sep = 0.5 * (leaf->p[order[isplit-1]][dir] + leaf->p[order[isplit]][dir]);
-
-          // add new point to one of the new leaves
-          if (p[dir] < node->sep)
-              leaf1->Add( leaf_index, p, pi );
-          else
-              leaf2->Add( leaf_index, p, pi );
-
-          ball_leaves.Free(leaf);
-          n_leaves++;
-          n_nodes+=2;
-        }
-    }
-
-    void DeleteElement (T pi)
-    {
-      // static Timer timer("BTree::DeleteElement"); RegionTimer rt(timer);
-      Leaf *leaf = leaf_index[pi];
-      auto & n_elements = leaf->n_elements;
-      auto & index = leaf->index;
-      auto & p = leaf->p;
-
-      for (auto i : IntRange(n_elements))
-        {
-          if(index[i] == pi)
-            {
-              n_elements--;
-              if(i!=n_elements)
-                {
-                  index[i] = index[n_elements];
-                  p[i] = p[n_elements];
-                }
-              return;
-            }
-        }
-    }
-  };
-
-  typedef BTree<3> TBoxTree;
-//   typedef BoxTree<3> TBoxTree;
-
   static const int deltetfaces[][3] = 
     { { 1, 2, 3 },
       { 2, 0, 3 },
@@ -486,7 +227,7 @@ namespace netgen
   void AddDelaunayPoint (PointIndex newpi, const Point3d & newp, 
 			 NgArray<DelaunayTet> & tempels, 
 			 Mesh & mesh,
-			 TBoxTree & tettree,
+			 BoxTree<3> & tettree,
 			 MeshNB & meshnb,
 			 NgArray<Point<3> > & centers, NgArray<double> & radi2,
 			 NgArray<int> & connected, NgArray<int> & treesearch, 
@@ -894,7 +635,7 @@ namespace netgen
     pmin2 = pmin2 + 0.1 * (pmin2 - pmax2);
     pmax2 = pmax2 + 0.1 * (pmax2 - pmin2);
 
-    TBoxTree tettree(pmin2, pmax2);
+    BoxTree<3> tettree(pmin2, pmax2);
 
 
     tempels.Append (startel);

libsrc/stlgeom/stlgeomchart.cpp

@@ -146,8 +146,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
       chartbound.Clear();
       chartbound.SetChart(&chart);
       
-      chartbound.BuildSearchTree();  // different !!!
-
       if (!found) throw Exception("Make Atlas, no starttrig found"); 
 
       //find surrounding trigs
@@ -324,7 +322,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
       innerchartpts.SetSize(innerchartpoints.Size());
       for (size_t i = 0; i < innerchartpoints.Size(); i++)
         innerchartpts[i] = GetPoint(innerchartpoints[i]);
-      // chartbound.BuildSearchTree();  // different !!!
       
       // NgProfiler::StopTimer (timer2);
       // NgProfiler::StartTimer (timer3);

libsrc/stlgeom/stltool.cpp

@@ -1095,6 +1095,9 @@ void STLBoundary ::AddTriangle(const STLTriangle & t)
   segs[1] = INDEX_2(t[1], t[2]);
   segs[2] = INDEX_2(t[2], t[0]);
 
+  if(!searchtree)
+      BuildSearchTree();
+
   for (auto seg : segs)
     {
       STLBoundarySeg bseg(seg[0], seg[1], geometry->GetPoints(), chart);
@@ -1312,21 +1315,18 @@ bool STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3
 
 void STLBoundary :: BuildSearchTree()
 {
-  delete searchtree;
-  
   Box<2> box2d(Box<2>::EMPTY_BOX);
   Box<3> box3d = geometry->GetBoundingBox();
+
   for (size_t i = 0; i < 8; i++)
     box2d.Add ( chart->Project2d (box3d.GetPointNr(i)));
 
-  // comment to enable searchtree:
+  searchtree = make_unique<BoxTree<2,INDEX_2>> (box2d);
-  // searchtree = new BoxTree<2,INDEX_2> (box2d);
+//   searchtree = nullptr;
-  searchtree = nullptr;
 }
 
 void STLBoundary :: DeleteSearchTree()
 {
-  delete searchtree;
   searchtree = nullptr;
 }
 
@@ -1349,68 +1349,49 @@ bool STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
   Line2d l1 (p2d1, p2d2);
 
   double eps = 1e-6;
-  bool ok = true;
+
+  auto hasIntersection = [&] (auto i2) NETGEN_LAMBDA_INLINE
+    {
+      const STLBoundarySeg & seg = boundary_ht[i2];
+
+      if (seg.IsSmoothEdge()) return false;
+      if (!box2d.Intersect (seg.BoundingBox())) return false;
+
+      const Point<2> & sp1 = seg.P2D1();
+      const Point<2> & sp2 = seg.P2D2();
+
+      Line2d l2 (sp1, sp2);
+      double lam1, lam2;
+
+      int err = CrossPointBarycentric (l1, l2, lam1, lam2);
+      bool in1 = (lam1 > eps) && (lam1 < 1-eps);
+      bool on1 = (lam1 > -eps) && (lam1 < 1 + eps);
+      bool in2 = (lam2 > eps) && (lam2 < 1-eps);
+      bool on2 = (lam2 > -eps) && (lam2 < 1 + eps);
+
+      if(!err && ((on1 && in2) || (on2 && in1)))
+          return true;
+      return false;
+    };
 
   if (searchtree)
     {
-      NgArrayMem<INDEX_2,100> pis;
+      bool has_intersection = false;
-      searchtree -> GetIntersecting (box2d.PMin(), box2d.PMax(), pis);
+      searchtree -> GetFirstIntersecting (box2d.PMin(), box2d.PMax(),
-      for (auto i2 : pis)
+              [&] (auto i2) NETGEN_LAMBDA_INLINE
-        {
+              {
-          const STLBoundarySeg & seg = boundary_ht[i2];
+                  has_intersection = hasIntersection(i2);
-          
+                  return has_intersection;
-          if (seg.IsSmoothEdge()) continue;
+              });
-          if (!box2d.Intersect (seg.BoundingBox())) continue;
+      return !has_intersection;
-          
-          const Point<2> & sp1 = seg.P2D1();
-          const Point<2> & sp2 = seg.P2D2();
-          
-          Line2d l2 (sp1, sp2);
-          double lam1, lam2;
-          
-          int err = CrossPointBarycentric (l1, l2, lam1, lam2);
-          bool in1 = (lam1 > eps) && (lam1 < 1-eps);
-          bool on1 = (lam1 > -eps) && (lam1 < 1 + eps);
-          bool in2 = (lam2 > eps) && (lam2 < 1-eps);
-          bool on2 = (lam2 > -eps) && (lam2 < 1 + eps);
-
-          if(!err && ((on1 && in2) || (on2 && in1)))
-            {
-              ok = false;
-              break;
-            }
-        }
     }
-  
   else
     {
       for(auto [i2, seg] : boundary_ht)
-      {
+        if(hasIntersection(i2))
-        if (seg.IsSmoothEdge()) continue;
+              return false;
-        if (!box2d.Intersect (seg.BoundingBox())) continue;
+      return true;
-        
-        const Point<2> & sp1 = seg.P2D1();
-        const Point<2> & sp2 = seg.P2D2();
-        
-        Line2d l2 (sp1, sp2);
-        double lam1, lam2;
-
-        int err = CrossPointBarycentric (l1, l2, lam1, lam2);
-        
-        bool in1 = (lam1 > eps) && (lam1 < 1-eps);
-        bool on1 = (lam1 > -eps) && (lam1 < 1 + eps);
-        bool in2 = (lam2 > eps) && (lam2 < 1-eps);
-        bool on2 = (lam2 > -eps) && (lam2 < 1 + eps);
-        if(!err && ((on1 && in2) || (on2 && in1)))
-          {
-            ok = false;
-            break;
-          }
-      }
-    
     }
-
-  return ok;
 }
 
 

libsrc/stlgeom/stltool.hpp

@@ -191,10 +191,10 @@ private:
   const STLChart * chart;
   // NgArray<STLBoundarySeg> boundary;
   ClosedHashTable<INDEX_2, STLBoundarySeg> boundary_ht;  
-  BoxTree<2,INDEX_2> * searchtree = nullptr;
+  unique_ptr<BoxTree<2,INDEX_2>> searchtree;
 public:
   STLBoundary(STLGeometry * ageometry);
-  ~STLBoundary() { delete searchtree; }
+  ~STLBoundary() {}
 
   void Clear() { /* boundary.SetSize(0); */ boundary_ht = ClosedHashTable<INDEX_2,STLBoundarySeg>(); }
   void SetChart (const STLChart * achart) { chart = achart; }