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

Delaunay for 2d mesh generation

Browse Source 
Squashed commit of the following:

commit 84f36ffeb409f5fddb389c75ee48b4872b516ae9
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 18:27:15 2020 +0200

    revert change in spline partitioning

commit d4aef23a22a9beb26c4453267c99dd7533174ced
Merge: 15a467aa 97dfecd0
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 17:59:00 2020 +0200

    Merge branch 'master' into delaunay2d

commit 15a467aa7f7cb09f9ea3d984905fe3da69f0b238
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 17:44:31 2020 +0200

    delaunay2d - fix trig orientation

commit be223412ad972722a51b64a5bccf7ca2bec566c8
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 17:29:46 2020 +0200

    fix delaunay swapping

commit 48b95ae2ee1cbabcfae79dfd1cb7af1fd69d77f3
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 17:27:51 2020 +0200

    testout only with debug settings

commit d82b7a7cecb6f65f42b79b666fc58d0116dc0365
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 16:37:10 2020 +0200

    delaunay only for large domains

commit 1f51eaca1ff7a3777e4f32ba9e35e48d496a2854
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 16:21:33 2020 +0200

    compress points in delaunay

commit 20a223f36f3912a208db80c717d9dd87851ba43f
Merge: 2446b746 4c15146d
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 14:31:14 2020 +0200

    Merge branch 'master' into delaunay2d

commit 2446b74687ee56633a86e748e85343919edbd5ad
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 14:22:01 2020 +0200

    optimize CalcPartition() and PartitionBoundary()

commit 3baa58833348a72f16853530a5d17e73424186df
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 16 12:24:17 2020 +0200

    MeshingParameters - delaunay2d option (default is off)

commit e79b113dde9b9c4c5b92239817c6058ca468c319
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Thu Oct 15 16:12:45 2020 +0200

    fix windows build error

commit 92c7b9c1ed4016458980bbc21c61dae07f4444c7
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Thu Oct 15 17:51:44 2020 +0200

    delaunay bugfix

commit 6880194107819cfb2d23206e7e0f48ff5aa3fc10
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Thu Oct 15 17:51:05 2020 +0200

    csg2d - fix bug with splines

commit 1d9baa299d49e1f6fa16f4368885601ed01c5de7
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 18:40:23 2020 +0200

    CSG2d - faster AddIntersections (search tree per loop)

commit 2679ef0dee10cdf486441af47ca4c081aa7eb50b
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 17:13:17 2020 +0200

    bounding box for Loop

commit 894c6345b737693e32cbda368e56f56764e11085
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 12:48:51 2020 +0200

    remove debug check, output in blockfill

commit 2b0a0892c41e746b12e5e852cdb138acd3d2c4e3
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 12:40:53 2020 +0200

    compress mesh after delaunay

commit 1de33c87eee3199d4d9b18544f66e53329b47a2f
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 12:37:07 2020 +0200

    revert change in improve2d

commit 41a60e89533e94b93b92202ac17852d3aee9acbb
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 12:25:07 2020 +0200

    cleanup delaunay2d

commit c16aae324969cd5a90748953019933690d013337
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 11:39:56 2020 +0200

    sunburst chart - tooltip formatting

commit 4d61e1fdeab302ba357904f22f951361935791f0
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Oct 14 11:03:37 2020 +0200

    delaunay seems to work

commit 8bd43f54d1efd6862f1b403cdb6c8ce9b5f7b3c6
Merge: 90ac7adb 25efdadd
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Tue Oct 13 12:08:01 2020 +0200

    Merge remote-tracking branch 'gitlab/master' into delaunay2d

commit 90ac7adb562cf2402345c5dfb4281bd097b5d62d
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Tue Oct 13 12:04:49 2020 +0200

    fix Loop::operator=

commit 1eb4f2de3b6576f503a073011a208fa8f609524e
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Tue Oct 13 12:04:13 2020 +0200

    more statistics in sunburst chart

commit db8b97ffbbc7db2a3413c4f8a5528eebe3488d57
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Tue Oct 13 11:17:28 2020 +0200

    more work on delaunay2d

commit eaa675f2351252b5fde423f241b10e231d1eb97e
Merge: 0eb9f9bd 8f837cb9
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Mon Oct 12 12:51:31 2020 +0200

    Merge remote-tracking branch 'gitlab/delaunay2d' into delaunay2d

commit 0eb9f9bd1c31a0e3c3c796c9280b1c1d007ace26
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Mon Oct 12 12:50:10 2020 +0200

    further csg2d optimization

commit 8f837cb9a281acca7c33159985da3b6992fe638f
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Tue Oct 6 19:02:31 2020 +0200

    csg2d - optimize Loop::operator=

commit 7bb4f16b886b20902b0d3563716055fc1734d47e
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Tue Oct 6 10:28:20 2020 +0200

    csg2d performance (search tree, inside tests)

commit 2c9ebce04d7989223327a1875e1b65bf180c95f5
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Fri Oct 2 16:33:24 2020 +0200

    [WIP] delaunay2d

commit 749df2311a3ac1976faaa9f0b60846709a2087b9
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Thu Oct 1 11:36:03 2020 +0200

    something

commit cda9fffde33a86b71467debb86848fdb9cfbf80c
Author: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date:   Wed Sep 30 12:06:53 2020 +0200

    delaunay2d - fix size of starting trig
This commit is contained in:
Matthias Hochsteger 2020-10-16 18:28:09 +02:00
parent 97dfecd040
commit 8ba6bad6fd
4 changed files with 656 additions and 204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
libsrc/geom2d/genmesh2d.cpp
View File

@ -392,10 +392,16 @@ namespace netgen
shared_ptr<Mesh> & mesh, shared_ptr<Mesh> & mesh,
MeshingParameters & mp) MeshingParameters & mp)
{ {
static Timer tall("MeshFromSpline2D"); RegionTimer rtall(tall);
static Timer t_h("SetH");
static Timer t_tensor("tensor domain meshing");
static Timer t_part_boundary("PartitionBoundary");
static Timer t_hpref("mark hpref points");
PrintMessage (1, "Generate Mesh from spline geometry"); PrintMessage (1, "Generate Mesh from spline geometry");
Box<2> bbox = geometry.GetBoundingBox (); Box<2> bbox = geometry.GetBoundingBox ();
t_h.Start();
if (bbox.Diam() < mp.maxh) if (bbox.Diam() < mp.maxh)
mp.maxh = bbox.Diam(); mp.maxh = bbox.Diam();
@ -412,11 +418,14 @@ namespace netgen
t_part_boundary.Start();
geometry.PartitionBoundary (mp, mp.maxh, *mesh); geometry.PartitionBoundary (mp, mp.maxh, *mesh);
t_part_boundary.Stop();
PrintMessage (3, "Boundary mesh done, np = ", mesh->GetNP()); PrintMessage (3, "Boundary mesh done, np = ", mesh->GetNP());
t_hpref.Start();
// marks mesh points for hp-refinement // marks mesh points for hp-refinement
for (int i = 0; i < geometry.GetNP(); i++) for (int i = 0; i < geometry.GetNP(); i++)
if (geometry.GetPoint(i).hpref) if (geometry.GetPoint(i).hpref)
@ -434,6 +443,7 @@ namespace netgen
} }
(*mesh)[mpi].Singularity(geometry.GetPoint(i).hpref); (*mesh)[mpi].Singularity(geometry.GetPoint(i).hpref);
} }
t_hpref.Stop();
int maxdomnr = 0; int maxdomnr = 0;
@ -459,6 +469,7 @@ namespace netgen
mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) ); mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) );
mesh->CalcLocalH(mp.grading); mesh->CalcLocalH(mp.grading);
t_h.Stop();
int bnp = mesh->GetNP(); // boundary points int bnp = mesh->GetNP(); // boundary points
auto BndPntRange = mesh->Points().Range(); auto BndPntRange = mesh->Points().Range();
@ -469,6 +480,7 @@ namespace netgen
for (int domnr = 1; domnr <= maxdomnr; domnr++) for (int domnr = 1; domnr <= maxdomnr; domnr++)
if (geometry.GetDomainTensorMeshing (domnr)) if (geometry.GetDomainTensorMeshing (domnr))
{ // tensor product mesh { // tensor product mesh
RegionTimer rt(t_tensor);
NgArray<PointIndex, PointIndex::BASE> nextpi(bnp); NgArray<PointIndex, PointIndex::BASE> nextpi(bnp);
NgArray<int, PointIndex::BASE> si1(bnp), si2(bnp); NgArray<int, PointIndex::BASE> si1(bnp), si2(bnp);
@ -556,8 +568,10 @@ namespace netgen
static Timer t_domain("Mesh domain");
for (int domnr = 1; domnr <= maxdomnr; domnr++) for (int domnr = 1; domnr <= maxdomnr; domnr++)
{ {
RegionTimer rt(t_domain);
if (geometry.GetDomainTensorMeshing (domnr)) continue; if (geometry.GetDomainTensorMeshing (domnr)) continue;
double h = mp.maxh; double h = mp.maxh;
@ -573,16 +587,26 @@ namespace netgen
Meshing2 meshing (geometry, mp, Box<3> (pmin, pmax)); Meshing2 meshing (geometry, mp, Box<3> (pmin, pmax));
NgArray<int, PointIndex::BASE> compress(bnp); NgArray<int, PointIndex::BASE> compress(mesh->GetNP());
compress = -1; compress = -1;
int cnt = 0; int cnt = 0;
for (PointIndex pi : BndPntRange)
if ( (*mesh)[pi].GetLayer() == geometry.GetDomainLayer(domnr)) for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
{ {
meshing.AddPoint ((*mesh)[pi], pi); const auto & s = (*mesh)[si];
cnt++; if ( s.domin==domnr || s.domout==domnr )
compress[pi] = cnt; {
} for (auto pi : {s[0], s[1]})
{
if(compress[pi]==-1)
{
meshing.AddPoint((*mesh)[pi], pi);
cnt++;
compress[pi] = cnt;
}
}
}
}
PointGeomInfo gi; PointGeomInfo gi;
gi.trignum = 1; gi.trignum = 1;
@ -600,12 +624,15 @@ namespace netgen
} }
} }
// not complete, use at own risk ... if(mp.delaunay2d && cnt>100)
// meshing.Delaunay(*mesh, domnr, mp); meshing.Delaunay(*mesh, domnr, mp);
mp.checkoverlap = 0; else
auto res = meshing.GenerateMesh (*mesh, mp, h, domnr); {
if (res != 0) // mp.checkoverlap = 0;
throw NgException("meshing failed"); auto res = meshing.GenerateMesh (*mesh, mp, h, domnr);
if (res != 0)
throw NgException("meshing failed");
}
for (SurfaceElementIndex sei = oldnf; sei < mesh->GetNSE(); sei++) for (SurfaceElementIndex sei = oldnf; sei < mesh->GetNSE(); sei++)
(*mesh)[sei].SetIndex (domnr); (*mesh)[sei].SetIndex (domnr);

804
libsrc/meshing/delaunay2d.cpp
View File

@ -1,34 +1,52 @@
#include <mystdlib.h> #include <mystdlib.h>
#include "meshing.hpp" #include "meshing.hpp"
#include <geom2d/csg2d.hpp>
// not yet working .... // not yet working ....
namespace netgen namespace netgen
{ {
using ngcore::INT;
extern void Optimize2d (Mesh & mesh, MeshingParameters & mp);
static inline Point<2> P2( Point<3> p )
{
return {p[0], p[1]};
}
static inline Point<3> P3( Point<2> p )
{
return {p[0], p[1], 0};
}
class DelaunayTrig class DelaunayTrig
{ {
PointIndex pnums[3]; PointIndex pnums[3];
Point<3> c; Point<2> c;
public:
double r; double r;
double rad2; double rad2;
public: DelaunayTrig () = default;
DelaunayTrig () { ; }
DelaunayTrig (int p1, int p2, int p3) DelaunayTrig (int p1, int p2, int p3)
{ pnums[0] = p1; pnums[1] = p2; pnums[2] = p3; } {
pnums[0] = p1;
pnums[1] = p2;
pnums[2] = p3;
}
PointIndex & operator[] (int j) { return pnums[j]; } PointIndex & operator[] (int j) { return pnums[j]; }
const PointIndex & operator[] (int j) const { return pnums[j]; } const PointIndex & operator[] (int j) const { return pnums[j]; }
void CalcCenter (Mesh & mesh) void CalcCenter (Mesh & mesh)
{ {
Point<3> p1 = mesh[pnums[0]]; Point<2> p1 = P2(mesh[pnums[0]]);
Point<3> p2 = mesh[pnums[1]]; Point<2> p2 = P2(mesh[pnums[1]]);
Point<3> p3 = mesh[pnums[2]]; Point<2> p3 = P2(mesh[pnums[2]]);
Vec<3> v1 = p2-p1; Vec<2> v1 = p2-p1;
Vec<3> v2 = p3-p1; Vec<2> v2 = p3-p1;
Mat<2,2> mat, inv; Mat<2,2> mat, inv;
mat(0,0) = v1*v1; mat(0,0) = v1*v1;
mat(0,1) = v1*v2; mat(0,1) = v1*v2;
@ -45,9 +63,301 @@ namespace netgen
r = sqrt(rad2); r = sqrt(rad2);
} }
Point<3> Center() const { return c; } Point<2> Center() const { return c; }
double Radius2() const { return rad2; } 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)); } Box<2> BoundingBox() const { return Box<2> (c-Vec<2>(r,r), c+Vec<2>(r,r)); }
};
class DelaunayMesh
{
ngcore::ClosedHashTable<INT<2>, INT<2>> edge_to_trig;
Array<DelaunayTrig> trigs;
unique_ptr<DelaunayTree<2>> tree;
Mesh & mesh;
Array<int> closeels;
Array<int> intersecting;
Array<INT<2>> edges;
int GetNeighbour( int eli, int edge )
{
auto p0 = trigs[eli][(edge+1)%3];
auto p1 = trigs[eli][(edge+2)%3];
if(p1<p0)
Swap(p0,p1);
INT<2> hash = {p0,p1};
if(!edge_to_trig.Used(hash))
return -1;
auto i2 = edge_to_trig.Get({p0,p1});
return i2[0] == eli ? i2[1] : i2[0];
}
void SetNeighbour( int eli, int edge )
{
auto p0 = trigs[eli][(edge+1)%3];
auto p1 = trigs[eli][(edge+2)%3];
if(p1<p0)
Swap(p0,p1);
INT<2> hash = {p0,p1};
if(!edge_to_trig.Used(hash))
edge_to_trig[hash] = {eli, -1};
else
{
auto i2 = edge_to_trig.Get({p0,p1});
if(i2[0]==-1)
i2[0] = eli;
else
{
if(i2[1]==-1)
i2[1] = eli;
}
edge_to_trig[hash] = i2;
}
}
void UnsetNeighbours( int eli )
{
for(int edge : Range(3))
{
auto p0 = trigs[eli][(edge+1)%3];
auto p1 = trigs[eli][(edge+2)%3];
if(p1<p0)
Swap(p0,p1);
INT<2> hash = {p0,p1};
auto i2 = edge_to_trig.Get({p0,p1});
if(i2[0]==eli)
i2[0] = i2[1];
i2[1] = -1;
edge_to_trig[hash] = i2;
}
}
void AppendTrig( int pi0, int pi1, int pi2 )
{
DelaunayTrig el;
el[0] = pi0;
el[1] = pi1;
el[2] = pi2;
el.CalcCenter(mesh);
trigs.Append(el);
int ti = trigs.Size()-1;
tree->Insert(el.BoundingBox(), ti);
for(int i : Range(3))
SetNeighbour(ti, i);
}
public:
DelaunayMesh( Mesh & mesh_, Box<2> box )
: mesh(mesh_)
{
Vec<2> vdiag = box.PMax()-box.PMin();
double w = vdiag(0);
double h = vdiag(1);
Point<2> p0 = box.PMin() + Vec<2> ( -3*h, -h);
Point<2> p1 = box.PMin() + Vec<2> (w+3*h, -h);
Point<2> p2 = box.Center() + Vec<2> (0, 1.5*h+0.5*w);
box.Add( p0 );
box.Add( p1 );
box.Add( p2 );
tree = make_unique<DelaunayTree<2>>(box);
auto pi0 = mesh.AddPoint (P3(p0));
auto pi1 = mesh.AddPoint (P3(p1));
auto pi2 = mesh.AddPoint (P3(p2));
AppendTrig(pi0, pi1, pi2);
}
void AddPoint( PointIndex pi_new )
{
Point<2> newp = P2(mesh[pi_new]);
intersecting.SetSize(0);
edges.SetSize(0);
int definitive_overlapping_trig = -1;
double minquot{1e20};
tree->GetFirstIntersecting (newp, newp, [&] (const auto i_trig)
{
const auto trig = trigs[i_trig];
double rad2 = trig.Radius2();
double d2 = Dist2 (trig.Center(), newp);
if (d2 >= rad2) return false;
if (d2 < 0.999 * rad2)
{
definitive_overlapping_trig = i_trig;
return true;
}
if (definitive_overlapping_trig == -1 || d2 < 0.99*minquot*rad2)
{
minquot = d2/rad2;
definitive_overlapping_trig = i_trig;
}
return false;
});
if(definitive_overlapping_trig==-1)
{
cout << "Error in tree - didnt find overlap, check all trigs again" << endl;
for(auto i_trig : trigs.Range())
{
const auto trig = trigs[i_trig];
if(trig[0]==-1)
continue;
double rad2 = trig.Radius2();
double d2 = Dist2 (trig.Center(), newp);
if (d2 < 0.999 * rad2)
{
definitive_overlapping_trig = i_trig;
break;
}
}
}
BitArray trig_visited(trigs.Size());
trig_visited.Clear();
if(definitive_overlapping_trig==-1)
throw Exception("point not in any circle "+ ToString(pi_new));
Array<int> trigs_to_visit;
trigs_to_visit.Append(definitive_overlapping_trig);
intersecting.Append(definitive_overlapping_trig);
trig_visited.SetBit(definitive_overlapping_trig);
while(trigs_to_visit.Size())
{
int ti = trigs_to_visit.Last();
trigs_to_visit.DeleteLast();
trig_visited.SetBit(ti);
auto & trig = trigs[ti];
for(auto ei : Range(3))
{
auto nb = GetNeighbour(ti, ei);
if(nb==-1)
continue;
if(trig_visited.Test(nb))
continue;
const auto & trig_nb = trigs[nb];
trig_visited.SetBit(nb);
bool is_intersecting = Dist2(newp, trig_nb.Center()) < trig_nb.Radius2()*(1+1e-12);
if(!is_intersecting)
{
const Point<2> p0 = P2(mesh[PointIndex (trig[(ei+1)%3])]);
const Point<2> p1 = P2(mesh[PointIndex (trig[(ei+2)%3])]);
const Point<2> p2 = P2(mesh[PointIndex (trig[ei])]);
auto v = p1-p0;
Vec<2> n = {-v[1], v[0]};
n /= n.Length();
double dist = n * (newp-p1);
double scal = n * (p2 - p1);
if (scal > 0) dist *= -1;
if (dist > -1e-10)
is_intersecting = true;
}
if(is_intersecting)
{
trigs_to_visit.Append(nb);
intersecting.Append(nb);
}
}
}
// find outer edges
for (auto j : intersecting)
{
const DelaunayTrig & trig = trigs[j];
for (int k = 0; k < 3; k++)
{
int p1 = trig[k];
int p2 = trig[(k+1)%3];
INT<2> edge{p1,p2};
edge.Sort();
bool found = false;
for (int l = 0; l < edges.Size(); l++)
if (edges[l] == edge)
{
edges.RemoveElement(l);
found = true;
break;
}
if (!found)
edges.Append (edge);
}
}
for (int j : intersecting)
{
UnsetNeighbours(j);
trigs[j][0] = -1;
trigs[j][1] = -1;
trigs[j][2] = -1;
}
for (auto edge : edges)
AppendTrig( edge[0], edge[1], pi_new );
for (int j : intersecting)
tree->DeleteElement (j);
static int counter=0;
if(0)
{
Mesh m;
m = mesh;
m.ClearSurfaceElements();
for (DelaunayTrig & trig : trigs)
{
if (trig[0] < 0) 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 (1);
m.AddSurfaceElement (el);
}
m.Save("meshes/mesh_"+ToString(counter++)+".vol.gz");
}
}
Array<DelaunayTrig> & GetElements() { return trigs; }
}; };
ostream & operator<< (ostream & ost, DelaunayTrig trig) ostream & operator<< (ostream & ost, DelaunayTrig trig)
@ -59,20 +369,18 @@ namespace netgen
void Meshing2 :: BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp) void Meshing2 :: BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp)
{ {
static int timer = NgProfiler::CreateTimer ("Meshing2::BlockFill"); static Timer timer("Meshing2::BlockFill");
static int timer1 = NgProfiler::CreateTimer ("Meshing2::BlockFill 1"); static Timer timer1("Meshing2::BlockFill 1");
static int timer2 = NgProfiler::CreateTimer ("Meshing2::BlockFill 2"); static Timer timer2("Meshing2::BlockFill 2");
static int timer3 = NgProfiler::CreateTimer ("Meshing2::BlockFill 3"); static Timer timer3("Meshing2::BlockFill 3");
static int timer4 = NgProfiler::CreateTimer ("Meshing2::BlockFill 4"); static Timer timer4("Meshing2::BlockFill 4");
NgProfiler::RegionTimer reg (timer); RegionTimer reg (timer);
NgProfiler::StartTimer (timer1); timer1.Start();
double filldist = mp.filldist; double filldist = mp.filldist;
cout << "blockfill local h" << endl; PrintMessage (6, "blockfill local h");
cout << "rel filldist = " << filldist << endl;
PrintMessage (3, "blockfill local h");
NgArray<Point<3> > npoints; NgArray<Point<3> > npoints;
@ -95,15 +403,12 @@ namespace netgen
} }
cout << "bbox = " << bbox << endl;
// Point<3> mpc = bbox.Center(); // Point<3> mpc = bbox.Center();
bbox.Increase (bbox.Diam()/2); bbox.Increase (bbox.Diam()/2);
Box<3> meshbox = bbox; Box<3> meshbox = bbox;
NgProfiler::StopTimer (timer1); timer1.Stop();
NgProfiler::StartTimer (timer2); timer2.Start();
LocalH loch2 (bbox, 1, 2); LocalH loch2 (bbox, 1, 2);
@ -147,15 +452,17 @@ namespace netgen
} }
while (changed); while (changed);
NgProfiler::StopTimer (timer2); timer2.Stop();
NgProfiler::StartTimer (timer3); timer3.Start();
if (debugparam.slowchecks) if (debugparam.slowchecks)
(*testout) << "Blockfill with points: " << endl; {
*testout << "loch = " << mesh.LocalHFunction() << endl; (*testout) << "Blockfill with points: " << endl;
*testout << "loch = " << mesh.LocalHFunction() << endl;
*testout << "npoints = " << endl << npoints << endl; *testout << "npoints = " << endl << npoints << endl;
}
for (int i = 1; i <= npoints.Size(); i++) for (int i = 1; i <= npoints.Size(); i++)
{ {
@ -179,8 +486,8 @@ namespace netgen
} }
} }
NgProfiler::StopTimer (timer3); timer3.Stop();
NgProfiler::StartTimer (timer4); timer4.Start();
// find outer points // find outer points
@ -224,7 +531,7 @@ namespace netgen
} }
} }
NgProfiler::StopTimer (timer4); timer4.Stop();
} }
@ -232,191 +539,302 @@ namespace netgen
void Meshing2 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp) void Meshing2 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp)
{ {
static int timer = NgProfiler::CreateTimer ("Meshing2::Delaunay - total"); static Timer timer("Meshing2::Delaunay");
static int timerstart = NgProfiler::CreateTimer ("Meshing2::Delaunay - start"); static Timer timer_addpoints("add points");
static int timerfinish = NgProfiler::CreateTimer ("Meshing2::Delaunay - finish"); RegionTimer reg (timer);
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);
PrintMessage (4, "2D Delaunay meshing");
auto first_point_blockfill = mesh.Points().Range().Next();
cout << "2D Delaunay meshing (in progress)" << endl;
BlockFillLocalH (mesh, mp); BlockFillLocalH (mesh, mp);
NgProfiler::StartTimer (timerstart); auto last_point_blockfill = mesh.Points().Range().Next();
// do the delaunay // Bounding box for starting trig in delaunay
Box<2> bbox (Box<2>::EMPTY_BOX);
// face bounding box:
Box<3> bbox (Box<3>::EMPTY_BOX);
for (int i = 0; i < adfront.GetNFL(); i++) for (int i = 0; i < adfront.GetNFL(); i++)
{ {
const FrontLine & line = adfront.GetLine(i); const FrontLine & line = adfront.GetLine(i);
bbox.Add (Point<3> (adfront.GetPoint (line.L()[0]))); bbox.Add (P2(Point<3> (adfront.GetPoint (line.L()[0]))));
bbox.Add (Point<3> (adfront.GetPoint (line.L()[1]))); bbox.Add (P2(Point<3> (adfront.GetPoint (line.L()[1]))));
} }
for (PointIndex pi : Range(first_point_blockfill, last_point_blockfill))
bbox.Add(P2(mesh[pi]));
for (int i = 0; i < mesh.LockedPoints().Size(); i++) for (int i = 0; i < mesh.LockedPoints().Size(); i++)
bbox.Add (mesh.Point (mesh.LockedPoints()[i])); bbox.Add (P2(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);
NgArray<DelaunayTrig> tempels;
startel.CalcCenter (mesh);
tempels.Append (startel);
searchtree.Insert(startel.BoundingBox(), 0);
NgArray<int> closeels;
NgArray<int> intersecting;
NgArray<INDEX_2> edges;
// reorder points
NgArray<PointIndex, PointIndex::BASE, PointIndex> mixed(old_points.Size());
int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
int prim;
Array<PointIndex> old_points;
BitArray add_point(mesh.Points().Size()+1);
add_point.Clear();
for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
{ {
int i = 0; const auto & s = mesh[si];
while (old_points.Size() % prims[i] == 0) i++; if ( s.domin==domainnr || s.domout==domainnr )
prim = prims[i]; {
add_point.SetBit(s[0]);
add_point.SetBit(s[1]);
}
} }
for (PointIndex pi : old_points) Mesh tempmesh;
mixed[pi] = PointIndex ( (prim * pi) % old_points.Size() + PointIndex::BASE ); tempmesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
tempmesh.AddFaceDescriptor (FaceDescriptor (2, 1, 0, 0));
tempmesh.AddFaceDescriptor (FaceDescriptor (3, 1, 0, 0));
NgProfiler::StopTimer (timerstart); Array<PointIndex, PointIndex> compress;
NgProfiler::StartTimer (timer1); Array<PointIndex, PointIndex> icompress(mesh.Points().Size());
for(auto pi : mesh.Points().Range())
for (PointIndex i1 : old_points) if(add_point.Test(pi))
{ {
PointIndex i = mixed[i1]; icompress[pi] = tempmesh.AddPoint(mesh[pi]);
compress.Append(pi);
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); for (PointIndex pi : Range(first_point_blockfill, last_point_blockfill))
NgProfiler::StartTimer (timerfinish);
for (DelaunayTrig & trig : tempels)
{ {
if (trig[0] < 0) continue; icompress[pi] = tempmesh.AddPoint(mesh[pi]);
compress.Append(pi);
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()) // DelaunayMesh adds surrounding trig (don't add the last 3 points to delaunay AGAIN!
*testout << pi << ": " << mesh[pi].Type() << endl; auto tempmesh_points = tempmesh.Points().Range();
NgProfiler::StopTimer (timerfinish); DelaunayMesh dmesh(tempmesh, bbox);
timer_addpoints.Start();
// // reorder points
// NgArray<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 );
for (auto pi : tempmesh_points)
dmesh.AddPoint(pi);
timer_addpoints.Stop();
for (auto seg : mesh.LineSegments())
{
if ( seg.domin == domainnr || seg.domout == domainnr )
{
if(seg.domin==domainnr)
seg.domout = 0;
if(seg.domout==domainnr)
seg.domin = 0;
seg[0] = icompress[seg[0]];
seg[1] = icompress[seg[1]];
tempmesh.AddSegment(seg);
}
}
for (auto & trig : dmesh.GetElements())
{
if (trig[0] < 0) continue;
Element2d el(trig[0], trig[1], trig[2]);
el.SetIndex (1);
tempmesh.AddSurfaceElement (el);
}
bool conforming = false;
while(!conforming)
{
conforming = true;
BitArray marked_points(tempmesh.Points().Size()+1);
marked_points = false;
// Check for trigs cutting a boundary edge (non-conforming mesh)
auto point_to_trigs = tempmesh.CreatePoint2SurfaceElementTable( 0 );
for (auto & seg : tempmesh.LineSegments())
{
int count_adjacent = 0;;
PointIndex pi0 = seg[0];
PointIndex pi1 = seg[1];
if(marked_points.Test(pi0)) continue;
if(marked_points.Test(pi1)) continue;
for(auto sei : point_to_trigs[pi0])
for( auto i : Range(3))
if(tempmesh[sei][i] == pi1)
count_adjacent++;
if(count_adjacent==2)
continue;
PointIndex pi2;
PointIndex pi3;
ArrayMem<SurfaceElementIndex, 2> cutting_trigs;
for(auto sei : point_to_trigs[pi0])
{
auto & el = tempmesh[sei];
pi2 = el[0] == pi0 ? el[1] : el[0];
pi3 = el[2] == pi0 ? el[1] : el[2];
double alpha, beta;
auto itype = intersect( P2(tempmesh[pi0]), P2(tempmesh[pi1]), P2(tempmesh[pi2]), P2(tempmesh[pi3]), alpha, beta );
if(itype == X_INTERSECTION)
{
cutting_trigs.Append(sei);
break;
}
}
if(cutting_trigs.Size()==0)
continue;
for(auto sei : point_to_trigs[pi2])
{
if(sei==cutting_trigs[0])
continue;
for(auto i : IntRange(3))
if(tempmesh[sei][i]==pi3)
cutting_trigs.Append(sei);
}
// Found two trigs cutting a boundary edge -> perform swap
if(cutting_trigs.Size()==2)
{
conforming = false;
if(marked_points.Test(pi2)) continue;
if(marked_points.Test(pi3)) continue;
auto & el0 = tempmesh[cutting_trigs[0]];
auto & el1 = tempmesh[cutting_trigs[1]];
pi1 = el1[0]+el1[1]+el1[2] - pi2-pi3;
if(marked_points.Test(pi1)) continue;
marked_points.SetBit(pi0);
marked_points.SetBit(pi1);
marked_points.SetBit(pi2);
marked_points.SetBit(pi3);
el0[0] = pi2;
el0[1] = pi1;
el0[2] = pi0;
el1[0] = pi3;
el1[1] = pi0;
el1[2] = pi1;
}
}
}
auto point_to_trigs = tempmesh.CreatePoint2SurfaceElementTable( 0 );
// Mark edges and trigs as inside or outside, starting with boundary edges
enum POSITION { UNKNOWN, BOUNDARY, INSIDE, OUTSIDE };
Array<POSITION, SurfaceElementIndex> trig_pos(tempmesh.SurfaceElements().Size());
ngcore::ClosedHashTable<INT<2>, POSITION> edge_pos(3*tempmesh.SurfaceElements().Size());
trig_pos = UNKNOWN;
for (auto & seg : tempmesh.LineSegments())
{
ArrayMem<SurfaceElementIndex, 2> els;
INT<2> edge{seg[0], seg[1]};
edge.Sort();
edge_pos[edge] = BOUNDARY;
for(auto sei : point_to_trigs[seg[0]])
for( auto i : Range(3))
if(tempmesh[sei][i] == seg[1])
els.Append(sei);
for(auto sei : els)
{
auto & el = tempmesh[sei];
PointIndex pi2 = el[0]+el[1]+el[2] - seg[0] - seg[1];
bool is_left = ::netgen::Area(P2(tempmesh[seg[0]]), P2(tempmesh[seg[1]]), P2(tempmesh[pi2]))>0.0;
POSITION pos;
if(is_left == (seg.domin==domainnr))
pos = INSIDE;
else
pos = OUTSIDE;
INT<2> e1{seg[0], pi2};
INT<2> e2{seg[1], pi2};
e1.Sort();
e2.Sort();
if(!edge_pos.Used(e1))
edge_pos[e1] = pos;
if(!edge_pos.Used(e2))
edge_pos[e2] = pos;
trig_pos[sei] = pos;
}
}
// Advance from boundary edges/trigs to all others
bool have_unknown_trigs = true;
while(have_unknown_trigs)
{
have_unknown_trigs = false;
for (auto sei : Range(tempmesh.SurfaceElements()))
{
auto & el = tempmesh[sei];
if(trig_pos[sei] == UNKNOWN)
{
have_unknown_trigs = true;
// any edge of unkown trig already marked?
for(auto i : IntRange(3))
{
INT<2> edge{el[(i+1)%3], el[(i+2)%3]};
edge.Sort();
if(edge_pos.Used(edge) && edge_pos[edge]!=BOUNDARY)
{
trig_pos[sei] = edge_pos[edge];
break;
}
}
}
// if we could mark the trig -> also mark all edges
if(trig_pos[sei] != UNKNOWN)
for(auto i : IntRange(3))
{
INT<2> edge{el[(i+1)%3], el[(i+2)%3]};
edge.Sort();
if(!edge_pos.Used(edge) || edge_pos[edge]==BOUNDARY)
edge_pos[edge] = trig_pos[sei];
}
}
}
// add inside trigs to actual mesh
for (auto sei : Range(tempmesh.SurfaceElements()))
{
if(trig_pos[sei] == INSIDE)
{
auto el = tempmesh[sei];
Vec<3> n = Cross (tempmesh[el[1]]-tempmesh[el[0]],
tempmesh[el[2]]-tempmesh[el[0]]);
if (n(2) < 0) Swap (el[1], el[2]);
el[0] = compress[el[0]];
el[1] = compress[el[1]];
el[2] = compress[el[2]];
el.SetIndex(domainnr);
mesh.AddSurfaceElement(el);
}
}
mesh.Compress();
} }
} }

4
libsrc/meshing/meshtype.hpp
View File

@ -1263,8 +1263,10 @@ namespace netgen
bool uselocalh = true; bool uselocalh = true;
/// grading for local h /// grading for local h
double grading = 0.3; double grading = 0.3;
/// use delaunay meshing /// use delaunay for 3d meshing
bool delaunay = true; bool delaunay = true;
/// use delaunay for 2d meshing
bool delaunay2d = false;
/// maximal mesh size /// maximal mesh size
double maxh = 1e10; double maxh = 1e10;
/// minimal mesh size /// minimal mesh size

5
libsrc/meshing/python_mesh.hpp
View File

@ -48,6 +48,9 @@ filldist: float = 0.1
delaunay: bool = True delaunay: bool = True
Use delaunay meshing. Use delaunay meshing.
delaunay2d : bool = True
Use delaunay meshing for 2d geometries.
Optimization Parameters Optimization Parameters
----------------------- -----------------------
@ -109,6 +112,8 @@ inline void CreateMPfromKwargs(MeshingParameters& mp, py::kwargs kwargs, bool th
mp.grading = py::cast<double>(kwargs.attr("pop")("grading")); mp.grading = py::cast<double>(kwargs.attr("pop")("grading"));
if(kwargs.contains("delaunay")) if(kwargs.contains("delaunay"))
mp.delaunay = py::cast<bool>(kwargs.attr("pop")("delaunay")); mp.delaunay = py::cast<bool>(kwargs.attr("pop")("delaunay"));
if(kwargs.contains("delaunay2d"))
mp.delaunay2d = py::cast<bool>(kwargs.attr("pop")("delaunay2d"));
if(kwargs.contains("maxh")) if(kwargs.contains("maxh"))
mp.maxh = py::cast<double>(kwargs.attr("pop")("maxh")); mp.maxh = py::cast<double>(kwargs.attr("pop")("maxh"));
if(kwargs.contains("minh")) if(kwargs.contains("minh"))