From ef64a5e7eb649924bde86e629cf8ff85da6b5d69 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 1 Oct 2019 11:22:30 +0200 Subject: [PATCH 1/3] Replace BoxTree with BTree --- libsrc/gprim/adtree.hpp | 381 ++++++++++++++++++++++++++++++------ libsrc/meshing/delaunay.cpp | 259 +----------------------- 2 files changed, 325 insertions(+), 315 deletions(-) diff --git a/libsrc/gprim/adtree.hpp b/libsrc/gprim/adtree.hpp index 6ab7669e..7e3686bc 100644 --- a/libsrc/gprim/adtree.hpp +++ b/libsrc/gprim/adtree.hpp @@ -741,90 +741,355 @@ public: const ADTree3 & Tree() const { return *tree; }; }; - - template + template class BoxTree { - T_ADTree<2*dim,T> * tree; - Point boxpmin, boxpmax; + 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( ClosedHashTable &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 + { + 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; + ClosedHashTable leaf_index; + + Point global_min, global_max; + double tol; + size_t n_leaves; + size_t n_nodes; + BlockAllocator ball_nodes; + BlockAllocator ball_leaves; + public: BoxTree (const Box & abox) + : BoxTree( abox.PMin(), abox.PMax() ) + { } + + BoxTree (const Point & pmin, const Point & pmax) + : global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf)) { - boxpmin = abox.PMin(); - boxpmax = abox.PMax(); + 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 + void GetFirstIntersecting (const Point & pmin, const Point & pmax, + TFunc func=[](auto pi){return false;}) const + { + // static Timer timer("BTree::GetIntersecting"); RegionTimer rt(timer); + // static Timer timer1("BTree::GetIntersecting-LinearSearch"); + ArrayMem stack; + ArrayMem 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); - tpmax(i) = tpmax(i+dim) = boxpmax(i); + tpmin(i) = global_min(i); + tpmax(i) = pmax(i)+tol; + + tpmin(i+dim) = pmin(i)-tol; + tpmax(i+dim) = global_max(i); } - tree = new T_ADTree<2*dim,T> (tpmin, tpmax); - } - - BoxTree (const Point & apmin, const Point & apmax) - { - boxpmin = apmin; - boxpmax = apmax; - Point<2*dim> tpmin, tpmax; - for (int i = 0; i < dim; i++) + + stack.SetSize(0); + stack.Append(&root); + dir_stack.SetSize(0); + dir_stack.Append(0); + + while(stack.Size()) { - tpmin(i) = tpmin(i+dim) = boxpmin(i); - tpmax(i) = tpmax(i+dim) = boxpmax(i); + 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); + } + } } - tree = new T_ADTree<2*dim,T> (tpmin, tpmax); } - - ~BoxTree () + + void GetIntersecting (const Point & pmin, const Point & pmax, + NgArray & pis) const { - delete tree; + pis.SetSize(0); + GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;}); } - - void Insert (const Point & bmin, const Point & 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 & box, T pi) { Insert (box.PMin(), box.PMax(), pi); } - - void DeleteElement (T pi) + + void Insert (const Point & pmin, const Point & pmax, T pi) { - tree->DeleteElement(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 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 = (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 GetIntersecting (const Point & pmin, const Point & pmax, - NgArray & pis) const + void DeleteElement (T pi) { - Point<2*dim> tpmin, tpmax; - double tol = Tolerance(); - for (size_t i = 0; i < dim; i++) + // static Timer timer("BTree::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; + + for (auto i : IntRange(n_elements)) { - tpmin(i) = boxpmin(i); - tpmax(i) = pmax(i)+tol; - - tpmin(i+dim) = pmin(i)-tol; - tpmax(i+dim) = boxpmax(i); + if(index[i] == pi) + { + n_elements--; + if(i!=n_elements) + { + index[i] = index[n_elements]; + p[i] = p[n_elements]; + } + return; + } } - - 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; }; }; +// template +// class BoxTree +// { +// T_ADTree<2*dim,T> * tree; +// Point boxpmin, boxpmax; +// public: +// BoxTree (const Box & 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 & apmin, const Point & 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 & bmin, const Point & 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 & box, T pi) +// { +// Insert (box.PMin(), box.PMax(), pi); +// } +// +// void DeleteElement (T pi) +// { +// tree->DeleteElement(pi); +// } +// +// void GetIntersecting (const Point & pmin, const Point & pmax, +// NgArray & 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; }; +// }; + } #endif diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 2de53c43..3cb18261 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -5,264 +5,9 @@ namespace netgen { - template - 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, 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; - Point global_min, global_max; - double tol; - size_t n_leaves; - size_t n_nodes; - BlockAllocator ball_nodes; - BlockAllocator ball_leaves; - - public: - - BTree (const Point & pmin, const Point & 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 - void GetFirstIntersecting (const Point & pmin, const Point & pmax, - TFunc func=[](auto pi){return false;}) const - { - // static Timer timer("BTree::GetIntersecting"); RegionTimer rt(timer); - // static Timer timer1("BTree::GetIntersecting-LinearSearch"); - static Array stack(1000); - static Array 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 & pmin, const Point & pmax, - NgArray & pis) const - { - pis.SetSize(0); - GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;}); - } - - void Insert (const Point & pmin, const Point & 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 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 = (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; +// typedef BTree<3> TBoxTree; + typedef BoxTree<3> TBoxTree; static const int deltetfaces[][3] = { { 1, 2, 3 }, From 7ced41e56f7ea57b0a3385e19e5762cf370aede7 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 1 Oct 2019 11:23:22 +0200 Subject: [PATCH 2/3] Use searchtree in STLBoundary --- libsrc/gprim/adtree.hpp | 173 ++++++++++++++++---------------- libsrc/meshing/delaunay.cpp | 8 +- libsrc/stlgeom/stlgeomchart.cpp | 3 - libsrc/stlgeom/stltool.cpp | 12 +-- libsrc/stlgeom/stltool.hpp | 4 +- 5 files changed, 97 insertions(+), 103 deletions(-) diff --git a/libsrc/gprim/adtree.hpp b/libsrc/gprim/adtree.hpp index 7e3686bc..11db1d78 100644 --- a/libsrc/gprim/adtree.hpp +++ b/libsrc/gprim/adtree.hpp @@ -741,97 +741,98 @@ public: const ADTree3 & Tree() const { return *tree; }; }; - template - class BoxTree +template +class BoxTree +{ +public: + // Number of entries per leaf + static constexpr int N = 100; + + struct Node; + + struct Leaf { - public: - // Number of entries per leaf - static constexpr int N = 100; + Point<2*dim> p[N]; + T index[N]; + int n_elements; - struct Node; + Leaf() : n_elements(0) + { } - struct Leaf - { - Point<2*dim> p[N]; - T index[N]; - int n_elements; - - Leaf() : n_elements(0) {} - - void Add( ClosedHashTable &leaf_index, const Point<2*dim> &ap, T aindex ) + void Add( ClosedHashTable &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 + struct Node + { + union { - union - { - Node *children[2]; - Leaf *leaf; - }; - double sep; - int level; - - Node() - : children{nullptr,nullptr} - { } - - ~Node() - { } - - Leaf *GetLeaf() const - { - return children[1] ? nullptr : leaf; - } + Node *children[2]; + Leaf *leaf; }; + double sep; + int level; - private: - Node root; - -// Array leaf_index; - ClosedHashTable leaf_index; - - Point global_min, global_max; - double tol; - size_t n_leaves; - size_t n_nodes; - BlockAllocator ball_nodes; - BlockAllocator ball_leaves; - - public: - BoxTree (const Box & abox) - : BoxTree( abox.PMin(), abox.PMax() ) + Node() + : children{nullptr,nullptr} { } - BoxTree (const Point & pmin, const Point & pmax) - : global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf)) + ~Node() + { } + + Leaf *GetLeaf() const + { + return children[1] ? nullptr : leaf; + } + }; + +private: + Node root; + + ClosedHashTable leaf_index; + + Point global_min, global_max; + double tol; + size_t n_leaves; + size_t n_nodes; + BlockAllocator ball_nodes; + BlockAllocator ball_leaves; + +public: + + BoxTree (const Point & pmin, const Point & 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() + BoxTree (const Box & box) + : BoxTree(box.PMin(), box.PMax()) + { } + + size_t GetNLeaves() { return n_leaves; } - size_t GetNNodes() + size_t GetNNodes() { return n_nodes; } - template - void GetFirstIntersecting (const Point & pmin, const Point & pmax, - TFunc func=[](auto pi){return false;}) const + template + void GetFirstIntersecting (const Point & pmin, const Point & pmax, + TFunc func=[](auto pi){return false;}) const { - // static Timer timer("BTree::GetIntersecting"); RegionTimer rt(timer); - // static Timer timer1("BTree::GetIntersecting-LinearSearch"); + // static Timer timer("BoxTree::GetIntersecting"); RegionTimer rt(timer); + // static Timer timer1("BoxTree::GetIntersecting-LinearSearch"); ArrayMem stack; ArrayMem dir_stack; @@ -842,7 +843,7 @@ public: { tpmin(i) = global_min(i); tpmax(i) = pmax(i)+tol; - + tpmin(i+dim) = pmin(i)-tol; tpmax(i+dim) = global_max(i); } @@ -862,7 +863,7 @@ public: if(Leaf *leaf = node->GetLeaf()) { -// RegionTimer rt1(timer1); + // RegionTimer rt1(timer1); for (auto i : IntRange(leaf->n_elements)) { bool intersect = true; @@ -875,7 +876,7 @@ public: if (p[d] < tpmin[d]) intersect = false; if(intersect) - if(func(leaf->index[i])) return; + if(func(leaf->index[i])) return; } } else @@ -896,21 +897,21 @@ public: } } - void GetIntersecting (const Point & pmin, const Point & pmax, - NgArray & pis) const + void GetIntersecting (const Point & pmin, const Point & pmax, + NgArray & pis) const { pis.SetSize(0); GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;}); } - void Insert (const Box & box, T pi) + void Insert (const Box & box, T pi) { Insert (box.PMin(), box.PMax(), pi); } - void Insert (const Point & pmin, const Point & pmax, T pi) + void Insert (const Point & pmin, const Point & pmax, T pi) { - // static Timer timer("BTree::Insert"); RegionTimer rt(timer); + // static Timer timer("BoxTree::Insert"); RegionTimer rt(timer); int dir = 0; Point<2*dim> p; for (auto i : IntRange(dim)) @@ -982,9 +983,9 @@ public: } } - void DeleteElement (T pi) + void DeleteElement (T pi) { - // static Timer timer("BTree::DeleteElement"); RegionTimer rt(timer); + // static Timer timer("BoxTree::DeleteElement"); RegionTimer rt(timer); Leaf *leaf = leaf_index[pi]; leaf_index.Delete(pi); auto & n_elements = leaf->n_elements; @@ -1005,7 +1006,7 @@ public: } } } - }; +}; // template // class BoxTree @@ -1025,7 +1026,7 @@ public: // } // tree = new T_ADTree<2*dim,T> (tpmin, tpmax); // } -// +// // BoxTree (const Point & apmin, const Point & apmax) // { // boxpmin = apmin; @@ -1038,36 +1039,36 @@ public: // } // tree = new T_ADTree<2*dim,T> (tpmin, tpmax); // } -// +// // ~BoxTree () // { // delete tree; // } -// +// // void Insert (const Point & bmin, const Point & 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 & box, T pi) // { // Insert (box.PMin(), box.PMax(), pi); // } -// -// void DeleteElement (T pi) +// +// void DeleteElement (T pi) // { // tree->DeleteElement(pi); // } -// -// void GetIntersecting (const Point & pmin, const Point & pmax, +// +// void GetIntersecting (const Point & pmin, const Point & pmax, // NgArray & pis) const // { // Point<2*dim> tpmin, tpmax; @@ -1076,15 +1077,15 @@ public: // { // 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; }; diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 3cb18261..67c455d2 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -5,10 +5,6 @@ namespace netgen { - -// typedef BTree<3> TBoxTree; - typedef BoxTree<3> TBoxTree; - static const int deltetfaces[][3] = { { 1, 2, 3 }, { 2, 0, 3 }, @@ -231,7 +227,7 @@ namespace netgen void AddDelaunayPoint (PointIndex newpi, const Point3d & newp, NgArray & tempels, Mesh & mesh, - TBoxTree & tettree, + BoxTree<3> & tettree, MeshNB & meshnb, NgArray > & centers, NgArray & radi2, NgArray & connected, NgArray & treesearch, @@ -639,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); diff --git a/libsrc/stlgeom/stlgeomchart.cpp b/libsrc/stlgeom/stlgeomchart.cpp index 972a6fa5..e4cf5b21 100644 --- a/libsrc/stlgeom/stlgeomchart.cpp +++ b/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); diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp index f1b931ba..a70df538 100644 --- a/libsrc/stlgeom/stltool.cpp +++ b/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 = new BoxTree<2,INDEX_2> (box2d); - searchtree = nullptr; + searchtree = make_unique> (box2d); +// searchtree = nullptr; } void STLBoundary :: DeleteSearchTree() { - delete searchtree; searchtree = nullptr; } diff --git a/libsrc/stlgeom/stltool.hpp b/libsrc/stlgeom/stltool.hpp index 9312167e..6c9235d8 100644 --- a/libsrc/stlgeom/stltool.hpp +++ b/libsrc/stlgeom/stltool.hpp @@ -191,10 +191,10 @@ private: const STLChart * chart; // NgArray boundary; ClosedHashTable boundary_ht; - BoxTree<2,INDEX_2> * searchtree = nullptr; + unique_ptr> searchtree; public: STLBoundary(STLGeometry * ageometry); - ~STLBoundary() { delete searchtree; } + ~STLBoundary() {} void Clear() { /* boundary.SetSize(0); */ boundary_ht = ClosedHashTable(); } void SetChart (const STLChart * achart) { chart = achart; } From 8bfccdf1df402298c718a01ed6d60dfbf989cbf5 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 1 Oct 2019 12:34:23 +0200 Subject: [PATCH 3/3] Use BoxTree::GetFirstIntersecting --- libsrc/stlgeom/stltool.cpp | 91 +++++++++++++++----------------------- 1 file changed, 36 insertions(+), 55 deletions(-) diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp index a70df538..b5c9f6f2 100644 --- a/libsrc/stlgeom/stltool.cpp +++ b/libsrc/stlgeom/stltool.cpp @@ -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 pis; - searchtree -> GetIntersecting (box2d.PMin(), box2d.PMax(), pis); - for (auto i2 : pis) - { - const STLBoundarySeg & seg = boundary_ht[i2]; - - if (seg.IsSmoothEdge()) continue; - if (!box2d.Intersect (seg.BoundingBox())) continue; - - 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; - } - } + bool has_intersection = false; + searchtree -> GetFirstIntersecting (box2d.PMin(), box2d.PMax(), + [&] (auto i2) NETGEN_LAMBDA_INLINE + { + has_intersection = hasIntersection(i2); + return has_intersection; + }); + return !has_intersection; } - else - { + { for(auto [i2, seg] : boundary_ht) - { - if (seg.IsSmoothEdge()) continue; - if (!box2d.Intersect (seg.BoundingBox())) continue; - - 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; - } - } - + if(hasIntersection(i2)) + return false; + return true; } - - return ok; }