diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 6d00f1fe..905f09b6 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -392,10 +392,16 @@ namespace netgen shared_ptr & mesh, 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"); Box<2> bbox = geometry.GetBoundingBox (); + t_h.Start(); if (bbox.Diam() < mp.maxh) mp.maxh = bbox.Diam(); @@ -412,11 +418,14 @@ namespace netgen + t_part_boundary.Start(); geometry.PartitionBoundary (mp, mp.maxh, *mesh); + t_part_boundary.Stop(); PrintMessage (3, "Boundary mesh done, np = ", mesh->GetNP()); + t_hpref.Start(); // marks mesh points for hp-refinement for (int i = 0; i < geometry.GetNP(); i++) if (geometry.GetPoint(i).hpref) @@ -434,6 +443,7 @@ namespace netgen } (*mesh)[mpi].Singularity(geometry.GetPoint(i).hpref); } + t_hpref.Stop(); int maxdomnr = 0; @@ -459,6 +469,7 @@ namespace netgen mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) ); mesh->CalcLocalH(mp.grading); + t_h.Stop(); int bnp = mesh->GetNP(); // boundary points auto BndPntRange = mesh->Points().Range(); @@ -469,6 +480,7 @@ namespace netgen for (int domnr = 1; domnr <= maxdomnr; domnr++) if (geometry.GetDomainTensorMeshing (domnr)) { // tensor product mesh + RegionTimer rt(t_tensor); NgArray nextpi(bnp); NgArray si1(bnp), si2(bnp); @@ -556,8 +568,10 @@ namespace netgen + static Timer t_domain("Mesh domain"); for (int domnr = 1; domnr <= maxdomnr; domnr++) { + RegionTimer rt(t_domain); if (geometry.GetDomainTensorMeshing (domnr)) continue; double h = mp.maxh; @@ -573,16 +587,26 @@ namespace netgen Meshing2 meshing (geometry, mp, Box<3> (pmin, pmax)); - NgArray compress(bnp); + NgArray compress(mesh->GetNP()); compress = -1; int cnt = 0; - for (PointIndex pi : BndPntRange) - if ( (*mesh)[pi].GetLayer() == geometry.GetDomainLayer(domnr)) - { - meshing.AddPoint ((*mesh)[pi], pi); - cnt++; - compress[pi] = cnt; - } + + for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++) + { + const auto & s = (*mesh)[si]; + if ( s.domin==domnr || s.domout==domnr ) + { + for (auto pi : {s[0], s[1]}) + { + if(compress[pi]==-1) + { + meshing.AddPoint((*mesh)[pi], pi); + cnt++; + compress[pi] = cnt; + } + } + } + } PointGeomInfo gi; gi.trignum = 1; @@ -600,12 +624,15 @@ namespace netgen } } - // not complete, use at own risk ... - // meshing.Delaunay(*mesh, domnr, mp); - mp.checkoverlap = 0; - auto res = meshing.GenerateMesh (*mesh, mp, h, domnr); - if (res != 0) - throw NgException("meshing failed"); + if(mp.delaunay2d && cnt>100) + meshing.Delaunay(*mesh, domnr, mp); + else + { + // mp.checkoverlap = 0; + auto res = meshing.GenerateMesh (*mesh, mp, h, domnr); + if (res != 0) + throw NgException("meshing failed"); + } for (SurfaceElementIndex sei = oldnf; sei < mesh->GetNSE(); sei++) (*mesh)[sei].SetIndex (domnr); diff --git a/libsrc/meshing/delaunay2d.cpp b/libsrc/meshing/delaunay2d.cpp index a61016e2..66b22708 100644 --- a/libsrc/meshing/delaunay2d.cpp +++ b/libsrc/meshing/delaunay2d.cpp @@ -1,34 +1,52 @@ #include #include "meshing.hpp" +#include + + // not yet working .... 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 { PointIndex pnums[3]; - Point<3> c; + Point<2> c; + public: double r; double rad2; - public: - DelaunayTrig () { ; } + DelaunayTrig () = default; 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]; } 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; + Point<2> p1 = P2(mesh[pnums[0]]); + Point<2> p2 = P2(mesh[pnums[1]]); + Point<2> p3 = P2(mesh[pnums[2]]); + Vec<2> v1 = p2-p1; + Vec<2> v2 = p3-p1; Mat<2,2> mat, inv; mat(0,0) = v1*v1; mat(0,1) = v1*v2; @@ -45,9 +63,301 @@ namespace netgen r = sqrt(rad2); } - Point<3> Center() const { return c; } + Point<2> 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)); } + Box<2> BoundingBox() const { return Box<2> (c-Vec<2>(r,r), c+Vec<2>(r,r)); } + }; + + class DelaunayMesh + { + ngcore::ClosedHashTable, INT<2>> edge_to_trig; + Array trigs; + unique_ptr> tree; + Mesh & mesh; + + Array closeels; + Array intersecting; + Array> edges; + + int GetNeighbour( int eli, int edge ) + { + auto p0 = trigs[eli][(edge+1)%3]; + auto p1 = trigs[eli][(edge+2)%3]; + if(p1 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 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 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>(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 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 & GetElements() { return trigs; } + }; ostream & operator<< (ostream & ost, DelaunayTrig trig) @@ -59,20 +369,18 @@ namespace netgen 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); + static Timer timer("Meshing2::BlockFill"); + static Timer timer1("Meshing2::BlockFill 1"); + static Timer timer2("Meshing2::BlockFill 2"); + static Timer timer3("Meshing2::BlockFill 3"); + static Timer timer4("Meshing2::BlockFill 4"); + RegionTimer reg (timer); - NgProfiler::StartTimer (timer1); + timer1.Start(); double filldist = mp.filldist; - cout << "blockfill local h" << endl; - cout << "rel filldist = " << filldist << endl; - PrintMessage (3, "blockfill local h"); + PrintMessage (6, "blockfill local h"); NgArray > npoints; @@ -95,15 +403,12 @@ namespace netgen } - cout << "bbox = " << bbox << endl; - - // Point<3> mpc = bbox.Center(); bbox.Increase (bbox.Diam()/2); Box<3> meshbox = bbox; - NgProfiler::StopTimer (timer1); - NgProfiler::StartTimer (timer2); + timer1.Stop(); + timer2.Start(); LocalH loch2 (bbox, 1, 2); @@ -147,15 +452,17 @@ namespace netgen } while (changed); - NgProfiler::StopTimer (timer2); - NgProfiler::StartTimer (timer3); + timer2.Stop(); + timer3.Start(); if (debugparam.slowchecks) - (*testout) << "Blockfill with points: " << endl; - *testout << "loch = " << mesh.LocalHFunction() << endl; - - *testout << "npoints = " << endl << npoints << endl; + { + (*testout) << "Blockfill with points: " << endl; + *testout << "loch = " << mesh.LocalHFunction() << endl; + + *testout << "npoints = " << endl << npoints << endl; + } for (int i = 1; i <= npoints.Size(); i++) { @@ -179,8 +486,8 @@ namespace netgen } } - NgProfiler::StopTimer (timer3); - NgProfiler::StartTimer (timer4); + timer3.Stop(); + timer4.Start(); // 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) { - 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); + static Timer timer("Meshing2::Delaunay"); + static Timer timer_addpoints("add points"); + RegionTimer reg (timer); + PrintMessage (4, "2D Delaunay meshing"); - - cout << "2D Delaunay meshing (in progress)" << endl; - + auto first_point_blockfill = mesh.Points().Range().Next(); BlockFillLocalH (mesh, mp); - NgProfiler::StartTimer (timerstart); + auto last_point_blockfill = mesh.Points().Range().Next(); - // do the delaunay - - - // face bounding box: - Box<3> bbox (Box<3>::EMPTY_BOX); + // Bounding box for starting trig in delaunay + Box<2> bbox (Box<2>::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]))); + bbox.Add (P2(Point<3> (adfront.GetPoint (line.L()[0])))); + 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++) - 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 tempels; - startel.CalcCenter (mesh); - - tempels.Append (startel); - searchtree.Insert(startel.BoundingBox(), 0); - - NgArray closeels; - NgArray intersecting; - NgArray edges; - - - - - // reorder points - NgArray mixed(old_points.Size()); - int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 }; - int prim; - + Array old_points; + BitArray add_point(mesh.Points().Size()+1); + add_point.Clear(); + for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++) { - int i = 0; - while (old_points.Size() % prims[i] == 0) i++; - prim = prims[i]; + const auto & s = mesh[si]; + if ( s.domin==domainnr || s.domout==domainnr ) + { + add_point.SetBit(s[0]); + add_point.SetBit(s[1]); + } } - for (PointIndex pi : old_points) - mixed[pi] = PointIndex ( (prim * pi) % old_points.Size() + PointIndex::BASE ); + Mesh tempmesh; + tempmesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0)); + tempmesh.AddFaceDescriptor (FaceDescriptor (2, 1, 0, 0)); + tempmesh.AddFaceDescriptor (FaceDescriptor (3, 1, 0, 0)); + + Array compress; + Array icompress(mesh.Points().Size()); + + for(auto pi : mesh.Points().Range()) + if(add_point.Test(pi)) + { + icompress[pi] = tempmesh.AddPoint(mesh[pi]); + compress.Append(pi); + } + + for (PointIndex pi : Range(first_point_blockfill, last_point_blockfill)) + { + icompress[pi] = tempmesh.AddPoint(mesh[pi]); + compress.Append(pi); + } + + // DelaunayMesh adds surrounding trig (don't add the last 3 points to delaunay AGAIN! + auto tempmesh_points = tempmesh.Points().Range(); + + DelaunayMesh dmesh(tempmesh, bbox); + + timer_addpoints.Start(); + +// // reorder points +// NgArray 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 (auto pi : tempmesh_points) + dmesh.AddPoint(pi); + + timer_addpoints.Stop(); - for (PointIndex i1 : old_points) + for (auto seg : mesh.LineSegments()) + { + if ( seg.domin == domainnr || seg.domout == domainnr ) { - 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); + 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); } + } - NgProfiler::StopTimer (timer1); - NgProfiler::StartTimer (timerfinish); + for (auto & trig : dmesh.GetElements()) + { + if (trig[0] < 0) continue; - for (DelaunayTrig & trig : tempels) + 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()) { - if (trig[0] < 0) continue; + int count_adjacent = 0;; + PointIndex pi0 = seg[0]; + PointIndex pi1 = seg[1]; + if(marked_points.Test(pi0)) continue; + if(marked_points.Test(pi1)) continue; - Point<3> c = Center (mesh[trig[0]], mesh[trig[1]], mesh[trig[2]]); - if (!adfront.Inside (Point<2> (c(0),c(1)))) continue; + for(auto sei : point_to_trigs[pi0]) + for( auto i : Range(3)) + if(tempmesh[sei][i] == pi1) + count_adjacent++; - 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]); + if(count_adjacent==2) + continue; - Element2d el(trig[0], trig[1], trig[2]); - el.SetIndex (domainnr); - mesh.AddSurfaceElement (el); + PointIndex pi2; + PointIndex pi3; + ArrayMem 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; + } } + } - for (PointIndex pi : mesh.Points().Range()) - *testout << pi << ": " << mesh[pi].Type() << endl; + auto point_to_trigs = tempmesh.CreatePoint2SurfaceElementTable( 0 ); - NgProfiler::StopTimer (timerfinish); + // Mark edges and trigs as inside or outside, starting with boundary edges + enum POSITION { UNKNOWN, BOUNDARY, INSIDE, OUTSIDE }; + Array trig_pos(tempmesh.SurfaceElements().Size()); + ngcore::ClosedHashTable, POSITION> edge_pos(3*tempmesh.SurfaceElements().Size()); + trig_pos = UNKNOWN; + + for (auto & seg : tempmesh.LineSegments()) + { + ArrayMem 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(); } } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index be4482d6..2f49ad70 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1263,8 +1263,10 @@ namespace netgen bool uselocalh = true; /// grading for local h double grading = 0.3; - /// use delaunay meshing + /// use delaunay for 3d meshing bool delaunay = true; + /// use delaunay for 2d meshing + bool delaunay2d = false; /// maximal mesh size double maxh = 1e10; /// minimal mesh size diff --git a/libsrc/meshing/python_mesh.hpp b/libsrc/meshing/python_mesh.hpp index 6deab5d2..110c26ee 100644 --- a/libsrc/meshing/python_mesh.hpp +++ b/libsrc/meshing/python_mesh.hpp @@ -48,6 +48,9 @@ filldist: float = 0.1 delaunay: bool = True Use delaunay meshing. +delaunay2d : bool = True + Use delaunay meshing for 2d geometries. + Optimization Parameters ----------------------- @@ -109,6 +112,8 @@ inline void CreateMPfromKwargs(MeshingParameters& mp, py::kwargs kwargs, bool th mp.grading = py::cast(kwargs.attr("pop")("grading")); if(kwargs.contains("delaunay")) mp.delaunay = py::cast(kwargs.attr("pop")("delaunay")); + if(kwargs.contains("delaunay2d")) + mp.delaunay2d = py::cast(kwargs.attr("pop")("delaunay2d")); if(kwargs.contains("maxh")) mp.maxh = py::cast(kwargs.attr("pop")("maxh")); if(kwargs.contains("minh"))