ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-24 21:10:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

cleanup, more parallel

Browse Source
This commit is contained in:
Matthias Hochsteger 2021-06-04 16:23:18 +02:00
parent 36d9ead3bc
commit ba148e8b3b
1 changed files with 70 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

134
libsrc/meshing/delaunay.cpp
View File

@ -737,10 +737,13 @@ namespace netgen
} }
void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray<DelaunayTet> & tempels, NgBitArray & badnode, int np ) void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray<DelaunayTet> & tempels, int np )
{ {
static Timer tdegenerated("Delaunay - remove degenerated"); RegionTimer rt(tdegenerated); static Timer tdegenerated("Delaunay - remove degenerated"); RegionTimer rt(tdegenerated);
NgBitArray badnode(points.Size());
badnode.Clear();
int ndeg = 0; int ndeg = 0;
for (int i = 1; i <= tempels.Size(); i++) for (int i = 1; i <= tempels.Size(); i++)
{ {
@ -791,7 +794,7 @@ namespace netgen
} }
// Remove flat tets containing two adjacent surface trigs // Remove flat tets containing two adjacent surface trigs
void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray<DelaunayTet> & tempels, NgArray<int> & openels, NgBitArray & badnode, int np ) void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray<DelaunayTet> & tempels, NgArray<int> & openels )
{ {
static Timer topenel("Delaunay - find openel"); RegionTimer rt(topenel); static Timer topenel("Delaunay - find openel"); RegionTimer rt(topenel);
@ -876,76 +879,82 @@ namespace netgen
table.Add(tri[2], openel_i); table.Add(tri[2], openel_i);
}, mesh.GetNP()); }, mesh.GetNP());
for (auto i : openels) ngcore::BitArray badnode(mesh.GetNP()+PointIndex::BASE);
{ badnode.Clear();
const Element2d & tri = mesh.OpenElement(i);
for( auto edge : Range(3) ) ngcore::ParallelForRange(openels.Size(), [&] (auto myrange) {
{ for (auto i_ : myrange)
auto pi0 = tri[edge]; {
auto pi1 = tri[(edge+1)%3]; auto i = openels[i_];
if(pi0>pi1) const Element2d & tri = mesh.OpenElement(i);
Swap(pi0, pi1);
// find other trig with edge pi0, pi1 for( auto edge : Range(3) )
int i_other = -1;
for(auto ii : p2sel[pi0])
{ {
if(ii==i) auto pi0 = tri[edge];
continue; auto pi1 = tri[(edge+1)%3];
auto & tri_other = mesh.OpenElement(ii); if(pi0>pi1)
if(tri_other[0]==pi1 || tri_other[1]==pi1 || tri_other[2]==pi1) Swap(pi0, pi1);
// find other trig with edge pi0, pi1
int i_other = -1;
for(auto ii : p2sel[pi0])
{ {
i_other = ii; if(ii==i)
break; continue;
} auto & tri_other = mesh.OpenElement(ii);
} if(tri_other[0]==pi1 || tri_other[1]==pi1 || tri_other[2]==pi1)
if(i_other>i)
{
auto & tri_other = mesh.OpenElement(i_other);
PointIndex pi2 = tri[(edge+2)%3];
PointIndex pi3 = tri_other[0]+tri_other[1]+tri_other[2] - pi0 - pi1;
if(pi2>pi3)
Swap(pi2, pi3);
// search for tet with edge pi2-pi3
for(auto ei : p2el[pi2])
{
auto & el = tempels[ei];
if(el[0]==pi3 || el[1]==pi3 || el[2]==pi3 || el[3]==pi3)
{ {
const Point3d & p1 = mesh[pi0]; i_other = ii;
const Point3d & p2 = mesh[pi1];
const Point3d & p3 = mesh[pi2];
const Point3d & p4 = mesh[pi3];
Vec3d v1(p1, p2);
Vec3d v2(p1, p3);
Vec3d v3(p1, p4);
Vec3d n = Cross (v1, v2);
double vol = n * v3;
double h = v1.Length() + v2.Length() + v3.Length();
if (fabs (vol) < 1e-4 * (h * h * h)) // old: 1e-12
{
badnode.Set(pi2);
badnode.Set(pi3);
}
break; break;
} }
} }
if(i_other>i)
{
auto & tri_other = mesh.OpenElement(i_other);
PointIndex pi2 = tri[(edge+2)%3];
PointIndex pi3 = tri_other[0]+tri_other[1]+tri_other[2] - pi0 - pi1;
if(pi2>pi3)
Swap(pi2, pi3);
// search for tet with edge pi2-pi3
for(auto ei : p2el[pi2])
{
auto & el = tempels[ei];
if(el[0]==pi3 || el[1]==pi3 || el[2]==pi3 || el[3]==pi3)
{
const Point3d & p1 = mesh[pi0];
const Point3d & p2 = mesh[pi1];
const Point3d & p3 = mesh[pi2];
const Point3d & p4 = mesh[pi3];
Vec3d v1(p1, p2);
Vec3d v2(p1, p3);
Vec3d v3(p1, p4);
Vec3d n = Cross (v1, v2);
double vol = n * v3;
double h = v1.Length() + v2.Length() + v3.Length();
if (fabs (vol) < 1e-4 * (h * h * h)) // old: 1e-12
{
badnode.SetBitAtomic(pi2);
badnode.SetBitAtomic(pi3);
}
break;
}
}
}
} }
} }
} });
for (int i = ne; i >= 1; i--) for (int i = ne; i >= 1; i--)
{ {
const DelaunayTet & el = tempels.Get(i); const DelaunayTet & el = tempels.Get(i);
if (badnode.Test(el[0]) || if (badnode[el[0]] ||
badnode.Test(el[1]) || badnode[el[1]] ||
badnode.Test(el[2]) || badnode[el[2]] ||
badnode.Test(el[3]) ) badnode[el[3]] )
tempels.DeleteElement(i); tempels.DeleteElement(i);
} }
} }
@ -1643,13 +1652,10 @@ namespace netgen
tempels.Append (el); tempels.Append (el);
} }
NgBitArray badnode(mesh.GetNP()); DelaunayRemoveDegenerated(mesh.Points(), tempels, np);
badnode.Clear();
DelaunayRemoveDegenerated(mesh.Points(), tempels, badnode, np);
NgArray<int> openels; NgArray<int> openels;
DelaunayRemoveTwoTriaTets(mesh, tempels, openels, badnode, np); DelaunayRemoveTwoTriaTets(mesh, tempels, openels);
DelaunayRemoveIntersecting(mesh, tempels, openels, pmin, pmax); DelaunayRemoveIntersecting(mesh, tempels, openels, pmin, pmax);
DelaunayRemoveOuter(mesh, tempels, adfront); DelaunayRemoveOuter(mesh, tempels, adfront);