ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-26 21:00:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

use STL-searchtree

Browse Source
This commit is contained in:
Joachim Schöberl 2019-09-30 18:21:01 +02:00
parent dd70e94143
commit 693bc66d49
2 changed files with 53 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
libsrc/meshing/smoothing2.cpp
View File

@ -730,11 +730,11 @@ namespace netgen
} }
static Timer timer("MeshSmoothing 2D"); static Timer timer("MeshSmoothing 2D");
static int timer1 = NgProfiler::CreateTimer ("MeshSmoothing 2D start"); // static int timer1 = NgProfiler::CreateTimer ("MeshSmoothing 2D start");
static int timer2 = NgProfiler::CreateTimer ("MeshSmoothing 2D - BFGS"); // static int timer2 = NgProfiler::CreateTimer ("MeshSmoothing 2D - BFGS");
RegionTimer reg (timer); RegionTimer reg (timer);
NgProfiler::StartTimer (timer1); // NgProfiler::StartTimer (timer1);
CheckMeshApproximation (mesh); CheckMeshApproximation (mesh);
@ -918,12 +918,18 @@ namespace netgen
int cnt = 0; int cnt = 0;
NgProfiler::StopTimer (timer1); // NgProfiler::StopTimer (timer1);
/* /*
for (PointIndex pi = PointIndex::BASE; pi < mesh.GetNP()+PointIndex::BASE; pi++) for (PointIndex pi = PointIndex::BASE; pi < mesh.GetNP()+PointIndex::BASE; pi++)
if (mesh[pi].Type() == SURFACEPOINT) if (mesh[pi].Type() == SURFACEPOINT)
*/ */
static Timer tloop("MeshSmooting 2D - loop");
static Timer tloop1("MeshSmooting 2D - loop 1");
static Timer tloop2("MeshSmooting 2D - loop 2");
static Timer tloop3("MeshSmooting 2D - loop 3");
tloop.Start();
for (int hi = 0; hi < icompress.Size(); hi++) for (int hi = 0; hi < icompress.Size(); hi++)
{ {
PointIndex pi = icompress[hi]; PointIndex pi = icompress[hi];
@ -941,7 +947,9 @@ namespace netgen
// if (elementsonpoint[pi].Size() == 0) continue; // if (elementsonpoint[pi].Size() == 0) continue;
if (elementsonpoint[hi].Size() == 0) continue; if (elementsonpoint[hi].Size() == 0) continue;
tloop1.Start();
ld.sp1 = mesh[pi]; ld.sp1 = mesh[pi];
// Element2d & hel = mesh[elementsonpoint[pi][0]]; // Element2d & hel = mesh[elementsonpoint[pi][0]];
@ -963,7 +971,9 @@ namespace netgen
ld.lochs.SetSize (0); ld.lochs.SetSize (0);
ld.loc_pnts2.SetSize (0); ld.loc_pnts2.SetSize (0);
ld.loc_pnts3.SetSize (0); ld.loc_pnts3.SetSize (0);
tloop1.Stop();
tloop2.Start();
for (int j = 0; j < elementsonpoint[hi].Size(); j++) for (int j = 0; j < elementsonpoint[hi].Size(); j++)
{ {
SurfaceElementIndex sei = elementsonpoint[hi][j]; SurfaceElementIndex sei = elementsonpoint[hi][j];
@ -987,7 +997,9 @@ namespace netgen
ld.lochs.Append (mesh.GetH(pmid)); ld.lochs.Append (mesh.GetH(pmid));
} }
} }
tloop2.Stop();
GetNormalVector (ld.surfi, ld.sp1, ld.gi1, ld.normal); GetNormalVector (ld.surfi, ld.sp1, ld.gi1, ld.normal);
ld.t1 = ld.normal.GetNormal (); ld.t1 = ld.normal.GetNormal ();
ld.t2 = Cross (ld.normal, ld.t1); ld.t2 = Cross (ld.normal, ld.t1);
@ -1014,8 +1026,9 @@ namespace netgen
x = 0; x = 0;
par.typx = 0.3*ld.lochs[0]; par.typx = 0.3*ld.lochs[0];
NgProfiler::StartTimer (timer2); // NgProfiler::StartTimer (timer2);
tloop3.Start();
if (mixed) if (mixed)
{ {
BFGS (x, surfminfj, par, 1e-6); BFGS (x, surfminfj, par, 1e-6);
@ -1024,8 +1037,8 @@ namespace netgen
{ {
BFGS (x, surfminf, par, 1e-6); BFGS (x, surfminf, par, 1e-6);
} }
tloop3.Stop();
NgProfiler::StopTimer (timer2); // NgProfiler::StopTimer (timer2);
Point3d origp = mesh[pi]; Point3d origp = mesh[pi];
int loci = 1; int loci = 1;

31
libsrc/stlgeom/stltool.cpp
View File

@ -828,7 +828,36 @@ STLTrigId STLChart :: ProjectNormal (Point<3> & p3d) const
int nt = GetNT(); int nt = GetNT();
double lamtol = 1e-6; double lamtol = 1e-6;
QuadraticFunction3d quadfun(p3d, GetNormal()); QuadraticFunction3d quadfun(p3d, GetNormal());
if (inner_searchtree)
{
NgArray<int> trigs;
Point<2> p2d = Project2d (p3d);
inner_searchtree->GetIntersecting(p2d, p2d, trigs);
for (STLTrigId i : trigs)
{
auto & trig = geometry->GetTriangle(i);
const Point<3> & c = trig.center;
if (quadfun.Eval(c) > sqr (trig.rad))
continue;
Point<3> p = p3d;
Vec<3> lam;
int err = trig.ProjectInPlain(geometry->GetPoints(), GetNormal(), p, lam);
bool inside = (err == 0 && lam(0) > -lamtol &&
lam(1) > -lamtol && (1-lam(0)-lam(1)) > -lamtol);
if (inside)
{
p3d = p;
return i;
}
}
}
for (int j = 1; j <= nt; j++) for (int j = 1; j <= nt; j++)
{ {
STLTrigId i = GetTrig1(j); STLTrigId i = GetTrig1(j);