From a651a2d97e0fe9d5212d39cc3b13fb6f2d109bfd Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Tue, 8 Oct 2019 18:35:20 +0200 Subject: [PATCH] EdgeSwapping() - some cleanup and parallelization of table building --- libsrc/meshing/improve2.cpp | 148 +++++++++++++++++------------------- libsrc/meshing/improve2.hpp | 2 +- 2 files changed, 72 insertions(+), 78 deletions(-) diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp index 799c33c5..af0a6edf 100644 --- a/libsrc/meshing/improve2.cpp +++ b/libsrc/meshing/improve2.cpp @@ -8,6 +8,7 @@ namespace netgen { + using ngcore::ParallelForRange; @@ -28,7 +29,7 @@ namespace netgen Array &swapped, const SurfaceElementIndex t1, const int o1, const int t, - NgArray &pdef, + Array &pdef, const bool check_only) { bool should; @@ -175,6 +176,7 @@ namespace netgen void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric) { static Timer timer("EdgeSwapping (2D)"); RegionTimer reg(timer); + static Timer timer_nb("EdgeSwapping-Find neighbors"); if (usemetric) PrintMessage (3, "Edgeswapping, metric"); else @@ -215,59 +217,68 @@ namespace netgen Array swapped(mesh.GetNSE()); - NgArray pdef(mesh.GetNP()); - NgArray pangle(mesh.GetNP()); + Array pdef(mesh.GetNP()); + Array pangle(mesh.GetNP()); - - // int e; - // double d; - // Vec3d nv1, nv2; - - // double loch(-1); static const double minangle[] = { 0, 1.481, 2.565, 3.627, 4.683, 5.736, 7, 9 }; - for (int i = 0; i < seia.Size(); i++) + if(faceindex == 0) { - const Element2d & sel = mesh[seia[i]]; - for (int j = 0; j < 3; j++) - pangle[sel[j]] = 0.0; + ParallelForRange( Range(pangle), [&] (auto myrange) + { + for (auto i : myrange) + pangle[i] = 0.0; + }); } - // pangle = 0; - - for (int i = 0; i < seia.Size(); i++) + else { - const Element2d & sel = mesh[seia[i]]; - for (int j = 0; j < 3; j++) - { - POINTTYPE typ = mesh[sel[j]].Type(); - if (typ == FIXEDPOINT || typ == EDGEPOINT) - { - pangle[sel[j]] += - Angle (mesh[sel[(j+1)%3]] - mesh[sel[j]], - mesh[sel[(j+2)%3]] - mesh[sel[j]]); - } - } + ParallelForRange( Range(seia), [&] (auto myrange) + { + for (auto i : myrange) + { + const Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < 3; j++) + pangle[sel[j]] = 0.0; + } + }); } - // for (PointIndex pi = PointIndex::BASE; - // pi < mesh.GetNP()+PointIndex::BASE; pi++) - - // pdef = 0; - for (int i = 0; i < seia.Size(); i++) - { - const Element2d & sel = mesh[seia[i]]; - for (int j = 0; j < 3; j++) - { - PointIndex pi = sel[j]; - if (mesh[pi].Type() == INNERPOINT || mesh[pi].Type() == SURFACEPOINT) - pdef[pi] = -6; - else - for (int j = 0; j < 8; j++) - if (pangle[pi] >= minangle[j]) - pdef[pi] = -1-j; + ParallelForRange( Range(seia), [&] (auto myrange) + { + for (auto i : myrange) + { + const Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < 3; j++) + { + POINTTYPE typ = mesh[sel[j]].Type(); + if (typ == FIXEDPOINT || typ == EDGEPOINT) + { + pangle[sel[j]] += + Angle (mesh[sel[(j+1)%3]] - mesh[sel[j]], + mesh[sel[(j+2)%3]] - mesh[sel[j]]); + } + } } - } + }); + + ParallelForRange( Range(seia), [&] (auto myrange) + { + for (auto i : myrange) + { + const Element2d & sel = mesh[seia[i]]; + for (int j = 0; j < 3; j++) + { + PointIndex pi = sel[j]; + if (mesh[pi].Type() == INNERPOINT || mesh[pi].Type() == SURFACEPOINT) + pdef[pi] = -6; + else + for (int j = 0; j < 8; j++) + if (pangle[pi] >= minangle[j]) + pdef[pi] = -1-j; + } + } + }); /* for (int i = 0; i < seia.Size(); i++) @@ -277,41 +288,23 @@ namespace netgen pdef[sel[j]]++; } */ - for (SurfaceElementIndex sei : seia) - for (PointIndex pi : mesh[sei].PNums<3>()) - pdef[pi]++; - - // for (int i = 0; i < seia.Size(); i++) - for (SurfaceElementIndex sei : seia) - for (int j = 0; j < 3; j++) + ParallelForRange( Range(seia), [&] (auto myrange) { - neighbors[sei].SetNr (j, -1); - neighbors[sei].SetOrientation (j, 0); - } - - /* - NgArray normals(mesh.GetNP()); - for (i = 1; i <= mesh.GetNSE(); i++) - { - Element2d & hel = mesh.SurfaceElement(i); - if (hel.GetIndex() == faceindex) - for (k = 1; k <= 3; k++) - { - int pi = hel.PNum(k); - SelectSurfaceOfPoint (mesh.Point(pi), hel.GeomInfoPi(k)); - int surfi = mesh.GetFaceDescriptor(faceindex).SurfNr(); - GetNormalVector (surfi, mesh.Point(pi), normals.Elem(pi)); - normals.Elem(pi) /= normals.Elem(pi).Length(); - } - } - */ - - /* - for (int i = 0; i < seia.Size(); i++) - { - const Element2d & sel = mesh[seia[i]]; - */ + for (auto i : myrange) + { + auto sei = seia[i]; + for (PointIndex pi : mesh[sei].template PNums<3>()) + AsAtomic(pdef[pi])++; + for (int j = 0; j < 3; j++) + { + neighbors[sei].SetNr (j, -1); + neighbors[sei].SetOrientation (j, 0); + } + } + }); + + timer_nb.Start(); for (SurfaceElementIndex sei : seia) { const Element2d & sel = mesh[sei]; @@ -358,6 +351,7 @@ namespace netgen } } } + timer_nb.Stop(); for (SurfaceElementIndex sei : seia) swapped[sei] = false; diff --git a/libsrc/meshing/improve2.hpp b/libsrc/meshing/improve2.hpp index 3d62ee45..03ff5559 100644 --- a/libsrc/meshing/improve2.hpp +++ b/libsrc/meshing/improve2.hpp @@ -37,7 +37,7 @@ public: const NgArray* > & from, NgArray* > & dest); bool EdgeSwapping (Mesh & mesh, const int usemetric, Array &neighbors, Array &swapped, - const SurfaceElementIndex t1, const int edge, const int t, NgArray &pdef, const bool check_only=false); + const SurfaceElementIndex t1, const int edge, const int t, Array &pdef, const bool check_only=false); void EdgeSwapping (Mesh & mesh, int usemetric); void CombineImprove (Mesh & mesh); void SplitImprove (Mesh & mesh);