Special search tree for Delaunay (commented out)

- float instead of double
- Array for leaves instead of HashTable (the values of the tree are
        contiguous integer numbers)

libsrc/meshing/delaunay.cpp

@@ -1,9 +1,289 @@
 #include <mystdlib.h>
 #include "meshing.hpp"
 
-
 namespace netgen
 {
+  template<int dim, typename T=INDEX>
+  class DelaunayTree
+  {
+  public:
+    // Number of entries per leaf
+    static constexpr int N = 100;
+
+    struct Node;
+
+    struct Leaf
+    {
+      Point<2*dim, float> p[N];
+      T index[N];
+      int n_elements;
+      int nr;
+
+      Leaf() : n_elements(0)
+      { }
+
+
+      void Add( Array<Leaf*> &leaves, Array<T> &leaf_index, const Point<2*dim> &ap, T aindex )
+        {
+          p[n_elements] = ap;
+          index[n_elements] = aindex;
+          n_elements++;
+          if(leaf_index.Size()<aindex+1)
+              leaf_index.SetSize(aindex+1);
+          leaf_index[aindex] = nr;
+        }
+    };
+
+    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*> leaves;
+    Array<T> 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:
+
+    DelaunayTree (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.leaf->nr = 0;
+        leaves.Append(root.leaf);
+        root.level = 0;
+        tol = 1e-7 * Dist(pmax, pmin);
+      }
+
+    DelaunayTree (const Box<dim> & box)
+        : DelaunayTree(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("DelaunayTree::GetIntersecting"); RegionTimer rt(timer);
+        // static Timer timer1("DelaunayTree::GetIntersecting-LinearSearch");
+        ArrayMem<const Node*, 100> stack;
+        ArrayMem<int, 100> dir_stack;
+
+
+        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 Box<dim> & box, T pi)
+      {
+        Insert (box.PMin(), box.PMax(), pi);
+      }
+
+    void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
+      {
+        // static Timer timer("DelaunayTree::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(leaves, 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<float, 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();
+
+            leaf1->nr = leaf->nr;
+            leaf2->nr = leaves.Size();
+            leaves.Append(leaf2);
+            leaves[leaf1->nr] = leaf1;
+
+            for (auto i : order.Range(isplit))
+                leaf1->Add(leaves, leaf_index, leaf->p[i], leaf->index[i] );
+            for (auto i : order.Range(isplit, N))
+                leaf2->Add(leaves, 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( leaves, leaf_index, p, pi );
+            else
+                leaf2->Add( leaves, leaf_index, p, pi );
+
+            ball_leaves.Free(leaf);
+            n_leaves++;
+            n_nodes+=2;
+          }
+      }
+
+    void DeleteElement (T pi)
+      {
+        // static Timer timer("DelaunayTree::DeleteElement"); RegionTimer rt(timer);
+        Leaf *leaf = leaves[leaf_index[pi]];
+        leaf_index[pi] = -1;
+        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 BoxTree<3> DTREE;
+  typedef DelaunayTree<3> DTREE;
 
   static const int deltetfaces[][3] = 
     { { 1, 2, 3 },
@@ -227,14 +507,14 @@ namespace netgen
   void AddDelaunayPoint (PointIndex newpi, const Point3d & newp, 
 			 NgArray<DelaunayTet> & tempels, 
 			 Mesh & mesh,
-			 BoxTree<3> & tettree,
+			 DTREE & tettree,
 			 MeshNB & meshnb,
 			 NgArray<Point<3> > & centers, NgArray<double> & radi2,
 			 NgArray<int> & connected, NgArray<int> & treesearch, 
 			 NgArray<int> & freelist, SphereList & list,
 			 IndexSet & insphere, IndexSet & closesphere)
   {
-    static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t);
+    static Timer t("Meshing3::AddDelaunayPoint");// RegionTimer reg(t);
     // static Timer tsearch("addpoint, search");
     // static Timer tinsert("addpoint, insert");
 
@@ -635,7 +915,7 @@ namespace netgen
     pmin2 = pmin2 + 0.1 * (pmin2 - pmax2);
     pmax2 = pmax2 + 0.1 * (pmax2 - pmin2);
 
-    BoxTree<3> tettree(pmin2, pmax2);
+    DTREE tettree(pmin2, pmax2);
 
 
     tempels.Append (startel);