From d0edaa57bb8aa62215c19893e9702b8a7ac56eae Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Wed, 16 Jun 2021 15:05:58 +0200 Subject: [PATCH] prepare LocalH tree before blockfill sequentially --- libsrc/meshing/meshfunc.cpp | 68 +++++++++++------- libsrc/meshing/meshing3.cpp | 134 +++++++++++++++++++++++------------- libsrc/meshing/meshing3.hpp | 2 + 3 files changed, 133 insertions(+), 71 deletions(-) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 8778daa3..5eb8801c 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -227,6 +227,26 @@ namespace netgen } } + void PrepareForBlockFillLocalH(MeshingData & md) + { + static Timer t("PrepareForBlockFillLocalH"); RegionTimer rt(t); + md.meshing = make_unique(nullptr); + + auto & mesh = md.mesh; + + mesh.CalcSurfacesOfNode(); + mesh.FindOpenElements(md.domain); + + for (PointIndex pi : mesh.Points().Range()) + md.meshing->AddPoint (mesh[pi], pi); + + for (int i = 1; i <= mesh.GetNOpenElements(); i++) + md.meshing->AddBoundaryElement (mesh.OpenElement(i)); + + if (mesh.HasLocalHFunction()) + md.meshing->PrepareBlockFillLocalH(mesh, md.mp); + } + void MeshDomain( MeshingData & md) { @@ -237,15 +257,7 @@ namespace netgen if(mp.only3D_domain_nr && mp.only3D_domain_nr != domain) return; - mesh.CalcSurfacesOfNode(); - mesh.FindOpenElements(domain); - md.meshing = make_unique(nullptr); - for (PointIndex pi : mesh.Points().Range()) - md.meshing->AddPoint (mesh[pi], pi); - - for (int i = 1; i <= mesh.GetNOpenElements(); i++) - md.meshing->AddBoundaryElement (mesh.OpenElement(i)); if (mp.delaunay && mesh.GetNOpenElements()) @@ -261,9 +273,18 @@ namespace netgen mesh.GetNE(), " elements"); } - // mesh.Save("before_findopenels.vol"); - // mesh.CalcSurfacesOfNode(); - // mesh.FindOpenElements(domain); + Box<3> domain_bbox( Box<3>::EMPTY_BOX ); + + for (auto & sel : mesh.SurfaceElements()) + { + if (sel.IsDeleted() ) continue; + + for (auto pi : sel.PNums()) + domain_bbox.Add (mesh[pi]); + } + domain_bbox.Increase (0.01 * domain_bbox.Diam()); + + mesh.FindOpenElements(domain); int cntsteps = 0; int meshed; @@ -273,7 +294,7 @@ namespace netgen if (multithread.terminate) break; - mesh.FindOpenElements(); + mesh.FindOpenElements(domain); PrintMessage (5, mesh.GetNOpenElements(), " open faces"); GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol"); cntsteps++; @@ -289,15 +310,14 @@ namespace netgen Array glob2loc(mesh.GetNP()); glob2loc = PointIndex::INVALID; + for (PointIndex pi : mesh.Points().Range()) + if (domain_bbox.IsIn (mesh[pi])) + glob2loc[pi] = meshing.AddPoint (mesh[pi], pi); + for (auto sel : mesh.OpenElements() ) { for(auto & pi : sel.PNums()) - { - if(!glob2loc[pi].IsValid()) - glob2loc[pi] = meshing.AddPoint (mesh[pi], pi); - pi = glob2loc[pi]; - } meshing.AddBoundaryElement (sel); } @@ -312,7 +332,7 @@ namespace netgen mesh.CalcSurfacesOfNode(); - mesh.FindOpenElements(); + mesh.FindOpenElements(domain); // teterrpow = 2; if (mesh.GetNOpenElements() != 0) @@ -358,7 +378,7 @@ namespace netgen { PrintMessage (3, "Check subdomain ", domain, " / ", mesh.GetNDomains()); - mesh.FindOpenElements(); + mesh.FindOpenElements(domain); bool res = (mesh.CheckConsistentBoundary() != 0); if (res) @@ -710,14 +730,12 @@ namespace netgen CloseOpenQuads( md[i] ); }); + for(auto & md_ : md) + PrepareForBlockFillLocalH(md_); + ParallelFor( md.Range(), [&](int i) { - // try { - MeshDomain(md[i]); - // } - // catch( const NgException & e ) { - // md[i].mesh.Save("error_"+ToString(i)+".vol"); - // } + MeshDomain(md[i]); }); MergeMeshes(mesh3d, md); diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp index c7424b51..8c4ee5aa 100644 --- a/libsrc/meshing/meshing3.cpp +++ b/libsrc/meshing/meshing3.cpp @@ -180,6 +180,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) // static int meshing3_timer_d = NgProfiler::CreateTimer ("Meshing3::GenerateMesh d"); // NgProfiler::RegionTimer reg (meshing3_timer); + cout << "start tet meshing with " << adfront->GetNP() << " points and " << adfront->GetNF() << " faces " << endl; NgArray locpoints; // local points NgArray locfaces; // local faces @@ -1094,11 +1095,10 @@ static int TestSameSide (const Point3d & p1, const Point3d & p2) */ - -void Meshing3 :: BlockFillLocalH (Mesh & mesh, +void Meshing3 :: PrepareBlockFillLocalH (Mesh & mesh, const MeshingParameters & mp) { - static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t); + static Timer t("Mesing3::PrepareBlockFillLocalH"); RegionTimer reg(t); double filldist = mp.filldist; @@ -1107,10 +1107,90 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, PrintMessage (3, "blockfill local h"); - NgArray > npoints; - adfront -> CreateTrees(); + double maxh = 0; + + for (int i = 1; i <= adfront->GetNF(); i++) + { + const MiniElement2d & el = adfront->GetFace(i); + for (int j = 1; j <= 3; j++) + { + const Point3d & p1 = adfront->GetPoint (el.PNumMod(j)); + const Point3d & p2 = adfront->GetPoint (el.PNumMod(j+1)); + + double hi = Dist (p1, p2); + if (hi > maxh) maxh = hi; + } + } + + + if (mp.maxh < maxh) maxh = mp.maxh; + + // auto loch_ptr = mesh.LocalHFunction().Copy(); + // auto & loch = *loch_ptr; + auto & loch = mesh.LocalHFunction(); + + bool changed; + static Timer t1("loop1"); + t1.Start(); + do + { + loch.ClearFlags(); + + static Timer tbox("adfront-bbox"); + tbox.Start(); + for (int i = 1; i <= adfront->GetNF(); i++) + { + const MiniElement2d & el = adfront->GetFace(i); + + Box<3> bbox (adfront->GetPoint (el[0])); + bbox.Add (adfront->GetPoint (el[1])); + bbox.Add (adfront->GetPoint (el[2])); + + + double filld = filldist * bbox.Diam(); + bbox.Increase (filld); + + loch.CutBoundary (bbox); // .PMin(), bbox.PMax()); + } + tbox.Stop(); + + // locadfront = adfront; + loch.FindInnerBoxes (adfront, NULL); + + npoints.SetSize(0); + loch.GetInnerPoints (npoints); + + changed = false; + for (int i = 1; i <= npoints.Size(); i++) + { + if (loch.GetH(npoints.Get(i)) > 1.5 * maxh) + { + loch.SetH (npoints.Get(i), maxh); + changed = true; + } + } + } + while (changed); + t1.Stop(); + + + +} + +void Meshing3 :: BlockFillLocalH (Mesh & mesh, + const MeshingParameters & mp) +{ + static Timer t("Mesing3::BlockFillLocalH"); RegionTimer reg(t); + // PrepareBlockFillLocalH(mesh, mp); + + double filldist = mp.filldist; + + // (*testout) << "blockfill local h" << endl; + // (*testout) << "rel filldist = " << filldist << endl; + PrintMessage (3, "blockfill local h"); + Box<3> bbox ( Box<3>::EMPTY_BOX ); double maxh = 0; @@ -1144,47 +1224,6 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, if (mp.maxh < maxh) maxh = mp.maxh; - auto loch_ptr = mesh.LocalHFunction().Copy(); - auto & loch = *loch_ptr; - - bool changed; - do - { - loch.ClearFlags(); - - for (int i = 1; i <= adfront->GetNF(); i++) - { - const MiniElement2d & el = adfront->GetFace(i); - - Box<3> bbox (adfront->GetPoint (el[0])); - bbox.Add (adfront->GetPoint (el[1])); - bbox.Add (adfront->GetPoint (el[2])); - - - double filld = filldist * bbox.Diam(); - bbox.Increase (filld); - - loch.CutBoundary (bbox); // .PMin(), bbox.PMax()); - } - - // locadfront = adfront; - loch.FindInnerBoxes (adfront, NULL); - - npoints.SetSize(0); - loch.GetInnerPoints (npoints); - - changed = false; - for (int i = 1; i <= npoints.Size(); i++) - { - if (loch.GetH(npoints.Get(i)) > 1.5 * maxh) - { - loch.SetH (npoints.Get(i), maxh); - changed = true; - } - } - } - while (changed); - if (debugparam.slowchecks) (*testout) << "Blockfill with points: " << endl; for (int i = 1; i <= npoints.Size(); i++) @@ -1211,6 +1250,8 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, // find outer points + static Timer tloch2("build loch2"); + tloch2.Start(); loch2.ClearFlags(); for (int i = 1; i <= adfront->GetNF(); i++) @@ -1248,6 +1289,7 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh, // loch2.CutBoundary (pmin, pmax); loch2.CutBoundary (Box<3> (pmin, pmax)); // pmin, pmax); } + tloch2.Stop(); // locadfront = adfront; loch2.FindInnerBoxes (adfront, NULL); diff --git a/libsrc/meshing/meshing3.hpp b/libsrc/meshing/meshing3.hpp index 65b153b9..2944c6c3 100644 --- a/libsrc/meshing/meshing3.hpp +++ b/libsrc/meshing/meshing3.hpp @@ -28,6 +28,7 @@ class Meshing3 NgArray problems; /// tolerance criterion double tolfak; + NgArray > npoints; public: /// Meshing3 (const string & rulefilename); @@ -63,6 +64,7 @@ public: /// void BlockFill (Mesh & mesh, double gh); /// + void PrepareBlockFillLocalH (Mesh & mesh, const MeshingParameters & mp); void BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp); /// uses points of adfront, and puts new elements into mesh