mirror of
https://github.com/NGSolve/netgen.git
synced 2025-01-12 22:20:35 +05:00
New BoxTree for Delaunay
This commit is contained in:
Matthias Hochsteger
parent
206aa3a6bd
commit
b2808109b6
@ -2,10 +2,231 @@
|
||||
#include "meshing.hpp"
|
||||
|
||||
|
||||
|
||||
namespace netgen
|
||||
{
|
||||
|
||||
template<int dim, typename T=INDEX>
|
||||
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> &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> leaf_index;
|
||||
Point<dim> global_min, global_max;
|
||||
double tol;
|
||||
public:
|
||||
|
||||
BTree (const Point<dim> & pmin, const Point<dim> & pmax)
|
||||
: global_min(pmin), global_max(pmax)
|
||||
{
|
||||
root.leaf = new Leaf();
|
||||
root.level = 0;
|
||||
tol = 1e-7 * Dist(pmax, pmin);
|
||||
}
|
||||
|
||||
void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
|
||||
NgArray<T> & pis) const
|
||||
{
|
||||
static Array<const Node*> stack(1000);
|
||||
static Array<int> dir_stack(1000);
|
||||
static NgArray<T> 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<dim> & pmin, const Point<dim> & 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<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 = 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<DelaunayTet> & tempels,
|
||||
Mesh & mesh,
|
||||
BoxTree<3> & tettree,
|
||||
TBoxTree & tettree,
|
||||
MeshNB & meshnb,
|
||||
NgArray<Point<3> > & centers, NgArray<double> & radi2,
|
||||
NgArray<int> & connected, NgArray<int> & 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);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user