ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-04-09 06:47:28 +05:00
Code Issues Packages Projects Releases Wiki Activity

LocalH::Find() utility function to find GradingBox

Browse Source 
Simplifies code and avoids searching for same grading box twice in SetH().
This commit is contained in:
Matthias Hochsteger 2025-03-26 15:23:18 +01:00
parent 3b79dbc8ff
commit f15ba64a90
2 changed files with 71 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

159
libsrc/meshing/localh.cpp
View File

@ -196,54 +196,41 @@ namespace netgen
fabs (p(1) - root->xmid[1]) > root->h2) fabs (p(1) - root->xmid[1]) > root->h2)
return; return;
if (GetH(p) <= 1.2 * h) return; GradingBox * box = Find(p);
if (box->HOpt() <= 1.2 * h) return;
GradingBox * box = root;
GradingBox * nbox = root;
GradingBox * ngb;
int childnr;
double x1[3], x2[3];
while (nbox)
{
box = nbox;
childnr = 0;
if (p(0) > box->xmid[0]) childnr += 1;
if (p(1) > box->xmid[1]) childnr += 2;
nbox = box->childs[childnr];
};
while (2 * box->h2 > h) while (2 * box->h2 > h)
{ {
childnr = 0; int childnr = 0;
if (p(0) > box->xmid[0]) childnr += 1; double x1[3], x2[3];
if (p(1) > box->xmid[1]) childnr += 2;
double h2 = box->h2; double h2 = box->h2;
if (childnr & 1) if (p(0) > box->xmid[0])
{ {
childnr += 1;
x1[0] = box->xmid[0]; x1[0] = box->xmid[0];
x2[0] = x1[0]+h2; // box->x2[0]; x2[0] = x1[0]+h2;
} }
else else
{ {
x2[0] = box->xmid[0]; x2[0] = box->xmid[0];
x1[0] = x2[0]-h2; // box->x1[0]; x1[0] = x2[0]-h2;
} }
if (childnr & 2) if (p(1) > box->xmid[1])
{ {
childnr += 2;
x1[1] = box->xmid[1]; x1[1] = box->xmid[1];
x2[1] = x1[1]+h2; // box->x2[1]; x2[1] = x1[1]+h2;
} }
else else
{ {
x2[1] = box->xmid[1]; x2[1] = box->xmid[1];
x1[1] = x2[1]-h2; // box->x1[1]; x1[1] = x2[1]-h2;
} }
x1[2] = x2[2] = 0; x1[2] = x2[2] = 0;
ngb = new GradingBox (x1, x2); auto ngb = new GradingBox (x1, x2);
box->childs[childnr] = ngb; box->childs[childnr] = ngb;
ngb->father = box; ngb->father = box;
@ -276,67 +263,52 @@ namespace netgen
fabs (p(2) - root->xmid[2]) > root->h2) fabs (p(2) - root->xmid[2]) > root->h2)
return; return;
if (GetH(p) <= 1.2 * h) return; GradingBox * box = Find(p);
if (box->HOpt() <= 1.2 * h) return;
GradingBox * box = root;
GradingBox * nbox = root;
GradingBox * ngb;
int childnr;
double x1[3], x2[3];
while (nbox)
{
box = nbox;
childnr = 0;
if (p(0) > box->xmid[0]) childnr += 1;
if (p(1) > box->xmid[1]) childnr += 2;
if (p(2) > box->xmid[2]) childnr += 4;
nbox = box->childs[childnr];
};
while (2 * box->h2 > h) while (2 * box->h2 > h)
{ {
childnr = 0; double x1[3], x2[3];
if (p(0) > box->xmid[0]) childnr += 1; int childnr = 0;
if (p(1) > box->xmid[1]) childnr += 2;
if (p(2) > box->xmid[2]) childnr += 4;
double h2 = box->h2; double h2 = box->h2;
if (childnr & 1)
if (p(0) > box->xmid[0])
{ {
childnr += 1;
x1[0] = box->xmid[0]; x1[0] = box->xmid[0];
x2[0] = x1[0]+h2; // box->x2[0]; x2[0] = x1[0]+h2;
} }
else else
{ {
x2[0] = box->xmid[0]; x2[0] = box->xmid[0];
x1[0] = x2[0]-h2; // box->x1[0]; x1[0] = x2[0]-h2;
} }
if (childnr & 2) if (p(1) > box->xmid[1])
{ {
childnr += 2;
x1[1] = box->xmid[1]; x1[1] = box->xmid[1];
x2[1] = x1[1]+h2; // box->x2[1]; x2[1] = x1[1]+h2;
} }
else else
{ {
x2[1] = box->xmid[1]; x2[1] = box->xmid[1];
x1[1] = x2[1]-h2; // box->x1[1]; x1[1] = x2[1]-h2;
} }
if (childnr & 4) if (p(2) > box->xmid[2])
{ {
childnr += 4;
x1[2] = box->xmid[2]; x1[2] = box->xmid[2];
x2[2] = x1[2]+h2; // box->x2[2]; x2[2] = x1[2]+h2;
} }
else else
{ {
x2[2] = box->xmid[2]; x2[2] = box->xmid[2];
x1[2] = x2[2]-h2; // box->x1[2]; x1[2] = x2[2]-h2;
} }
ngb = new GradingBox (x1, x2); auto ngb = new GradingBox (x1, x2);
box->childs[childnr] = ngb; box->childs[childnr] = ngb;
ngb->father = box; ngb->father = box;
@ -366,37 +338,7 @@ namespace netgen
double LocalH :: GetH (Point<3> x) const double LocalH :: GetH (Point<3> x) const
{ {
const GradingBox * box = root; return Find(x)->HOpt();
if (dimension == 2)
{
while (1)
{
int childnr = 0;
if (x(0) > box->xmid[0]) childnr += 1;
if (x(1) > box->xmid[1]) childnr += 2;
if (box->childs[childnr])
box = box->childs[childnr];
else
return box->hopt;
}
}
else
{
while (1)
{
int childnr = 0;
if (x(0) > box->xmid[0]) childnr += 1;
if (x(1) > box->xmid[1]) childnr += 2;
if (x(2) > box->xmid[2]) childnr += 4;
if (box->childs[childnr])
box = box->childs[childnr];
else
return box->hopt;
}
}
} }
@ -488,6 +430,41 @@ namespace netgen
} }
GradingBox * LocalH :: Find (Point<3> p) const
{
GradingBox * box = root;
if (dimension == 2)
{
while (1)
{
int childnr = 0;
if (p(0) > box->xmid[0]) childnr += 1;
if (p(1) > box->xmid[1]) childnr += 2;
if (box->childs[childnr])
box = box->childs[childnr];
else
return box;
}
}
else
{
while (1)
{
int childnr = 0;
if (p(0) > box->xmid[0]) childnr += 1;
if (p(1) > box->xmid[1]) childnr += 2;
if (p(2) > box->xmid[2]) childnr += 4;
if (box->childs[childnr])
box = box->childs[childnr];
else
return box;
}
}
return nullptr;
}
void LocalH :: FindInnerBoxes (const AdFront3 & adfront, void LocalH :: FindInnerBoxes (const AdFront3 & adfront,
int (*testinner)(const Point3d & p1)) int (*testinner)(const Point3d & p1))

3
libsrc/meshing/localh.hpp
View File

@ -54,6 +54,7 @@ namespace netgen
Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); } Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); }
double H2() const { return h2; } double H2() const { return h2; }
double HOpt() const { return hopt; }
bool HasChilds() const bool HasChilds() const
{ {
@ -119,6 +120,8 @@ namespace netgen
void CutBoundary (const Box<3> & box) void CutBoundary (const Box<3> & box)
{ CutBoundaryRec (box.PMin(), box.PMax(), root); } { CutBoundaryRec (box.PMin(), box.PMax(), root); }
GradingBox * Find(Point<3> p) const;
/// find inner boxes /// find inner boxes
void FindInnerBoxes (const class AdFront3 & adfront, void FindInnerBoxes (const class AdFront3 & adfront,