netgen/libsrc/meshing/delaunay.cpp
Matthias Hochsteger fe37f7d0b2 Emscripten fixes
2023-09-05 17:13:29 +02:00

1675 lines
40 KiB
C++

#include <mystdlib.h>
#include "meshing.hpp"
namespace netgen
{
// typedef BoxTree<3> DTREE;
typedef DelaunayTree<3> DTREE;
static const int deltetfaces[][3] =
{ { 1, 2, 3 },
{ 2, 0, 3 },
{ 0, 1, 3 },
{ 1, 0, 2 } };
class DelaunayTet
{
PointIndex pnums[4];
int nb[4] = {0,0,0,0};
public:
DelaunayTet () = default;
DelaunayTet (const Element & el)
{
for (int i = 0; i < 4; i++)
pnums[i] = el[i];
}
PointIndex & operator[] (int i) { return pnums[i]; }
PointIndex operator[] (int i) const { return pnums[i]; }
int & NB(int i) { return nb[i]; }
int NB(int i) const { return nb[i]; }
int FaceNr (INDEX_3 & face) const // which face nr is it ?
{
for (int i = 0; i < 3; i++)
if (pnums[i] != face.I1() &&
pnums[i] != face.I2() &&
pnums[i] != face.I3())
return i;
return 3;
}
INDEX_3 GetFace (int i) const
{
return INDEX_3 (pnums[deltetfaces[i][0]],
pnums[deltetfaces[i][1]],
pnums[deltetfaces[i][2]]);
}
void GetFace (int i, Element2d & face) const
{
// face.SetType(TRIG);
face[0] = pnums[deltetfaces[i][0]];
face[1] = pnums[deltetfaces[i][1]];
face[2] = pnums[deltetfaces[i][2]];
}
};
/*
Table to maintain neighbour elements
*/
class MeshNB
{
// face nodes -> one element
INDEX_3_CLOSED_HASHTABLE<int> faces;
//
NgArray<DelaunayTet> & tets;
public:
// estimated number of points
MeshNB (NgArray<DelaunayTet> & atets, int np)
: faces(200), tets(atets)
{ ; }
// add element with 4 nodes
void Add (int elnr);
// delete element with 4 nodes
void Delete (int elnr)
{
DelaunayTet & el = tets.Elem(elnr);
for (int i = 0; i < 4; i++)
faces.Set (el.GetFace(i).Sort(), el.NB(i));
}
// get neighbour of element elnr in direction fnr
int GetNB (int elnr, int fnr)
{
return tets.Get(elnr).NB(fnr);
}
//
void ResetFaceHT (int size)
{
faces.SetSize (size);
}
};
void MeshNB :: Add (int elnr)
{
DelaunayTet & el = tets.Elem(elnr);
for (int i = 0; i < 4; i++)
{
INDEX_3 i3 = INDEX_3::Sort (el.GetFace(i));
int posnr;
if (!faces.PositionCreate (i3, posnr))
{
// face already in use
int othertet = faces.GetData (posnr);
el.NB(i) = othertet;
if (othertet)
{
int fnr = tets.Get(othertet).FaceNr (i3);
tets.Elem(othertet).NB(fnr) = elnr;
}
}
else
{
faces.SetData (posnr, elnr);
el.NB(i) = 0;
}
}
}
/*
connected lists of cosphereical elements
*/
class SphereList
{
NgArray<int> links;
public:
SphereList ()
{ ; }
void AddElement (int elnr)
{
if (elnr > links.Size())
links.Append (1);
links.Elem(elnr) = elnr;
}
void DeleteElement (int elnr)
{
links.Elem(elnr) = 0;
}
void ConnectElement (int eli, int toi)
{
links.Elem (eli) = links.Get (toi);
links.Elem (toi) = eli;
}
void GetList (int eli, NgArray<int> & linked) const;
template <typename TFUNC>
void IterateList (int eli, TFUNC func)
{
int pi = eli;
do
{
func(pi);
pi = links.Get(pi);
}
while (pi != eli);
}
};
void SphereList :: GetList (int eli, NgArray<int> & linked) const
{
linked.SetSize (0);
int pi = eli;
do
{
#ifdef DELAUNAY_DEBUG
if (pi <= 0 || pi > links.Size())
{
cerr << "link, error " << endl;
cerr << "pi = " << pi << " linked.s = " << linked.Size() << endl;
exit(1);
}
if (linked.Size() > links.Size())
{
cerr << "links have loop" << endl;
exit(1);
}
#endif
linked.Append (pi);
pi = links.Get(pi);
}
while (pi != eli);
}
void AddDelaunayPoint (PointIndex newpi, const Point3d & newp,
NgArray<DelaunayTet> & tempels,
Mesh & mesh,
DTREE & tettree,
MeshNB & meshnb,
NgArray<Point<3> > & centers, NgArray<double> & radi2,
NgArray<int> & connected, NgArray<int> & treesearch,
NgArray<int> & freelist, SphereList & list,
IndexSet & insphere, IndexSet & closesphere, Array<DelaunayTet> & newels)
{
static Timer t("Meshing3::AddDelaunayPoint", NoTracing, NoTiming); RegionTimer reg(t);
static Timer tsearch("addpoint, search", NoTracing, NoTiming);
static Timer tfind("addpoint, find all tets", NoTracing, NoTiming);
static Timer tnewtets("addpoint, build new tets", NoTracing, NoTiming);
static Timer tinsert("addpoint, insert", NoTracing, NoTiming);
/*
find any sphere, such that newp is contained in
*/
// DelaunayTet el;
int cfelind = -1;
const Point<3> * pp[4];
Point<3> pc;
Point3d tpmin, tpmax;
/*
// stop search if intersecting point is close enough to center
tettree.GetFirstIntersecting (newp, newp, treesearch, [&](const auto pi)
{
double quot = Dist2 (centers.Get(pi), newp);
return (quot < 0.917632 * radi2.Get(pi));
} );
// tettree.GetIntersecting (newp, newp, treesearch);
double quot,minquot(1e20);
for (auto jjj : treesearch)
{
quot = Dist2 (centers.Get(jjj), newp) / radi2.Get(jjj);
if((cfelind == -1 || quot < 0.99*minquot) && quot < 1)
{
minquot = quot;
el = tempels.Get(jjj);
cfelind = jjj;
if(minquot < 0.917632)
break;
}
}
*/
tsearch.Start();
double minquot{1e20};
tettree.GetFirstIntersecting
(newp, newp, [&](const auto pi)
{
double rad2 = radi2.Get(pi);
double d2 = Dist2 (centers.Get(pi), newp); // / radi2.Get(pi);
if (d2 >= rad2) return false;
// if (d2 < 0.917632 * rad2)
if (d2 < 0.99 * rad2)
{
cfelind = pi;
return true;
}
if (cfelind == -1 || d2 < 0.99*minquot*rad2)
{
minquot = d2/rad2;
cfelind = pi;
}
return false;
} );
tsearch.Stop();
if (cfelind == -1)
{
PrintWarning ("Delaunay, point not in any sphere");
return;
}
tfind.Start();
/*
insphere: point is in sphere -> delete element
closesphere: point is close to sphere -> considered for same center
*/
// save overestimate
insphere.SetMaxIndex (2 * tempels.Size() + 5 * mesh.GetNP());
closesphere.SetMaxIndex (2 * tempels.Size() + 5 * mesh.GetNP());
insphere.Clear();
closesphere.Clear();
insphere.Add (cfelind);
bool changed = true;
int nstarti = 1, starti;
while (changed)
{
changed = false;
starti = nstarti;
nstarti = insphere.GetArray().Size()+1;
// if point in sphere, then it is also closesphere
for (int j = starti; j < nstarti; j++)
{
int helind = insphere.GetArray().Get(j);
if (!closesphere.IsIn (helind))
closesphere.Add (helind);
}
// add connected spheres to insphere - list
for (int j = starti; j < nstarti; j++)
{
list.IterateList (insphere.GetArray().Get(j),
[&](int celind)
{
if (tempels.Get(celind)[0] != PointIndex(PointIndex::INVALID) &&
!insphere.IsIn (celind))
{
changed = true;
insphere.Add (celind);
}
});
}
// check neighbour-tets
for (int j = starti; j < nstarti; j++)
for (int k = 0; k < 4; k++)
{
int helind = insphere.GetArray().Get(j);
int nbind = meshnb.GetNB (helind, k);
if (nbind && !insphere.IsIn (nbind) )
{
double d2 = Dist2 (centers.Get(nbind), newp);
if (d2 < radi2.Get(nbind) * (1+1e-8) )
closesphere.Add (nbind);
if (d2 < radi2.Get(nbind) * (1 + 1e-12))
{
// point is in sphere -> remove tet
insphere.Add (nbind);
changed = true;
}
else
{
INDEX_3 i3 = tempels.Get(helind).GetFace (k);
const Point<3> & p1 = mesh[PointIndex (i3.I1())];
const Point<3> & p2 = mesh[PointIndex (i3.I2())];
const Point<3> & p3 = mesh[PointIndex (i3.I3())];
Vec<3> n = Cross (p2-p1, p3-p1);
n /= n.Length();
double dist = n * (newp-p1);
double scal = n * (mesh.Point (tempels.Get(helind)[k])-p1);
if (scal > 0) dist *= -1;
if (dist > -1e-10) // 1e-10
{
insphere.Add (nbind);
changed = true;
}
}
}
}
} // while (changed)
tfind.Stop();
tnewtets.Start();
newels.SetSize(0);
Element2d face(TRIG);
for (int celind : insphere.GetArray())
for (int k = 0; k < 4; k++)
{
int nbind = meshnb.GetNB (celind, k);
if (!nbind || !insphere.IsIn (nbind))
{
tempels.Get (celind).GetFace (k, face);
// Element newel(TET);
DelaunayTet newel;
for (int l = 0; l < 3; l++)
newel[l] = face[l];
newel[3] = newpi;
newels.Append (newel);
#ifdef DEBUG_DELAUNAY
Vec<3> v1 = mesh[face[1]] - mesh[face[0]];
Vec<3> v2 = mesh[face[2]] - mesh[face[0]];
Vec<3> n = Cross (v1, v2);
n.Normalize();
if (n * Vec3d(mesh.Point (face[0]),
mesh.Point (tempels.Get(celind)[k]))
> 0)
n *= -1;
double hval = n * ( newp - mesh[face[0]]);
if (hval > -1e-12)
{
cerr << "vec to outer" << endl;
(*testout) << "vec to outer, hval = " << hval << endl;
(*testout) << "v1 x v2 = " << Cross (v1, v2) << endl;
(*testout) << "facep: "
<< mesh.Point (face[0]) << " "
<< mesh.Point (face[1]) << " "
<< mesh.Point (face[2]) << endl;
}
#endif
}
}
meshnb.ResetFaceHT (10*insphere.GetArray().Size()+1);
for (auto celind : insphere.GetArray())
{
meshnb.Delete (celind);
list.DeleteElement (celind);
for (int k = 0; k < 4; k++)
tempels.Elem(celind)[k] = PointIndex::INVALID;
tettree.DeleteElement (celind);
freelist.Append (celind);
}
bool hasclose = false;
for (int ind : closesphere.GetArray())
{
if (!insphere.IsIn(ind) &&
fabs (Dist2 (centers.Get (ind), newp) - radi2.Get(ind)) < 1e-8 )
hasclose = true;
}
/*
for (int j = 1; j <= newels.Size(); j++)
{
const auto & newel = newels.Get(j);
*/
for (const auto & newel : newels)
{
int nelind;
if (!freelist.Size())
{
tempels.Append (newel);
nelind = tempels.Size();
}
else
{
nelind = freelist.Last();
freelist.DeleteLast();
tempels.Elem(nelind) = newel;
}
meshnb.Add (nelind);
list.AddElement (nelind);
for (int k = 0; k < 4; k++)
pp[k] = &mesh.Point (newel[k]);
if (CalcSphereCenter (&pp[0], pc) )
{
#ifdef DEBUG_DELAUNAY
PrintSysError ("Delaunay: New tet is flat");
(*testout) << "new tet is flat" << endl;
for (int k = 0; k < 4; k++)
(*testout) << newel[k] << " ";
(*testout) << endl;
for (int k = 0; k < 4; k++)
(*testout) << *pp[k-1] << " ";
(*testout) << endl;
#endif
;
}
double r2 = Dist2 (*pp[0], pc);
if (hasclose)
/*
for (int k = 1; k <= closesphere.GetArray().Size(); k++)
{
int csameind = closesphere.GetArray().Get(k);
*/
for (int csameind : closesphere.GetArray())
{
if (!insphere.IsIn(csameind) &&
fabs (r2 - radi2.Get(csameind)) < 1e-10 &&
Dist2 (pc, centers.Get(csameind)) < 1e-20)
{
pc = centers.Get(csameind);
r2 = radi2.Get(csameind);
list.ConnectElement (nelind, csameind);
break;
}
}
if (centers.Size() < nelind)
{
centers.Append (pc);
radi2.Append (r2);
}
else
{
centers.Elem(nelind) = pc;
radi2.Elem(nelind) = r2;
}
closesphere.Add (nelind);
tpmax = tpmin = *pp[0];
for (int k = 1; k <= 3; k++)
{
tpmin.SetToMin (*pp[k]);
tpmax.SetToMax (*pp[k]);
}
tpmax = tpmax + 0.01 * (tpmax - tpmin);
tinsert.Start();
tettree.Insert (tpmin, tpmax, nelind);
tinsert.Stop();
}
tnewtets.Stop();
}
void Delaunay1 (Mesh & mesh, const MeshingParameters & mp, AdFront3 * adfront,
NgArray<DelaunayTet> & tempels,
int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax)
{
static Timer t("Meshing3::Delaunay1"); RegionTimer reg(t);
NgArray<Point<3>> centers;
NgArray<double> radi2;
Box<3> bbox(Box<3>::EMPTY_BOX);
for (auto & face : adfront->Faces())
for (PointIndex pi : face.Face().PNums())
bbox.Add (mesh.Point(pi));
for (PointIndex pi : mesh.LockedPoints())
bbox.Add (mesh.Point (pi));
pmin = bbox.PMin();
pmax = bbox.PMax();
Vec<3> vdiag = pmax-pmin;
// double r1 = vdiag.Length();
double r1 = sqrt (3.0) * max3(vdiag(0), vdiag(1), vdiag(2));
vdiag = Vec<3> (r1, r1, r1);
//double r2;
Point<3> pmin2 = pmin - 8 * vdiag;
Point<3> pmax2 = pmax + 8 * vdiag;
Point<3> cp1(pmin2), cp2(pmax2), cp3(pmax2), cp4(pmax2);
cp2(0) = pmin2(0);
cp3(1) = pmin2(1);
cp4(2) = pmin2(2);
size_t np = mesh.GetNP();
startel[0] = mesh.AddPoint (cp1);
startel[1] = mesh.AddPoint (cp2);
startel[2] = mesh.AddPoint (cp3);
startel[3] = mesh.AddPoint (cp4);
// flag points to use for Delaunay:
Array<bool, PointIndex> usep(np);
usep = false;
for (auto & face : adfront->Faces())
for (PointIndex pi : face.Face().PNums())
usep[pi] = true;
for (size_t i = oldnp + PointIndex::BASE;
i < np + PointIndex::BASE; i++)
usep[i] = true;
for (PointIndex pi : mesh.LockedPoints())
usep[pi] = true;
NgArray<int> freelist;
int cntp = 0;
MeshNB meshnb (tempels, mesh.GetNP() + 5);
SphereList list;
pmin2 = pmin2 + 0.1 * (pmin2 - pmax2);
pmax2 = pmax2 + 0.1 * (pmax2 - pmin2);
DTREE tettree(pmin2, pmax2);
tempels.Append (startel);
meshnb.Add (1);
list.AddElement (1);
NgArray<int> connected, treesearch;
Box<3> tbox(Box<3>::EMPTY_BOX);
for (size_t k = 0; k < 4; k++)
tbox.Add (mesh.Point(startel[k]));
Point<3> tpmin = tbox.PMin();
Point<3> tpmax = tbox.PMax();
tpmax = tpmax + 0.01 * (tpmax - tpmin);
tettree.Insert (tpmin, tpmax, 1);
Point<3> pc;
const Point<3> * pp[4];
for (int k = 0; k < 4; k++)
pp[k] = &mesh.Point (startel[k]);
CalcSphereCenter (&pp[0], pc);
centers.Append (pc);
radi2.Append (Dist2 (*pp[0], pc));
IndexSet insphere(mesh.GetNP());
IndexSet closesphere(mesh.GetNP());
// "random" reordering of points (speeds a factor 3 - 5 !!!)
NgArray<PointIndex, PointIndex::BASE, PointIndex> mixed(np);
// int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
// int prims[] = { 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
int prims[] = { 211, 223, 227, 229, 233, 239, 241, 251, 257, 263 };
int prim;
{
int i = 0;
while (np % prims[i] == 0) i++;
prim = prims[i];
}
// for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++)
for (PointIndex pi : mesh.Points().Range().Modify(0, -4))
mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE );
Array<DelaunayTet> newels;
// for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++)
for (PointIndex pi : mesh.Points().Range().Modify(0, -4))
{
if (pi % 1000 == 0)
{
if (pi % 10000 == 0)
PrintDot ('+');
else
PrintDot ('.');
}
multithread.percent = 100.0 * pi / np;
if (multithread.terminate)
break;
PointIndex newpi = mixed[pi];
if (!usep[newpi])
continue;
cntp++;
const MeshPoint & newp = mesh[newpi];
AddDelaunayPoint (newpi, newp, tempels, mesh,
tettree, meshnb, centers, radi2,
connected, treesearch, freelist, list, insphere, closesphere, newels);
}
for (int i = tempels.Size(); i >= 1; i--)
if (tempels.Get(i)[0] <= 0)
tempels.DeleteElement (i);
PrintDot ('\n');
PrintMessage (3, "Points: ", cntp);
PrintMessage (3, "Elements: ", tempels.Size());
PrintMessage (3, "Tree data entries per element: ", 1.0*tettree.N*tettree.GetNLeaves() / tempels.Size());
PrintMessage (3, "Tree nodes per element: ", 1.0*tettree.GetNNodes() / tempels.Size());
// (*mycout) << cntp << " / " << tempels.Size() << " points/elements" << endl;
/*
cout << "tempels: ";
tempels.PrintMemInfo(cout);
cout << "Searchtree: ";
tettree.Tree().PrintMemInfo(cout);
cout << "MeshNB: ";
meshnb.PrintMemInfo(cout);
*/
}
void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray<DelaunayTet> & tempels, int np )
{
static Timer tdegenerated("Delaunay - remove degenerated"); RegionTimer rt(tdegenerated);
NgBitArray badnode(points.Size());
badnode.Clear();
int ndeg = 0;
for (int i = 1; i <= tempels.Size(); i++)
{
Element el(4);
for (int j = 0; j < 4; j++)
el[j] = tempels.Elem(i)[j];
// Element & el = tempels.Elem(i);
const Point3d & lp1 = points[el[0]];
const Point3d & lp2 = points[el[1]];
const Point3d & lp3 = points[el[2]];
const Point3d & lp4 = points[el[3]];
Vec3d v1(lp1, lp2);
Vec3d v2(lp1, lp3);
Vec3d v3(lp1, lp4);
Vec3d n = Cross (v1, v2);
double vol = n * v3;
double h = v1.Length() + v2.Length() + v3.Length();
if (fabs (vol) < 1e-8 * (h * h * h) &&
(el[0] <= np && el[1] <= np &&
el[2] <= np && el[3] <= np) ) // old: 1e-12
{
badnode.Set(el[0]);
badnode.Set(el[1]);
badnode.Set(el[2]);
badnode.Set(el[3]);
ndeg++;
(*testout) << "vol = " << vol << " h = " << h << endl;
}
if (vol > 0)
Swap (el[2], el[3]);
}
auto ne = tempels.Size();
for (int i = ne; i >= 1; i--)
{
const DelaunayTet & el = tempels.Get(i);
if (badnode.Test(el[0]) ||
badnode.Test(el[1]) ||
badnode.Test(el[2]) ||
badnode.Test(el[3]) )
tempels.DeleteElement(i);
}
PrintMessage (3, ndeg, " degenerated elements removed");
}
// Remove flat tets containing two adjacent surface trigs
void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray<DelaunayTet> & tempels, NgArray<int> & openels )
{
static Timer topenel("Delaunay - find openel"); RegionTimer rt(topenel);
// find surface triangles which are no face of any tet
BitArray bnd_points( mesh.GetNP() + PointIndex::BASE );
bnd_points.Clear();
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
{
const Element2d & tri = mesh.OpenElement(i);
bnd_points.SetBit(tri[0]);
bnd_points.SetBit(tri[1]);
bnd_points.SetBit(tri[2]);
}
auto ne = tempels.Size();
Array<int> tets_with_3_bnd_points(ne);
atomic<int> cnt = 0;
// table of tets with >= 2 boundary points, store in extra array tets with >=3 boundary points
auto p2el = ngcore::CreateSortedTable<int, PointIndex>( ne,
[&](auto & table, int ei)
{
const auto & el = tempels[ei];
int num_bnd_points = 0;
for( auto i : Range(4) )
if(bnd_points[el[i]])
num_bnd_points++;
if(num_bnd_points>1)
{
table.Add (el[0], ei);
table.Add (el[1], ei);
table.Add (el[2], ei);
table.Add (el[3], ei);
}
// table creator is running this code 2 times, only store tets on last run
if(table.GetMode()==3 && num_bnd_points>2)
tets_with_3_bnd_points[cnt++] = ei;
}, mesh.GetNP());
tets_with_3_bnd_points.SetSize(cnt);
static Timer t1("Build face table"); t1.Start();
ngcore::ClosedHashTable< ngcore::INT<3>, int > face_table( 4*cnt + 3 );
for(auto ei : tets_with_3_bnd_points)
for(auto j : Range(4))
{
auto i3_ = tempels[ei].GetFace (j);
ngcore::INT<3> i3 = {i3_[0], i3_[1], i3_[2]};
if(bnd_points[i3[0]] && bnd_points[i3[1]] && bnd_points[i3[2]])
{
i3.Sort();
face_table.Set( i3, true );
}
}
t1.Stop();
static Timer t2("check faces"); t2.Start();
openels.SetSize(0);
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
{
const Element2d & tri = mesh.OpenElement(i);
ngcore::INT<3> i3(tri[0], tri[1], tri[2]);
i3.Sort();
if(!face_table.Used(i3))
openels.Append(i);
}
t2.Stop();
auto p2sel = ngcore::CreateSortedTable<int, PointIndex>( Range(openels.Size()),
[&](auto & table, int i)
{
auto openel_i = openels[i];
const Element2d & tri = mesh.OpenElement(openel_i);
table.Add(tri[0], openel_i);
table.Add(tri[1], openel_i);
table.Add(tri[2], openel_i);
}, mesh.GetNP());
ngcore::BitArray badnode(mesh.GetNP()+PointIndex::BASE);
badnode.Clear();
ngcore::ParallelForRange(openels.Size(), [&] (auto myrange) {
for (auto i_ : myrange)
{
auto i = openels[i_];
const Element2d & tri = mesh.OpenElement(i);
for( auto edge : Range(3) )
{
auto pi0 = tri[edge];
auto pi1 = tri[(edge+1)%3];
if(pi0>pi1)
Swap(pi0, pi1);
// find other trig with edge pi0, pi1
int i_other = -1;
for(auto ii : p2sel[pi0])
{
if(ii==i)
continue;
auto & tri_other = mesh.OpenElement(ii);
if(tri_other[0]==pi1 || tri_other[1]==pi1 || tri_other[2]==pi1)
{
i_other = ii;
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--)
{
const DelaunayTet & el = tempels.Get(i);
if (badnode[el[0]] ||
badnode[el[1]] ||
badnode[el[2]] ||
badnode[el[3]] )
tempels.DeleteElement(i);
}
}
void DelaunayRemoveIntersecting( const Mesh & mesh, NgArray<DelaunayTet> & tempels, NgArray<int> & openels, Point3d pmin, Point3d pmax )
{
static Timer trem_intersect("Delaunay - remove intersecting"); RegionTimer rt(trem_intersect);
// find intersecting:
PrintMessage (3, "Remove intersecting");
if (openels.Size())
{
BoxTree<3> setree(pmin, pmax);
/*
cout << "open elements in search tree: " << openels.Size() << endl;
cout << "pmin, pmax = " << pmin << " - " << pmax << endl;
*/
for (int i = 1; i <= openels.Size(); i++)
{
int fnr;
fnr = openels.Get(i);
if (fnr)
{
const Element2d & tri = mesh.OpenElement(fnr);
Point3d ltpmin (mesh.Point(tri[0]));
Point3d ltpmax (ltpmin);
for (int k = 2; k <= 3; k++)
{
ltpmin.SetToMin (mesh.Point (tri.PNum(k)));
ltpmax.SetToMax (mesh.Point (tri.PNum(k)));
}
setree.Insert (ltpmin, ltpmax, fnr);
}
}
NgArray<int> neartrias;
for (int i = 1; i <= tempels.Size(); i++)
{
const Point<3> *pp[4];
int tetpi[4];
DelaunayTet & el = tempels.Elem(i);
int intersect = 0;
for (int j = 0; j < 4; j++)
{
pp[j] = &mesh.Point(el[j]);
tetpi[j] = el[j];
}
Point3d tetpmin(*pp[0]);
Point3d tetpmax(tetpmin);
for (int j = 1; j < 4; j++)
{
tetpmin.SetToMin (*pp[j]);
tetpmax.SetToMax (*pp[j]);
}
tetpmin = tetpmin + 0.01 * (tetpmin - tetpmax);
tetpmax = tetpmax + 0.01 * (tetpmax - tetpmin);
setree.GetIntersecting (tetpmin, tetpmax, neartrias);
// for (j = 1; j <= mesh.GetNSE(); j++)
// {
for (int jj = 1; jj <= neartrias.Size(); jj++)
{
int j = neartrias.Get(jj);
const Element2d & tri = mesh.OpenElement(j);
const Point<3> *tripp[3];
int tripi[3];
for (int k = 1; k <= 3; k++)
{
tripp[k-1] = &mesh.Point (tri.PNum(k));
tripi[k-1] = tri.PNum(k);
}
if (IntersectTetTriangle (&pp[0], &tripp[0], tetpi, tripi))
{
/*
int il1, il2;
(*testout) << "intersect !" << endl;
(*testout) << "triind: ";
for (il1 = 0; il1 < 3; il1++)
(*testout) << " " << tripi[il1];
(*testout) << endl;
(*testout) << "tetind: ";
for (il2 = 0; il2 < 4; il2++)
(*testout) << " " << tetpi[il2];
(*testout) << endl;
(*testout) << "trip: ";
for (il1 = 0; il1 < 3; il1++)
(*testout) << " " << *tripp[il1];
(*testout) << endl;
(*testout) << "tetp: ";
for (il2 = 0; il2 < 4; il2++)
(*testout) << " " << *pp[il2];
(*testout) << endl;
*/
intersect = 1;
break;
}
}
if (intersect)
{
tempels.DeleteElement(i);
i--;
}
}
}
}
void DelaunayRemoveOuter( const Mesh & mesh, NgArray<DelaunayTet> & tempels, AdFront3 * adfront )
{
static Timer trem_outer("Delaunay - remove outer"); RegionTimer rt(trem_outer);
PrintMessage (3, "Remove outer");
// find connected tets (with no face between, and no hole due
// to removed intersecting tets.
// INDEX_3_HASHTABLE<INDEX_2> innerfaces(np);
INDEX_3_HASHTABLE<int> boundaryfaces(mesh.GetNOpenElements()/3+1);
/*
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
{
const Element2d & tri = mesh.OpenElement(i);
INDEX_3 i3 (tri[0], tri[1], tri[2]);
i3.Sort();
boundaryfaces.PrepareSet (i3);
}
*/
for (const Element2d & tri : mesh.OpenElements())
{
INDEX_3 i3 (tri[0], tri[1], tri[2]);
i3.Sort();
boundaryfaces.PrepareSet (i3);
}
boundaryfaces.AllocateElements();
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
{
const Element2d & tri = mesh.OpenElement(i);
INDEX_3 i3 (tri[0], tri[1], tri[2]);
i3.Sort();
boundaryfaces.Set (i3, 1);
}
/*
for (int i = 0; i < tempels.Size(); i++)
for (int j = 0; j < 4; j++)
tempels[i].NB(j) = 0;
*/
for (auto & el : tempels)
for (int j = 0; j < 4; j++)
el.NB(j) = 0;
TABLE<int,PointIndex::BASE> elsonpoint(mesh.GetNP());
/*
for (int i = 0; i < tempels.Size(); i++)
{
const DelaunayTet & el = tempels[i];
*/
for (const DelaunayTet & el : tempels)
{
INDEX_4 i4(el[0], el[1], el[2], el[3]);
i4.Sort();
elsonpoint.IncSizePrepare (i4.I1());
elsonpoint.IncSizePrepare (i4.I2());
}
elsonpoint.AllocateElementsOneBlock();
for (int i = 0; i < tempels.Size(); i++)
{
const DelaunayTet & el = tempels[i];
INDEX_4 i4(el[0], el[1], el[2], el[3]);
i4.Sort();
elsonpoint.Add (i4.I1(), i+1);
elsonpoint.Add (i4.I2(), i+1);
}
// cout << "elsonpoint mem: ";
// elsonpoint.PrintMemInfo(cout);
INDEX_3_CLOSED_HASHTABLE<INDEX_2> faceht(100);
Element2d hel(TRIG);
// for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
for (PointIndex pi : mesh.Points().Range())
{
faceht.SetSize (4 * elsonpoint[pi].Size());
for (int ii = 0; ii < elsonpoint[pi].Size(); ii++)
{
int i = elsonpoint[pi][ii];
const DelaunayTet & el = tempels.Get(i);
for (int j = 1; j <= 4; j++)
{
el.GetFace (j-1, hel);
hel.Invert();
hel.NormalizeNumbering();
if (hel[0] == pi)
{
INDEX_3 i3(hel[0], hel[1], hel[2]);
if (!boundaryfaces.Used (i3))
{
if (faceht.Used (i3))
{
INDEX_2 i2 = faceht.Get(i3);
tempels.Elem(i).NB(j-1) = i2.I1();
tempels.Elem(i2.I1()).NB(i2.I2()-1) = i;
}
else
{
hel.Invert();
hel.NormalizeNumbering();
INDEX_3 i3i(hel[0], hel[1], hel[2]);
INDEX_2 i2(i, j);
faceht.Set (i3i, i2);
}
}
}
}
}
}
/*
for (i = 1; i <= tempels.Size(); i++)
{
const DelaunayTet & el = tempels.Get(i);
for (j = 1; j <= 4; j++)
{
INDEX_3 i3;
Element2d face;
el.GetFace1 (j, face);
for (int kk = 1; kk <= 3; kk++)
i3.I(kk) = face.PNum(kk);
i3.Sort();
if (!boundaryfaces.Used (i3))
{
if (innerfaces.Used(i3))
{
INDEX_2 i2;
i2 = innerfaces.Get(i3);
i2.I2() = i;
innerfaces.Set (i3, i2);
}
else
{
INDEX_2 i2;
i2.I1() = i;
i2.I2() = 0;
innerfaces.Set (i3, i2);
}
}
}
}
*/
/*
(*testout) << "nb elements:" << endl;
for (i = 1; i <= tempels.Size(); i++)
{
(*testout) << i << " ";
for (j = 1; j <= 4; j++)
(*testout) << tempels.Get(i).NB1(j) << " ";
(*testout) << endl;
}
(*testout) << "pairs:" << endl;
for (i = 1; i <= innerfaces.GetNBags(); i++)
for (j = 1; j <= innerfaces.GetBagSize(i); j++)
{
INDEX_3 i3;
INDEX_2 i2;
innerfaces.GetData (i, j, i3, i2);
(*testout) << i2 << endl;
}
*/
/*
cout << "innerfaces: ";
innerfaces.PrintMemInfo (cout);
*/
// cout << "boundaryfaces: ";
// boundaryfaces.PrintMemInfo (cout);
PrintMessage (5, "tables filled");
auto ne = tempels.Size();
NgBitArray inner(ne), outer(ne);
inner.Clear();
outer.Clear();
NgArray<int> elstack;
/*
int starti = 0;
for (i = 1; i <= ne; i++)
{
const Element & el = tempels.Get(i);
for (j = 1; j <= 4; j++)
for (k = 1; k <= 4; k++)
if (el.PNum(j) == startel.PNum(k))
{
outer.Set(i);
starti = i;
}
}
*/
int lowest_undefined_el = 1;
while (1)
{
int inside;
bool done = 1;
int i;
for (i = lowest_undefined_el; i <= ne; i++)
if (!inner.Test(i) && !outer.Test(i))
{
lowest_undefined_el = i+1;
done = 0;
break;
}
if (done) break;
const DelaunayTet & el = tempels.Get(i);
const Point3d & p1 = mesh.Point (el[0]);
const Point3d & p2 = mesh.Point (el[1]);
const Point3d & p3 = mesh.Point (el[2]);
const Point3d & p4 = mesh.Point (el[3]);
Point3d ci = Center (p1, p2, p3, p4);
inside = adfront->Inside (ci);
/*
cout << "startel: " << i << endl;
cout << "inside = " << inside << endl;
cout << "ins2 = " << adfront->Inside (Center (ci, p1)) << endl;
cout << "ins3 = " << adfront->Inside (Center (ci, p2)) << endl;
*/
elstack.SetSize(0);
elstack.Append (i);
while (elstack.Size())
{
int ei = elstack.Last();
elstack.DeleteLast();
if (!inner.Test(ei) && !outer.Test(ei))
{
if (inside)
inner.Set(ei);
else
outer.Set(ei);
for (int j = 1; j <= 4; j++)
{
INDEX_3 i3 = tempels.Get(ei).GetFace(j-1);
/*
Element2d face;
tempels.Get(ei).GetFace(j, face);
for (int kk = 1; kk <= 3; kk++)
i3.I(kk) = face.PNum(kk);
*/
i3.Sort();
if (tempels.Get(ei).NB(j-1))
elstack.Append (tempels.Get(ei).NB(j-1));
/*
if (innerfaces.Used(i3))
{
INDEX_2 i2 = innerfaces.Get(i3);
int other = i2.I1() + i2.I2() - ei;
if (other != tempels.Get(ei).NB1(j))
cerr << "different1 !!" << endl;
if (other)
{
elstack.Append (other);
}
}
else
if (tempels.Get(ei).NB1(j))
cerr << "different2 !!" << endl;
*/
}
}
}
}
// check outer elements
if (debugparam.slowchecks)
{
for (int i = 1; i <= ne; i++)
{
const DelaunayTet & el = tempels.Get(i);
const Point3d & p1 = mesh.Point (el[0]);
const Point3d & p2 = mesh.Point (el[1]);
const Point3d & p3 = mesh.Point (el[2]);
const Point3d & p4 = mesh.Point (el[3]);
Point3d ci = Center (p1, p2, p3, p4);
// if (adfront->Inside (ci) != adfront->Inside (Center (ci, p1)))
// cout << "ERROR: outer test unclear !!!" << endl;
if (inner.Test(i) != adfront->Inside (ci))
{
/*
cout << "ERROR: outer test wrong !!!"
<< "inner = " << int(inner.Test(i))
<< "outer = " << int(outer.Test(i))
<< endl;
cout << "Vol = " << Determinant(Vec3d(p1, p2),
Vec3d(p1, p3),
Vec3d(p1, p4)) << endl;
*/
for (int j = 1; j <= 4; j++)
{
Point3d hp;
switch (j)
{
case 1: hp = Center (ci, p1); break;
case 2: hp = Center (ci, p2); break;
case 3: hp = Center (ci, p3); break;
case 4: hp = Center (ci, p4); break;
}
// cout << "inside(" << hp << ") = " << adfront->Inside(hp) << endl;
}
}
if (adfront->Inside(ci))
outer.Clear(i);
else
outer.Set(i);
}
}
/*
// find bug in innerfaces
tempmesh.DeleteVolumeElements();
for (i = 1; i <= innerfaces.GetNBags(); i++)
for (j = 1; j <= innerfaces.GetBagSize(i); j++)
{
INDEX_3 i3;
INDEX_2 i2;
innerfaces.GetData (i, j, i3, i2);
if (i2.I2())
{
if (outer.Test(i2.I1()) != outer.Test(i2.I2()))
{
tempmesh.AddVolumeElement (tempels.Get(i2.I1()));
tempmesh.AddVolumeElement (tempels.Get(i2.I2()));
cerr << "outer flag different for connected els" << endl;
}
}
}
cout << "Check intersectiong once more" << endl;
for (i = 1; i <= openels.Size(); i++)
{
tempmesh.SurfaceElement(2*openels.Get(i)).SetIndex(2);
tempmesh.SurfaceElement(2*openels.Get(i)-1).SetIndex(2);
}
// for (i = 1; i <= tempmesh.GetNE(); i++)
// for (j = 1; j <= tempmesh.GetNSE(); j++)
i = 6; j = 403;
if (i <= tempmesh.GetNE() && j <= tempmesh.GetNSE())
if (tempmesh.SurfaceElement(j).GetIndex()==2)
{
const Element & el = tempmesh.VolumeElement(i);
const Element2d & sel = tempmesh.SurfaceElement(j);
const Point3d *tripp[3];
const Point3d *pp[4];
int tetpi[4], tripi[3];
for (k = 1; k <= 4; k++)
{
pp[k-1] = &tempmesh.Point(el.PNum(k));
tetpi[k-1] = el.PNum(k);
}
for (k = 1; k <= 3; k++)
{
tripp[k-1] = &tempmesh.Point (sel.PNum(k));
tripi[k-1] = sel.PNum(k);
}
(*testout) << "Check Triangle " << j << ":";
for (k = 1; k <= 3; k++)
(*testout) << " " << sel.PNum(k);
for (k = 1; k <= 3; k++)
(*testout) << " " << tempmesh.Point(sel.PNum(k));
(*testout) << endl;
(*testout) << "Check Tet " << i << ":";
for (k = 1; k <= 4; k++)
(*testout) << " " << el.PNum(k);
for (k = 1; k <= 4; k++)
(*testout) << " " << tempmesh.Point(el.PNum(k));
(*testout) << endl;
if (IntersectTetTriangle (&pp[0], &tripp[0], tetpi, tripi))
{
cout << "Intesection detected !!" << endl;
}
}
tempmesh.Save ("temp.vol");
// end bug search
*/
for (int i = ne; i >= 1; i--)
{
if (outer.Test(i))
tempels.DeleteElement(i);
}
// mesh.points.SetSize(mesh.points.Size()-4);
PrintMessage (5, "outer removed");
}
void Meshing3 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp)
{
static Timer t("Meshing3::Delaunay"); RegionTimer reg(t);
PrintMessage (1, "Delaunay meshing");
PrintMessage (3, "number of points: ", mesh.GetNP());
// PushStatus ("Delaunay meshing");
NgArray<DelaunayTet> tempels;
Point3d pmin, pmax;
DelaunayTet startel;
int oldnp = mesh.GetNP();
if (mp.blockfill)
{
BlockFillLocalH (mesh, mp);
PrintMessage (3, "number of points: ", mesh.GetNP());
}
int np = mesh.GetNP();
Delaunay1 (mesh, mp, adfront, tempels, oldnp, startel, pmin, pmax);
{
// improve delaunay - mesh by swapping !!!!
Mesh tempmesh;
tempmesh.GetMemoryTracer().SetName("delaunay-tempmesh");
for (auto & meshpoint : mesh.Points())
tempmesh.AddPoint (meshpoint);
for (auto & tempel : tempels)
{
Element el(4);
for (int j = 0; j < 4; j++)
el[j] = tempel[j];
el.SetIndex (1);
const Point<3> & lp1 = mesh.Point (el[0]);
const Point<3> & lp2 = mesh.Point (el[1]);
const Point<3> & lp3 = mesh.Point (el[2]);
const Point<3> & lp4 = mesh.Point (el[3]);
Vec<3> v1 = lp2-lp1;
Vec<3> v2 = lp3-lp1;
Vec<3> v3 = lp4-lp1;
Vec<3> n = Cross (v1, v2);
double vol = n * v3;
if (vol > 0) swap (el[2], el[3]);
tempmesh.AddVolumeElement (el);
}
tempels.DeleteAll();
MeshQuality3d (tempmesh);
tempmesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
tempmesh.AddFaceDescriptor (FaceDescriptor (2, 1, 0, 0));
for (int i = 1; i <= mesh.GetNOpenElements(); i++)
{
Element2d sel = mesh.OpenElement(i);
sel.SetIndex(1);
tempmesh.AddSurfaceElement (sel);
swap (sel[1], sel[2]);
tempmesh.AddSurfaceElement (sel);
}
for (int i = 1; i <= 4; i++)
{
Element2d self(TRIG);
self.SetIndex (1);
startel.GetFace (i-1, self);
tempmesh.AddSurfaceElement (self);
}
// for (i = mesh.GetNP() - 3; i <= mesh.GetNP(); i++)
// tempmesh.AddLockedPoint (i);
for (auto pi : tempmesh.Points().Range())
tempmesh.AddLockedPoint (pi);
// tempmesh.PrintMemInfo(cout);
// tempmesh.Save ("tempmesh.vol");
{
MeshOptimize3d meshopt(mp);
tempmesh.Compress();
tempmesh.FindOpenElements ();
#ifndef EMSCRIPTEN
RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
#endif // EMSCRIPTEN
for (auto i : Range(10))
{
PrintMessage (5, "Num open: ", tempmesh.GetNOpenElements());
if(i%5==0)
tempmesh.FreeOpenElementsEnvironment (1);
meshopt.SwapImprove(tempmesh, OPT_CONFORM);
}
tempmesh.Compress();
}
MeshQuality3d (tempmesh);
tempels.SetSize(tempmesh.GetNE());
tempels.SetSize(0);
for (auto & el : tempmesh.VolumeElements())
tempels.Append (el);
}
DelaunayRemoveDegenerated(mesh.Points(), tempels, np);
NgArray<int> openels;
DelaunayRemoveTwoTriaTets(mesh, tempels, openels);
DelaunayRemoveIntersecting(mesh, tempels, openels, pmin, pmax);
DelaunayRemoveOuter(mesh, tempels, adfront);
for (int i = 0; i < tempels.Size(); i++)
{
Element el(4);
for (int j = 0; j < 4; j++)
el[j] = tempels[i][j];
mesh.AddVolumeElement (el);
}
mesh.FindOpenElements(domainnr);
mesh.Compress();
// PopStatus ();
}
}