From 39acabe406741eee4733dd132e83517d56299c0e Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 10:46:56 +0200 Subject: [PATCH 1/8] split delaunay postprocessing code into smaller funtions --- libsrc/meshing/delaunay.cpp | 336 +++++++++++++++++------------------- 1 file changed, 160 insertions(+), 176 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index bfed1841..01c69dbb 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -737,137 +737,13 @@ namespace netgen } - - - - - void Meshing3 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp) + void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray & tempels ) { - static Timer t("Meshing3::Delaunay"); RegionTimer reg(t); - - int np, ne; + static Timer tdegenerated("Delaunay - remove degenerated"); RegionTimer rt(tdegenerated); - PrintMessage (1, "Delaunay meshing"); - PrintMessage (3, "number of points: ", mesh.GetNP()); - PushStatus ("Delaunay meshing"); + auto np = points.Size(); - - NgArray tempels; - Point3d pmin, pmax; - - DelaunayTet startel; - - int oldnp = mesh.GetNP(); - if (mp.blockfill) - { - BlockFillLocalH (mesh, mp); - PrintMessage (3, "number of points: ", mesh.GetNP()); - } - - 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"); - - { - RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0); - for (int i = 1; i <= 4; i++) - { - tempmesh.FindOpenElements (); - - PrintMessage (5, "Num open: ", tempmesh.GetNOpenElements()); - tempmesh.CalcSurfacesOfNode (); - - tempmesh.FreeOpenElementsEnvironment (1); - - MeshOptimize3d meshopt(mp); - // tempmesh.CalcSurfacesOfNode(); - meshopt.SwapImprove(tempmesh, OPT_CONFORM); - } - } - - MeshQuality3d (tempmesh); - - tempels.SetSize(tempmesh.GetNE()); - tempels.SetSize(0); - for (auto & el : tempmesh.VolumeElements()) - tempels.Append (el); - } - - - - // remove degenerated - static Timer tdegenerated("Delaunay - remove degenerated"); - tdegenerated.Start(); - - NgBitArray badnode(mesh.GetNP()); + NgBitArray badnode(np); badnode.Clear(); int ndeg = 0; for (int i = 1; i <= tempels.Size(); i++) @@ -876,10 +752,10 @@ namespace netgen for (int j = 0; j < 4; j++) el[j] = tempels.Elem(i)[j]; // Element & el = tempels.Elem(i); - const Point3d & lp1 = mesh.Point (el[0]); - const Point3d & lp2 = mesh.Point (el[1]); - const Point3d & lp3 = mesh.Point (el[2]); - const Point3d & lp4 = mesh.Point (el[3]); + 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); @@ -903,7 +779,7 @@ namespace netgen Swap (el[2], el[3]); } - ne = tempels.Size(); + auto ne = tempels.Size(); for (int i = ne; i >= 1; i--) { const DelaunayTet & el = tempels.Get(i); @@ -916,11 +792,16 @@ namespace netgen PrintMessage (3, ndeg, " degenerated elements removed"); - tdegenerated.Stop(); + } + // Remove flat tets containing two adjacent surface trigs + NgArray DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray & tempels ) + { static Timer topenel("Delaunay - find openel"); topenel.Start(); + auto np = mesh.GetNP(); + // find surface triangles which are no face of any tet INDEX_3_HASHTABLE openeltab(mesh.GetNOpenElements()+3); @@ -1010,6 +891,8 @@ namespace netgen // cout << "tetedges:"; // tetedges.PrintMemInfo (cout); + NgBitArray badnode(np); + badnode.Clear(); for (INDEX_2_HASHTABLE::Iterator it = twotrias.Begin(); it != twotrias.End(); it++) @@ -1038,36 +921,7 @@ namespace netgen } } - /* - for (i = 1; i <= twotrias.GetNBags(); i++) - for (j = 1; j <= twotrias.GetBagSize (i); j++) - { - INDEX_2 hi2, hi3; - twotrias.GetData (i, j, hi2, hi3); - hi3.Sort(); - if (tetedges.Used (hi3)) - { - const Point3d & p1 = mesh.Point (hi2.I1()); - const Point3d & p2 = mesh.Point (hi2.I2()); - const Point3d & p3 = mesh.Point (hi3.I1()); - const Point3d & p4 = mesh.Point (hi3.I2()); - Vec3d v1(p1, p2); - Vec3d v2(p1, p3); - Vec3d v3(p1, p4); - Vec3d n = Cross (v1, v2); - double vol = n * v3; - - double h = v1.Length() + v2.Length() + v3.Length(); - if (fabs (vol) < 1e-4 * (h * h * h)) // old: 1e-12 - { - badnode.Set(hi3.I1()); - badnode.Set(hi3.I2()); - } - } - } - */ - - ne = tempels.Size(); + auto ne = tempels.Size(); for (int i = ne; i >= 1; i--) { const DelaunayTet & el = tempels.Get(i); @@ -1081,9 +935,12 @@ namespace netgen topenel.Stop(); - static Timer trem_intersect("Delaunay - remove intersecting"); - trem_intersect.Start(); + return openels; + } + void DelaunayRemoveIntersecting( const Mesh & mesh, NgArray & tempels, NgArray & openels, Point3d pmin, Point3d pmax ) + { + static Timer trem_intersect("Delaunay - remove intersecting"); RegionTimert(trem_intersect); // find intersecting: PrintMessage (3, "Remove intersecting"); @@ -1198,12 +1055,11 @@ namespace netgen } } } - + } - trem_intersect.Stop(); - - static Timer trem_outer("Delaunay - remove outer"); - trem_outer.Start(); + void DelaunayRemoveOuter( const Mesh & mesh, NgArray & tempels, AdFront3 * adfront ) + { + static Timer trem_outer("Delaunay - remove outer"); RegionTimer rt(trem_outer); PrintMessage (3, "Remove outer"); @@ -1393,7 +1249,7 @@ namespace netgen PrintMessage (5, "tables filled"); - ne = tempels.Size(); + auto ne = tempels.Size(); NgBitArray inner(ne), outer(ne); inner.Clear(); outer.Clear(); @@ -1649,6 +1505,138 @@ namespace netgen // 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); + + int np, ne; + + PrintMessage (1, "Delaunay meshing"); + PrintMessage (3, "number of points: ", mesh.GetNP()); + PushStatus ("Delaunay meshing"); + + + NgArray tempels; + Point3d pmin, pmax; + + DelaunayTet startel; + + int oldnp = mesh.GetNP(); + if (mp.blockfill) + { + BlockFillLocalH (mesh, mp); + PrintMessage (3, "number of points: ", mesh.GetNP()); + } + + 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 (); + RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0); + 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); + auto openels = DelaunayRemoveTwoTriaTets(mesh, tempels); + DelaunayRemoveIntersecting(mesh, tempels, openels, pmin, pmax); + DelaunayRemoveOuter(mesh, tempels, adfront); + for (int i = 0; i < tempels.Size(); i++) { Element el(4); @@ -1657,10 +1645,6 @@ namespace netgen mesh.AddVolumeElement (el); } - PrintMessage (5, "outer removed"); - - trem_outer.Stop(); - mesh.FindOpenElements(domainnr); mesh.Compress(); From 17af3d009117bc4516a247a92abc7ef48bb56976 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 10:51:43 +0200 Subject: [PATCH 2/8] Timers, cleanup in delaunay --- libsrc/meshing/delaunay.cpp | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 01c69dbb..df1a7902 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -795,17 +795,16 @@ namespace netgen } // Remove flat tets containing two adjacent surface trigs - NgArray DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray & tempels ) + void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray & tempels, NgArray & openels ) { - static Timer topenel("Delaunay - find openel"); - topenel.Start(); + static Timer topenel("Delaunay - find openel"); RegionTimer rt(topenel); auto np = mesh.GetNP(); // find surface triangles which are no face of any tet INDEX_3_HASHTABLE openeltab(mesh.GetNOpenElements()+3); - NgArray openels; + openels.SetSize(0); for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & tri = mesh.OpenElement(i); @@ -837,9 +836,6 @@ namespace netgen } - - - // find open triangle with close edge (from halfening of surface squares) INDEX_2_HASHTABLE twotrias(mesh.GetNOpenElements()+5); @@ -931,16 +927,11 @@ namespace netgen badnode.Test(el[3]) ) tempels.DeleteElement(i); } - - - topenel.Stop(); - - return openels; } void DelaunayRemoveIntersecting( const Mesh & mesh, NgArray & tempels, NgArray & openels, Point3d pmin, Point3d pmax ) { - static Timer trem_intersect("Delaunay - remove intersecting"); RegionTimert(trem_intersect); + static Timer trem_intersect("Delaunay - remove intersecting"); RegionTimer rt(trem_intersect); // find intersecting: PrintMessage (3, "Remove intersecting"); @@ -1633,7 +1624,9 @@ namespace netgen DelaunayRemoveDegenerated(mesh.Points(), tempels); - auto openels = DelaunayRemoveTwoTriaTets(mesh, tempels); + + NgArray openels; + DelaunayRemoveTwoTriaTets(mesh, tempels, openels); DelaunayRemoveIntersecting(mesh, tempels, openels, pmin, pmax); DelaunayRemoveOuter(mesh, tempels, adfront); From 7623289c278dce073ed33068c0a199c2731d629b Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 10:47:24 +0200 Subject: [PATCH 3/8] Timer to AdFront3::Inside --- libsrc/meshing/adfront3.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/libsrc/meshing/adfront3.cpp b/libsrc/meshing/adfront3.cpp index 31e1248a..588f1d23 100644 --- a/libsrc/meshing/adfront3.cpp +++ b/libsrc/meshing/adfront3.cpp @@ -787,6 +787,7 @@ void AdFront3 :: SetStartFront (int /* baseelnp */) bool AdFront3 :: Inside (const Point<3> & p) const { + static Timer timer("AdFront3::Inside"); RegionTimer rt(timer); int cnt; Vec3d n, v1, v2; DenseMatrix a(3), ainv(3); From 6c37ce33b0c932b79d097292b43001958d0736fa Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 11:23:37 +0200 Subject: [PATCH 4/8] CreatePoint2ElementTable with optional points bitarray --- libsrc/meshing/delaunay.cpp | 14 +++++++------- libsrc/meshing/meshclass.cpp | 34 +++++++++++++++++++++++++++------- libsrc/meshing/meshclass.hpp | 2 +- 3 files changed, 35 insertions(+), 15 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index df1a7902..88c6cd7e 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -1598,20 +1598,20 @@ namespace netgen // tempmesh.Save ("tempmesh.vol"); { - MeshOptimize3d meshopt(mp); - tempmesh.Compress(); - tempmesh.FindOpenElements (); RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0); - for (auto i : Range(10)) + for (int i = 1; i <= 4; i++) { + tempmesh.FindOpenElements (); + PrintMessage (5, "Num open: ", tempmesh.GetNOpenElements()); + tempmesh.CalcSurfacesOfNode (); - if(i%5==0) - tempmesh.FreeOpenElementsEnvironment (1); + tempmesh.FreeOpenElementsEnvironment (1); + MeshOptimize3d meshopt(mp); + // tempmesh.CalcSurfacesOfNode(); meshopt.SwapImprove(tempmesh, OPT_CONFORM); } - tempmesh.Compress(); } MeshQuality3d (tempmesh); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 436b3c0b..c9b5c7ed 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -6545,14 +6545,34 @@ namespace netgen } - Table Mesh :: CreatePoint2ElementTable() const + Table Mesh :: CreatePoint2ElementTable(std::optional points) const { - return ngcore::CreateSortedTable( volelements.Range(), - [&](auto & table, ElementIndex ei) - { - for (PointIndex pi : (*this)[ei].PNums()) - table.Add (pi, ei); - }, GetNP()); + if(points) + { + const auto & free_points = *points; + return ngcore::CreateSortedTable( volelements.Range(), + [&](auto & table, ElementIndex ei) + { + const auto & el = (*this)[ei]; + if(el.IsDeleted()) + return; + + for (PointIndex pi : el.PNums()) + if(free_points[pi]) + table.Add (pi, ei); + }, GetNP()); + } + else + return ngcore::CreateSortedTable( volelements.Range(), + [&](auto & table, ElementIndex ei) + { + const auto & el = (*this)[ei]; + if(el.IsDeleted()) + return; + + for (PointIndex pi : el.PNums()) + table.Add (pi, ei); + }, GetNP()); } Table Mesh :: CreatePoint2SurfaceElementTable( int faceindex ) const diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index e3e61477..3c3311ed 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -772,7 +772,7 @@ namespace netgen - Table CreatePoint2ElementTable() const; + Table CreatePoint2ElementTable(std::optional points = std::nullopt) const; Table CreatePoint2SurfaceElementTable( int faceindex=0 ) const; DLL_HEADER bool PureTrigMesh (int faceindex = 0) const; From a2cc1028493272ca1f267d70a2f2c31524fdbc5b Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 13:31:44 +0200 Subject: [PATCH 5/8] delaunay - stay consistent with code on master --- libsrc/meshing/delaunay.cpp | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 88c6cd7e..045fd1ba 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -737,14 +737,10 @@ namespace netgen } - void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray & tempels ) + void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray & tempels, NgBitArray & badnode, int np ) { static Timer tdegenerated("Delaunay - remove degenerated"); RegionTimer rt(tdegenerated); - auto np = points.Size(); - - NgBitArray badnode(np); - badnode.Clear(); int ndeg = 0; for (int i = 1; i <= tempels.Size(); i++) { @@ -795,12 +791,10 @@ namespace netgen } // Remove flat tets containing two adjacent surface trigs - void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray & tempels, NgArray & openels ) + void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray & tempels, NgArray & openels, NgBitArray & badnode, int np ) { static Timer topenel("Delaunay - find openel"); RegionTimer rt(topenel); - auto np = mesh.GetNP(); - // find surface triangles which are no face of any tet INDEX_3_HASHTABLE openeltab(mesh.GetNOpenElements()+3); @@ -887,9 +881,6 @@ namespace netgen // cout << "tetedges:"; // tetedges.PrintMemInfo (cout); - NgBitArray badnode(np); - badnode.Clear(); - for (INDEX_2_HASHTABLE::Iterator it = twotrias.Begin(); it != twotrias.End(); it++) { @@ -1622,11 +1613,13 @@ namespace netgen tempels.Append (el); } + NgBitArray badnode(mesh.GetNP()); + badnode.Clear(); - DelaunayRemoveDegenerated(mesh.Points(), tempels); + DelaunayRemoveDegenerated(mesh.Points(), tempels, badnode, np); NgArray openels; - DelaunayRemoveTwoTriaTets(mesh, tempels, openels); + DelaunayRemoveTwoTriaTets(mesh, tempels, openels, badnode, np); DelaunayRemoveIntersecting(mesh, tempels, openels, pmin, pmax); DelaunayRemoveOuter(mesh, tempels, adfront); From 9ddf2424e238d7d895c8ff5bd7face8a331b2eed Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 15:47:50 +0200 Subject: [PATCH 6/8] rewrite of DelaunayRemoveTwoTriaTets (much faster now) --- libsrc/meshing/delaunay.cpp | 230 ++++++++++++++++++++---------------- 1 file changed, 130 insertions(+), 100 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index 045fd1ba..e5abc65a 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -796,119 +796,149 @@ namespace netgen 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 tets_with_3_bnd_points(ne); + atomic cnt = 0; + + // table of tets with >= 2 boundary points, store in extra array tets with >=3 boundary points + auto p2el = ngcore::CreateSortedTable( 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(); - INDEX_3_HASHTABLE openeltab(mesh.GetNOpenElements()+3); openels.SetSize(0); for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & tri = mesh.OpenElement(i); - INDEX_3 i3(tri[0], tri[1], tri[2]); + ngcore::INT<3> i3(tri[0], tri[1], tri[2]); i3.Sort(); - openeltab.Set (i3, i); + if(!face_table.Used(i3)) + openels.Append(i); } - for (int i = 1; i <= tempels.Size(); i++) + t2.Stop(); + + auto p2sel = ngcore::CreateSortedTable( 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()); + + for (auto i : openels) { - for (int j = 0; j < 4; j++) - { - INDEX_3 i3 = tempels.Get(i).GetFace (j); - i3.Sort(); - if (openeltab.Used(i3)) - openeltab.Set (i3, 0); - } - } - - // and store them in openels - for (int i = 1; i <= openeltab.GetNBags(); i++) - for (int j = 1; j <= openeltab.GetBagSize(i); j++) - { - INDEX_3 i3; - int fnr; - openeltab.GetData (i, j, i3, fnr); - if (fnr) - openels.Append (fnr); - } + 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 open triangle with close edge (from halfening of surface squares) - - INDEX_2_HASHTABLE twotrias(mesh.GetNOpenElements()+5); - // for (i = 1; i <= mesh.GetNOpenElements(); i++) - for (int ii = 1; ii <= openels.Size(); ii++) - { - int i = openels.Get(ii); - const Element2d & el = mesh.OpenElement(i); - for (int j = 1; j <= 3; j++) - { - INDEX_2 hi2 (el.PNumMod (j), el.PNumMod(j+1)); - hi2.Sort(); - if (twotrias.Used(hi2)) - { - INDEX_2 hi3; - hi3 = twotrias.Get (hi2); - hi3.I2() = el.PNumMod (j+2); - twotrias.Set (hi2, hi3); - } - else - { - INDEX_2 hi3(el.PNumMod (j+2), 0); - twotrias.Set (hi2, hi3); - } - } + // 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.Set(pi2); + badnode.Set(pi3); + } + break; + } + } + } + } } - INDEX_2_HASHTABLE tetedges(tempels.Size() + 5); - for (int i = 1; i <= tempels.Size(); i++) - { - const DelaunayTet & el = tempels.Get(i); - INDEX_2 i2; - for (int j = 1; j <= 6; j++) - { - switch (j) - { - case 1: i2.I1()=el[0]; i2.I2()=el[1]; break; - case 2: i2.I1()=el[0]; i2.I2()=el[2]; break; - case 3: i2.I1()=el[0]; i2.I2()=el[3]; break; - case 4: i2.I1()=el[1]; i2.I2()=el[2]; break; - case 5: i2.I1()=el[1]; i2.I2()=el[3]; break; - case 6: i2.I1()=el[2]; i2.I2()=el[3]; break; - default: i2.I1()=i2.I2()=0; break; - } - i2.Sort(); - tetedges.Set (i2, 1); - } - } - // cout << "tetedges:"; - // tetedges.PrintMemInfo (cout); - - for (INDEX_2_HASHTABLE::Iterator it = twotrias.Begin(); - it != twotrias.End(); it++) - { - INDEX_2 hi2, hi3; - twotrias.GetData (it, hi2, hi3); - hi3.Sort(); - if (tetedges.Used (hi3)) - { - const Point3d & p1 = mesh.Point ( PointIndex (hi2.I1())); - const Point3d & p2 = mesh.Point ( PointIndex (hi2.I2())); - const Point3d & p3 = mesh.Point ( PointIndex (hi3.I1())); - const Point3d & p4 = mesh.Point ( PointIndex (hi3.I2())); - Vec3d v1(p1, p2); - Vec3d v2(p1, p3); - Vec3d v3(p1, p4); - Vec3d n = Cross (v1, v2); - double vol = n * v3; - - double h = v1.Length() + v2.Length() + v3.Length(); - if (fabs (vol) < 1e-4 * (h * h * h)) // old: 1e-12 - { - badnode.Set(hi3.I1()); - badnode.Set(hi3.I2()); - } - } - } - - auto ne = tempels.Size(); for (int i = ne; i >= 1; i--) { const DelaunayTet & el = tempels.Get(i); From 36d9ead3bc39b87e4aa2ecf7c24dc4905aade034 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 15:57:15 +0200 Subject: [PATCH 7/8] cmake - log output on failure in gitlab-ci --- cmake/SuperBuild.cmake | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cmake/SuperBuild.cmake b/cmake/SuperBuild.cmake index 6e9600df..79f136ab 100644 --- a/cmake/SuperBuild.cmake +++ b/cmake/SuperBuild.cmake @@ -166,6 +166,10 @@ else() set(NETGEN_BUILD_COMMAND ${CMAKE_COMMAND} --build ${CMAKE_CURRENT_BINARY_DIR}/netgen --config ${CMAKE_BUILD_TYPE}) endif() +if(DEFINED ENV{CI} AND WIN32) + set(log_output LOG_BUILD ON LOG_MERGED_STDOUTERR ON LOG_OUTPUT_ON_FAILURE ON) +endif() + ExternalProject_Add (netgen DEPENDS ${NETGEN_DEPENDENCIES} SOURCE_DIR ${PROJECT_SOURCE_DIR} @@ -174,6 +178,7 @@ ExternalProject_Add (netgen BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/netgen BUILD_COMMAND ${NETGEN_BUILD_COMMAND} STEP_TARGETS build + ${log_output} ) # Check if the git submodules (i.e. pybind11) are up to date From ba148e8b3b01f43e06925499195cbdc3ace21cce Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 4 Jun 2021 16:23:18 +0200 Subject: [PATCH 8/8] cleanup, more parallel --- libsrc/meshing/delaunay.cpp | 134 +++++++++++++++++++----------------- 1 file changed, 70 insertions(+), 64 deletions(-) diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index e5abc65a..8bee0775 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -737,10 +737,13 @@ namespace netgen } - void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray & tempels, NgBitArray & badnode, int np ) + void DelaunayRemoveDegenerated( const Mesh::T_POINTS & points, NgArray & 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++) { @@ -791,7 +794,7 @@ namespace netgen } // Remove flat tets containing two adjacent surface trigs - void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray & tempels, NgArray & openels, NgBitArray & badnode, int np ) + void DelaunayRemoveTwoTriaTets( const Mesh & mesh, NgArray & tempels, NgArray & openels ) { static Timer topenel("Delaunay - find openel"); RegionTimer rt(topenel); @@ -876,76 +879,82 @@ namespace netgen table.Add(tri[2], openel_i); }, mesh.GetNP()); - for (auto i : openels) - { - const Element2d & tri = mesh.OpenElement(i); + ngcore::BitArray badnode(mesh.GetNP()+PointIndex::BASE); + badnode.Clear(); - for( auto edge : Range(3) ) - { - auto pi0 = tri[edge]; - auto pi1 = tri[(edge+1)%3]; - if(pi0>pi1) - Swap(pi0, pi1); + ngcore::ParallelForRange(openels.Size(), [&] (auto myrange) { + for (auto i_ : myrange) + { + auto i = openels[i_]; + const Element2d & tri = mesh.OpenElement(i); - // find other trig with edge pi0, pi1 - int i_other = -1; - for(auto ii : p2sel[pi0]) + for( auto edge : Range(3) ) { - if(ii==i) - continue; - auto & tri_other = mesh.OpenElement(ii); - if(tri_other[0]==pi1 || tri_other[1]==pi1 || tri_other[2]==pi1) + 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]) { - 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) + if(ii==i) + continue; + auto & tri_other = mesh.OpenElement(ii); + if(tri_other[0]==pi1 || tri_other[1]==pi1 || tri_other[2]==pi1) { - 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.Set(pi2); - badnode.Set(pi3); - } + 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.Test(el[0]) || - badnode.Test(el[1]) || - badnode.Test(el[2]) || - badnode.Test(el[3]) ) + if (badnode[el[0]] || + badnode[el[1]] || + badnode[el[2]] || + badnode[el[3]] ) tempels.DeleteElement(i); } } @@ -1643,13 +1652,10 @@ namespace netgen tempels.Append (el); } - NgBitArray badnode(mesh.GetNP()); - badnode.Clear(); - - DelaunayRemoveDegenerated(mesh.Points(), tempels, badnode, np); + DelaunayRemoveDegenerated(mesh.Points(), tempels, np); NgArray openels; - DelaunayRemoveTwoTriaTets(mesh, tempels, openels, badnode, np); + DelaunayRemoveTwoTriaTets(mesh, tempels, openels); DelaunayRemoveIntersecting(mesh, tempels, openels, pmin, pmax); DelaunayRemoveOuter(mesh, tempels, adfront);