From b2808109b640dd8184c4de4a0e3e7713d2e0f669 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Wed, 25 Sep 2019 15:04:27 +0200 Subject: [PATCH 1/2] New BoxTree for Delaunay --- libsrc/meshing/delaunay.cpp | 227 +++++++++++++++++++++++++++++++++++- 1 file changed, 224 insertions(+), 3 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 9ebacaca..8842150f 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -2,10 +2,231 @@ #include "meshing.hpp" - namespace netgen { + template + class BTree + { + // 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, 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 + { + double sep; + Node *left, *right; + Leaf *leaf; + int level; + + Node() = default; + ~Node() + { + delete leaf; + delete left; + delete right; + } + }; + + Node root; + Array leaf_index; + Point global_min, global_max; + double tol; + public: + + BTree (const Point & pmin, const Point & pmax) + : global_min(pmin), global_max(pmax) + { + root.leaf = new Leaf(); + root.level = 0; + tol = 1e-7 * Dist(pmax, pmin); + } + + void GetIntersecting (const Point & pmin, const Point & pmax, + NgArray & pis) const + { + static Array stack(1000); + static Array dir_stack(1000); + static NgArray ref_pis; + + pis.SetSize(0); + + 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->leaf) + { + 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]) + continue; + for (int d = dim; d < 2*dim; d++) + if (p[d] < tpmin[d]) + continue; + pis.Append (leaf->index[i]); + } + } + else + { + int newdir = dir+1; + if(newdir==2*dim) newdir = 0; + if (tpmin[dir] <= node->sep) + { + stack.Append(node->left); + dir_stack.Append(newdir); + } + if (tpmax[dir] >= node->sep) + { + stack.Append(node->right); + dir_stack.Append(newdir); + } + } + } + } + + void Insert (const Point & pmin, const Point & pmax, T pi) + { + 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->leaf; + + // search correct leaf to add point + while(!leaf) + { + node = p[dir] < node->sep ? node->left : node->right; + dir++; + if(dir==2*dim) dir = 0; + leaf = node->leaf; + } + + // 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 coords(n_elements); + ArrayMem 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 = new Leaf(); + Leaf *leaf2 = new 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 = new Node(); + node1->leaf = leaf1; + node1->level = node->level+1; + + Node *node2 = new Node(); + node2->leaf = leaf2; + node2->level = node->level+1; + + node->left = node1; + node->right = node2; + node->leaf = nullptr; + 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 ); + + delete leaf; + } + } + + void DeleteElement (T pi) + { + 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 }, @@ -214,7 +435,7 @@ namespace netgen void AddDelaunayPoint (PointIndex newpi, const Point3d & newp, NgArray & tempels, Mesh & mesh, - BoxTree<3> & tettree, + TBoxTree & tettree, MeshNB & meshnb, NgArray > & centers, NgArray & radi2, NgArray & connected, NgArray & treesearch, @@ -586,7 +807,7 @@ namespace netgen pmin2 = pmin2 + 0.1 * (pmin2 - pmax2); pmax2 = pmax2 + 0.1 * (pmax2 - pmin2); - BoxTree<3> tettree(pmin2, pmax2); + TBoxTree tettree(pmin2, pmax2); tempels.Append (startel); From 7b1c05f12c065c2fa0bdf15f37444aa24ec16799 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Thu, 26 Sep 2019 14:02:37 +0200 Subject: [PATCH 2/2] Save memory in BTree, count number of leaves and nodes --- libsrc/meshing/delaunay.cpp | 65 +++++++++++++++++++++++++++---------- 1 file changed, 48 insertions(+), 17 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 8842150f..0ff41dc7 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -8,6 +8,7 @@ namespace netgen template class BTree { + public: // Number of entries per leaf static constexpr int N = 100; @@ -34,40 +35,67 @@ namespace netgen struct Node { + union + { + Node *children[2]; + Leaf *leaf; + }; double sep; - Node *left, *right; - Leaf *leaf; int level; - Node() = default; + Node() + : children{nullptr,nullptr} + { } + ~Node() { - delete leaf; - delete left; - delete right; + if(children[1]) + { + delete children[0]; + delete children[1]; + } + else + delete leaf; } + + Leaf *GetLeaf() const + { + return children[1] ? nullptr : leaf; + } }; + private: Node root; Array leaf_index; Point global_min, global_max; double tol; + size_t n_leaves; + size_t n_nodes; public: BTree (const Point & pmin, const Point & pmax) - : global_min(pmin), global_max(pmax) + : global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1) { root.leaf = new Leaf(); root.level = 0; tol = 1e-7 * Dist(pmax, pmin); } + size_t GetNLeaves() + { + return n_leaves; + } + + size_t GetNNodes() + { + return n_nodes; + } + void GetIntersecting (const Point & pmin, const Point & pmax, NgArray & pis) const { static Array stack(1000); static Array dir_stack(1000); - static NgArray ref_pis; pis.SetSize(0); @@ -95,7 +123,7 @@ namespace netgen int dir = dir_stack.Last(); dir_stack.DeleteLast(); - if(Leaf *leaf=node->leaf) + if(Leaf *leaf = node->GetLeaf()) { for (auto i : IntRange(leaf->n_elements)) { @@ -117,12 +145,12 @@ namespace netgen if(newdir==2*dim) newdir = 0; if (tpmin[dir] <= node->sep) { - stack.Append(node->left); + stack.Append(node->children[0]); dir_stack.Append(newdir); } if (tpmax[dir] >= node->sep) { - stack.Append(node->right); + stack.Append(node->children[1]); dir_stack.Append(newdir); } } @@ -140,15 +168,15 @@ namespace netgen } Node * node = &root; - Leaf * leaf = node->leaf; + Leaf * leaf = node->GetLeaf(); // search correct leaf to add point while(!leaf) { - node = p[dir] < node->sep ? node->left : node->right; + node = p[dir] < node->sep ? node->children[0] : node->children[1]; dir++; if(dir==2*dim) dir = 0; - leaf = node->leaf; + leaf = node->GetLeaf(); } // add point to leaf @@ -186,9 +214,8 @@ namespace netgen node2->leaf = leaf2; node2->level = node->level+1; - node->left = node1; - node->right = node2; - node->leaf = nullptr; + 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 @@ -198,6 +225,8 @@ namespace netgen leaf2->Add( leaf_index, p, pi ); delete leaf; + n_leaves++; + n_nodes+=2; } } @@ -891,6 +920,8 @@ namespace netgen PrintMessage (3, "Points: ", cntp); PrintMessage (3, "Elements: ", tempels.Size()); + PrintMessage (3, "Tree data entries per element: ", 1.0*tettree.N*tettree.GetNLeaves() / tempels.Size()); + PrintMessage (3, "Tree nodes per element: ", 1.0*tettree.GetNNodes() / tempels.Size()); // (*mycout) << cntp << " / " << tempels.Size() << " points/elements" << endl; /*