From c075cfb9c6ad2ac20a9796e81061413141824ee9 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Wed, 4 Sep 2019 18:00:47 +0200 Subject: [PATCH] Parallelize ImproveMesh --- libsrc/meshing/meshclass.hpp | 1 + libsrc/meshing/smoothing3.cpp | 164 +++++++++++++++++++++++++++++++++- 2 files changed, 161 insertions(+), 4 deletions(-) diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index e8871e11..4b91fab3 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -542,6 +542,7 @@ namespace netgen DLL_HEADER void DoArchive (Archive & archive); /// DLL_HEADER void ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY); + DLL_HEADER void ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY); /// void ImproveMeshJacobian (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY, const NgBitArray * usepoint = NULL); diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index 1c8adec1..7d4854b8 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -5,10 +5,13 @@ #include #endif #include +#include +#include namespace netgen { + using namespace ngcore; double MinFunctionSum :: Func (const Vector & x) const @@ -301,7 +304,8 @@ namespace netgen public: Mesh::T_POINTS & points; const Array & elements; - TABLE elementsonpoint; + TABLE &elementsonpoint; + bool own_elementsonpoint; const MeshingParameters & mp; PointIndex actpind; double h; @@ -310,10 +314,12 @@ namespace netgen PointFunction (Mesh::T_POINTS & apoints, const Array & aelements, const MeshingParameters & amp); - virtual ~PointFunction () { ; } + PointFunction (const PointFunction & pf); + virtual ~PointFunction () { if(own_elementsonpoint) delete &elementsonpoint; } virtual void SetPointIndex (PointIndex aactpind); void SetLocalH (double ah) { h = ah; } double GetLocalH () const { return h; } + const TABLE & GetPointToElementTable() { return elementsonpoint; }; virtual double PointFunctionValue (const Point<3> & pp) const; virtual double PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const; virtual double PointFunctionValueDeriv (const Point<3> & pp, const Vec<3> & dir, double & deriv) const; @@ -322,11 +328,16 @@ namespace netgen }; + PointFunction :: PointFunction (const PointFunction & pf) + : points(pf.points), elements(pf.elements), elementsonpoint(pf.elementsonpoint), own_elementsonpoint(false), mp(pf.mp) + { } + PointFunction :: PointFunction (Mesh::T_POINTS & apoints, const Array & aelements, const MeshingParameters & amp) - : points(apoints), elements(aelements), elementsonpoint(apoints.Size()), mp(amp) + : points(apoints), elements(aelements), elementsonpoint(* new TABLE(apoints.Size())), own_elementsonpoint(true), mp(amp) { + static Timer tim("PointFunction - build elementsonpoint table"); RegionTimer reg(tim); for (int i = 0; i < elements.Size(); i++) if (elements[i].NP() == 4) for (int j = 0; j < elements[i].NP(); j++) @@ -1359,7 +1370,7 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal) -void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) +void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL goal) { static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t); @@ -1480,6 +1491,151 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) } } +void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) +{ + static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t); + static Timer tcoloring("coloring"); + static Timer tcalcbadmax("Calc badmax"); + static Timer topt("optimize"); + static Timer trange("range"); + + // return ImproveMeshSequential(mp, goal); + + (*testout) << "Improve Mesh" << "\n"; + PrintMessage (3, "ImproveMesh"); + + int np = GetNP(); + int ne = GetNE(); + + PointFunction pf_glob(points, volelements, mp); + + auto & elementsonpoint = pf_glob.GetPointToElementTable(); + + const auto & getDofs = [&] (int i) + { + i += PointIndex::BASE; + return FlatArray(elementsonpoint[i].Size(), &elementsonpoint[i][0]); + }; + + Array colors(points.Size()); + + tcoloring.Start(); + int ncolors = ngcore::ComputeColoring( colors, ne, getDofs ); + TableCreator creator(ncolors); + for ( ; !creator.Done(); creator++) + { + ParallelForRange( Range(colors), [&](auto myrange) + { + for(auto i : myrange) + creator.Add(colors[i], i); + }); + } + + auto color_table = creator.MoveTable(); + tcoloring.Stop(); + + if (goal == OPT_QUALITY) + { + double bad1 = CalcTotalBad (points, volelements, mp); + (*testout) << "Total badness = " << bad1 << endl; + PrintMessage (5, "Total badness = ", bad1); + } + + + (*testout) << setprecision(8); + + NgArray pointh (points.Size()); + + if(lochfunc) + { + for (PointIndex pi : points.Range()) + pointh[pi] = GetH(points[pi]); + } + else + { + pointh = 0; + for (Element & el : VolumeElements()) + { + double h = pow(el.Volume(points),1./3.); + for (PointIndex pi : el.PNums()) + if (h > pointh[pi]) + pointh[pi] = h; + } + } + + const char * savetask = multithread.task; + multithread.task = "Smooth Mesh"; + + topt.Start(); + int counter = 0; + for (int color : Range(color_table.Size())) + { + if (multithread.terminate) + throw NgException ("Meshing stopped"); + + ParallelForRange( Range(color_table[color].Size()), [&](auto myrange) + { + RegionTracer reg(ngcore::TaskManager::GetThreadId(), trange, myrange.Size()); + Vector x(3); + + PointFunction pf{pf_glob}; + + Opti3FreeMinFunction freeminf(pf); + + OptiParameters par; + par.maxit_linsearch = 20; + par.maxit_bfgs = 20; + + for (auto i : myrange) + { + PointIndex pi(color_table[color][i]+PointIndex::BASE); + if ( (*this)[pi].Type() == INNERPOINT ) + { + counter++; + + double lh = pointh[pi]; + pf.SetLocalH (lh); + par.typx = lh; + + freeminf.SetPoint (points[pi]); + pf.SetPointIndex (pi); + + x = 0; + int pok; + pok = freeminf.Func (x) < 1e10; + + if (!pok) + { + pok = pf.MovePointToInner (); + + freeminf.SetPoint (points[pi]); + pf.SetPointIndex (pi); + } + + if (pok) + { + //*testout << "start BFGS, pok" << endl; + BFGS (x, freeminf, par); + //*testout << "BFGS complete, pok" << endl; + points[pi](0) += x(0); + points[pi](1) += x(1); + points[pi](2) += x(2); + } + } + } + }, 4*ngcore::TaskManager::GetNumThreads()); + } + topt.Stop(); + + multithread.task = savetask; + + if (goal == OPT_QUALITY) + { + double bad1 = CalcTotalBad (points, volelements, mp); + (*testout) << "Total badness = " << bad1 << endl; + PrintMessage (5, "Total badness = ", bad1); + } +}