netgen/libsrc/meshing/delaunay2d.cpp
2017-11-10 13:22:20 +01:00

423 lines
11 KiB
C++

#include <mystdlib.h>
#include "meshing.hpp"
// not yet working ....
namespace netgen
{
class DelaunayTrig
{
PointIndex pnums[3];
Point<3> c;
double r;
double rad2;
public:
DelaunayTrig () { ; }
DelaunayTrig (int p1, int p2, int p3)
{ pnums[0] = p1; pnums[1] = p2; pnums[2] = p3; }
PointIndex & operator[] (int j) { return pnums[j]; }
const PointIndex & operator[] (int j) const { return pnums[j]; }
void CalcCenter (Mesh & mesh)
{
Point<3> p1 = mesh[pnums[0]];
Point<3> p2 = mesh[pnums[1]];
Point<3> p3 = mesh[pnums[2]];
Vec<3> v1 = p2-p1;
Vec<3> v2 = p3-p1;
Mat<2,2> mat, inv;
mat(0,0) = v1*v1;
mat(0,1) = v1*v2;
mat(1,0) = v2*v1;
mat(1,1) = v2*v2;
Vec<2> rhs, sol;
rhs(0) = 0.5 * v1*v1;
rhs(1) = 0.5 * v2*v2;
CalcInverse (mat, inv);
sol = inv * rhs;
c = p1 + sol(0) * v1 + sol(1) * v2;
rad2 = Dist2(c, p1);
r = sqrt(rad2);
}
Point<3> Center() const { return c; }
double Radius2() const { return rad2; }
Box<3> BoundingBox() const { return Box<3> (c-Vec<3>(r,r,0.1), c+Vec<3>(r,r,0.1)); }
};
ostream & operator<< (ostream & ost, DelaunayTrig trig)
{
ost << trig[0] << "-" << trig[1] << "-" << trig[2] << endl;
return ost;
}
void Meshing2 :: BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp)
{
static int timer = NgProfiler::CreateTimer ("Meshing2::BlockFill");
static int timer1 = NgProfiler::CreateTimer ("Meshing2::BlockFill 1");
static int timer2 = NgProfiler::CreateTimer ("Meshing2::BlockFill 2");
static int timer3 = NgProfiler::CreateTimer ("Meshing2::BlockFill 3");
static int timer4 = NgProfiler::CreateTimer ("Meshing2::BlockFill 4");
NgProfiler::RegionTimer reg (timer);
NgProfiler::StartTimer (timer1);
double filldist = mp.filldist;
cout << "blockfill local h" << endl;
cout << "rel filldist = " << filldist << endl;
PrintMessage (3, "blockfill local h");
Array<Point<3> > npoints;
// adfront -> CreateTrees();
Box<3> bbox ( Box<3>::EMPTY_BOX );
double maxh = 0;
for (int i = 0; i < adfront->GetNFL(); i++)
{
const FrontLine & line = adfront->GetLine (i);
const Point<3> & p1 = adfront->GetPoint(line.L().I1());
const Point<3> & p2 = adfront->GetPoint(line.L().I2());
maxh = max (maxh, Dist (p1, p2));
bbox.Add (p1);
bbox.Add (p2);
}
cout << "bbox = " << bbox << endl;
// Point<3> mpc = bbox.Center();
bbox.Increase (bbox.Diam()/2);
Box<3> meshbox = bbox;
NgProfiler::StopTimer (timer1);
NgProfiler::StartTimer (timer2);
LocalH loch2 (bbox, 1, 2);
if (mp.maxh < maxh) maxh = mp.maxh;
bool changed;
do
{
mesh.LocalHFunction().ClearFlags();
for (int i = 0; i < adfront->GetNFL(); i++)
{
const FrontLine & line = adfront->GetLine(i);
Box<3> bbox (adfront->GetPoint (line.L().I1()));
bbox.Add (adfront->GetPoint (line.L().I2()));
double filld = filldist * bbox.Diam();
bbox.Increase (filld);
mesh.LocalHFunction().CutBoundary (bbox);
}
mesh.LocalHFunction().FindInnerBoxes (adfront, NULL);
npoints.SetSize(0);
mesh.LocalHFunction().GetInnerPoints (npoints);
changed = false;
for (int i = 0; i < npoints.Size(); i++)
{
if (mesh.LocalHFunction().GetH(npoints[i]) > 1.2 * maxh)
{
mesh.LocalHFunction().SetH (npoints[i], maxh);
changed = true;
}
}
}
while (changed);
NgProfiler::StopTimer (timer2);
NgProfiler::StartTimer (timer3);
if (debugparam.slowchecks)
(*testout) << "Blockfill with points: " << endl;
*testout << "loch = " << mesh.LocalHFunction() << endl;
*testout << "npoints = " << endl << npoints << endl;
for (int i = 1; i <= npoints.Size(); i++)
{
if (meshbox.IsIn (npoints.Get(i)))
{
PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
adfront->AddPoint (npoints.Get(i), gpnum);
if (debugparam.slowchecks)
{
(*testout) << npoints.Get(i) << endl;
Point<2> p2d (npoints.Get(i)(0), npoints.Get(i)(1));
if (!adfront->Inside(p2d))
{
cout << "add outside point" << endl;
(*testout) << "outside" << endl;
}
}
}
}
NgProfiler::StopTimer (timer3);
NgProfiler::StartTimer (timer4);
// find outer points
loch2.ClearFlags();
for (int i = 0; i < adfront->GetNFL(); i++)
{
const FrontLine & line = adfront->GetLine(i);
Box<3> bbox (adfront->GetPoint (line.L().I1()));
bbox.Add (adfront->GetPoint (line.L().I2()));
loch2.SetH (bbox.Center(), bbox.Diam());
}
for (int i = 0; i < adfront->GetNFL(); i++)
{
const FrontLine & line = adfront->GetLine(i);
Box<3> bbox (adfront->GetPoint (line.L().I1()));
bbox.Add (adfront->GetPoint (line.L().I2()));
bbox.Increase (filldist * bbox.Diam());
loch2.CutBoundary (bbox);
}
loch2.FindInnerBoxes (adfront, NULL);
// outer points : smooth mesh-grading
npoints.SetSize(0);
loch2.GetOuterPoints (npoints);
for (int i = 1; i <= npoints.Size(); i++)
{
if (meshbox.IsIn (npoints.Get(i)))
{
PointIndex gpnum = mesh.AddPoint (npoints.Get(i));
adfront->AddPoint (npoints.Get(i), gpnum);
}
}
NgProfiler::StopTimer (timer4);
}
void Meshing2 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp)
{
static int timer = NgProfiler::CreateTimer ("Meshing2::Delaunay - total");
static int timerstart = NgProfiler::CreateTimer ("Meshing2::Delaunay - start");
static int timerfinish = NgProfiler::CreateTimer ("Meshing2::Delaunay - finish");
static int timer1 = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental");
static int timer1a = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental a");
static int timer1b = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental b");
static int timer1c = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental c");
static int timer1d = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental d");
NgProfiler::RegionTimer reg (timer);
cout << "2D Delaunay meshing (in progress)" << endl;
BlockFillLocalH (mesh, mp);
NgProfiler::StartTimer (timerstart);
// do the delaunay
// face bounding box:
Box<3> bbox (Box<3>::EMPTY_BOX);
for (int i = 0; i < adfront->GetNFL(); i++)
{
const FrontLine & line = adfront->GetLine(i);
bbox.Add (Point<3> (adfront->GetPoint (line.L()[0])));
bbox.Add (Point<3> (adfront->GetPoint (line.L()[1])));
}
for (int i = 0; i < mesh.LockedPoints().Size(); i++)
bbox.Add (mesh.Point (mesh.LockedPoints()[i]));
cout << "bbox = " << bbox << endl;
// external point
Vec<3> vdiag = bbox.PMax()-bbox.PMin();
auto old_points = mesh.Points().Range();
DelaunayTrig startel;
startel[0] = mesh.AddPoint (bbox.PMin() + Vec<3> (-8*vdiag(0), -8*vdiag(1), 0));
startel[1] = mesh.AddPoint (bbox.PMin() + Vec<3> (+8*vdiag(0), -8*vdiag(1), 0));
startel[2] = mesh.AddPoint (bbox.PMin() + Vec<3> (0, 8*vdiag(1), 0));
Box<3> hbox;
hbox.Set (mesh[startel[0]]);
hbox.Add (mesh[startel[1]]);
hbox.Add (mesh[startel[2]]);
Point<3> hp = mesh[startel[0]];
hp(2) = 1; hbox.Add (hp);
hp(2) = -1; hbox.Add (hp);
BoxTree<3> searchtree(hbox);
Array<DelaunayTrig> tempels;
startel.CalcCenter (mesh);
tempels.Append (startel);
searchtree.Insert(startel.BoundingBox(), 0);
Array<int> closeels;
Array<int> intersecting;
Array<INDEX_2> edges;
// reorder points
Array<PointIndex, PointIndex::BASE, PointIndex> mixed(old_points.Size());
int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
int prim;
{
int i = 0;
while (old_points.Size() % prims[i] == 0) i++;
prim = prims[i];
}
for (PointIndex pi : old_points)
mixed[pi] = PointIndex ( (prim * pi) % old_points.Size() + PointIndex::BASE );
NgProfiler::StopTimer (timerstart);
NgProfiler::StartTimer (timer1);
for (PointIndex i1 : old_points)
{
PointIndex i = mixed[i1];
NgProfiler::StartTimer (timer1a);
Point<3> newp = mesh[i];
intersecting.SetSize(0);
edges.SetSize(0);
searchtree.GetIntersecting (newp, newp, closeels);
// for (int jj = 0; jj < closeels.Size(); jj++)
// for (int j = 0; j < tempels.Size(); j++)
for (int j : closeels)
{
if (tempels[j][0] < 0) continue;
Point<3> c = tempels[j].Center();
double r2 = tempels[j].Radius2();
bool inside = Dist2(mesh[i], c) < r2;
if (inside) intersecting.Append (j);
}
NgProfiler::StopTimer (timer1a);
NgProfiler::StartTimer (timer1b);
// find outer edges
for (auto j : intersecting)
{
const DelaunayTrig & trig = tempels[j];
for (int k = 0; k < 3; k++)
{
int p1 = trig[k];
int p2 = trig[(k+1)%3];
INDEX_2 edge(p1,p2);
edge.Sort();
bool found = false;
for (int l = 0; l < edges.Size(); l++)
if (edges[l] == edge)
{
edges.Delete(l);
found = true;
break;
}
if (!found) edges.Append (edge);
}
}
NgProfiler::StopTimer (timer1b);
NgProfiler::StartTimer (timer1c);
/*
for (int j = intersecting.Size()-1; j >= 0; j--)
tempels.Delete (intersecting[j]);
*/
for (int j : intersecting)
{
searchtree.DeleteElement (j);
tempels[j][0] = -1;
tempels[j][1] = -1;
tempels[j][2] = -1;
}
NgProfiler::StopTimer (timer1c);
NgProfiler::StartTimer (timer1d);
for (auto edge : edges)
{
DelaunayTrig trig (edge[0], edge[1], i);
trig.CalcCenter (mesh);
tempels.Append (trig);
searchtree.Insert(trig.BoundingBox(), tempels.Size()-1);
}
NgProfiler::StopTimer (timer1d);
}
NgProfiler::StopTimer (timer1);
NgProfiler::StartTimer (timerfinish);
for (DelaunayTrig & trig : tempels)
{
if (trig[0] < 0) continue;
Point<3> c = Center (mesh[trig[0]], mesh[trig[1]], mesh[trig[2]]);
if (!adfront->Inside (Point<2> (c(0),c(1)))) continue;
Vec<3> n = Cross (mesh[trig[1]]-mesh[trig[0]],
mesh[trig[2]]-mesh[trig[0]]);
if (n(2) < 0) Swap (trig[1], trig[2]);
Element2d el(trig[0], trig[1], trig[2]);
el.SetIndex (domainnr);
mesh.AddSurfaceElement (el);
}
for (PointIndex pi : mesh.Points().Range())
*testout << pi << ": " << mesh[pi].Type() << endl;
NgProfiler::StopTimer (timerfinish);
}
}