Merge branch 'parallel_improve_mesh_2d' into 'master'

Parallel improve mesh 2d

See merge request jschoeberl/netgen!270

libsrc/core/taskmanager.hpp

@@ -15,6 +15,7 @@
 
 #include "array.hpp"
 #include "paje_trace.hpp"
+#include "profiler.hpp"
 
 namespace ngcore
 {
@@ -1059,6 +1060,7 @@ public:
   template <typename Tmask>
   int ComputeColoring( FlatArray<int> colors, size_t ndofs, Tmask const & getDofs)
   {
+    static Timer timer("ComputeColoring - "+Demangle(typeid(Tmask).name())); RegionTimer rt(timer);
     static_assert(sizeof(unsigned int)==4, "Adapt type of mask array");
     auto n = colors.Size();
 

libsrc/meshing/improve2.cpp

@@ -177,8 +177,8 @@ namespace netgen
     else
       PrintMessage (3, "Edgeswapping, topological");
 
-    static int timerstart = NgProfiler::CreateTimer ("EdgeSwapping 2D start");
+    static Timer timerstart("EdgeSwapping 2D start");
-    NgProfiler::StartTimer (timerstart);
+    timerstart.Start();
 
     Array<SurfaceElementIndex> seia;
     bool mixed = false;
@@ -310,7 +310,7 @@ namespace netgen
     for (SurfaceElementIndex sei : seia)
       swapped[sei] = false;
     
-    NgProfiler::StopTimer (timerstart);
+    timerstart.Stop();
   
 
 
@@ -377,15 +377,15 @@ namespace netgen
       }
 
 
-    static int timer = NgProfiler::CreateTimer ("Combineimprove 2D");
+    static Timer timer ("Combineimprove 2D");
-    NgProfiler::RegionTimer reg (timer);
+    RegionTimer reg (timer);
 
-    static int timerstart = NgProfiler::CreateTimer ("Combineimprove 2D start");
+    static Timer timerstart ("Combineimprove 2D start");
-    NgProfiler::StartTimer  (timerstart);
+    timerstart.Start();
 
 
-    static int timerstart1 = NgProfiler::CreateTimer ("Combineimprove 2D start1");
+    static Timer timerstart1 ("Combineimprove 2D start1");
-    NgProfiler::StartTimer  (timerstart1);
+    timerstart1.Start();
 
     
     Array<SurfaceElementIndex> seia;
@@ -416,7 +416,7 @@ namespace netgen
     Array<bool,PointIndex> fixed(np);
     fixed = false;
 
-    NgProfiler::StopTimer  (timerstart1);
+    timerstart1.Stop();
 
     /*
     for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
@@ -468,7 +468,7 @@ namespace netgen
 	  }
       }
 
-    NgProfiler::StopTimer  (timerstart);
+    timerstart.Stop();
 
     for (int i = 0; i < seia.Size(); i++)
       {

libsrc/meshing/smoothing2.cpp

@@ -709,119 +709,95 @@ namespace netgen
 
   void MeshOptimize2d :: ImproveMesh (Mesh & mesh, const MeshingParameters & mp)
   {
-    if (!faceindex)
+    static Timer timer("MeshSmoothing 2D"); RegionTimer reg (timer);
-      {
-	PrintMessage (3, "Smoothing");
 
-	for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
+    PrintMessage (3, "Smoothing");
-	  {
-	    ImproveMesh (mesh, mp);
-	    if (multithread.terminate)
-	      throw NgException ("Meshing stopped");
-	  }
-	faceindex = 0;
-	return;
-      }
-
-    static Timer timer("MeshSmoothing 2D");
-    // static int timer1 = NgProfiler::CreateTimer ("MeshSmoothing 2D start");
-    // static int timer2 = NgProfiler::CreateTimer ("MeshSmoothing 2D - BFGS");
-
-    RegionTimer reg (timer);
-    // NgProfiler::StartTimer (timer1);
 
     CheckMeshApproximation (mesh);
 
-    Opti2dLocalData ld;
+    int ncolors;
-
+    Array<int> colors;
-
+    bool mixed = false;
-    Array<SurfaceElementIndex> seia;
+    auto elementsonpoint = mesh.CreatePoint2SurfaceElementTable( faceindex );
-    mesh.GetSurfaceElementsOfFace (faceindex, seia);
-    bool mixed = 0;
-    for (auto sei : seia)
-      if (mesh[sei].GetNP() != 3)
-	{
-	  mixed = true;
-	  break;
-	}
-
-    Vector x(2);
-
     NgArray<MeshPoint, PointIndex::BASE> savepoints(mesh.GetNP());
 
-    ld.uselocalh = mp.uselocalh;
+    Table<PointIndex> color_table;
+    if(faceindex)
+      {
+        Array<SurfaceElementIndex> seia;
+        mesh.GetSurfaceElementsOfFace (faceindex, seia);
+        for (auto sei : seia)
+          if (mesh[sei].GetNP() != 3)
+            {
+              mixed = true;
+              break;
+            }
 
-    NgArray<int, PointIndex::BASE> compress(mesh.GetNP());
+        Array<int, PointIndex> compress(mesh.GetNP());
-    NgArray<PointIndex> icompress;
+        NgArray<PointIndex> icompress;
-    for (int i = 0; i < seia.Size(); i++)
+        for (int i = 0; i < seia.Size(); i++)
-      {
+          {
-	const Element2d & el = mesh[seia[i]];
+            const Element2d & el = mesh[seia[i]];
-	for (int j = 0; j < el.GetNP(); j++)
+            for (int j = 0; j < el.GetNP(); j++)
-	  compress[el[j]] = -1;
+              compress[el[j]] = -1;
-      }
+          }
-    for (int i = 0; i < seia.Size(); i++)
+        for (int i = 0; i < seia.Size(); i++)
-      {
+          {
-	const Element2d & el = mesh[seia[i]];
+            const Element2d & el = mesh[seia[i]];
-	for (int j = 0; j < el.GetNP(); j++)
+            for (int j = 0; j < el.GetNP(); j++)
-	  if (compress[el[j]] == -1)
+              if (compress[el[j]] == -1)
-	    {
+                {
-	      compress[el[j]] = icompress.Size();
+                  compress[el[j]] = icompress.Size();
-	      icompress.Append(el[j]);
+                  icompress.Append(el[j]);
-	    }
+                }
-      }
+          }
-    NgArray<int> cnta(icompress.Size());
-    cnta = 0;
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & el = mesh[seia[i]];
-	for (int j = 0; j < el.GetNP(); j++)
-	  cnta[compress[el[j]]]++;
-      }
-    TABLE<SurfaceElementIndex> elementsonpoint(cnta);
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & el = mesh[seia[i]];
-	for (int j = 0; j < el.GetNP(); j++)
-	  elementsonpoint.Add (compress[el[j]], seia[i]);
-      }
 
+        const auto & getDofs = [&] (int i)
+          {
+            return elementsonpoint[icompress[i]];
+          };
+
+        colors.SetSize(icompress.Size());
+
+        ncolors = ngcore::ComputeColoring( colors, mesh.GetNSE(), getDofs );
+
+        TableCreator<PointIndex> creator(ncolors);
+        for ( ; !creator.Done(); creator++)
+            ParallelForRange( Range(colors), [&](auto myrange)
+                    {
+                    for(auto i : myrange)
+                    creator.Add(colors[i], icompress[i]);
+                    });
+
+        color_table = creator.MoveTable();
+      }
+    else
+      {
+        for (auto & se : mesh.SurfaceElements())
+          if (se.GetNP() != 3)
+            {
+              mixed = true;
+              break;
+            }
+        const auto & getDofs = [&] (int i)
+          {
+            return elementsonpoint[i+PointIndex::BASE];
+          };
+
+        colors.SetSize(mesh.GetNP());
+        ncolors = ngcore::ComputeColoring( colors, mesh.GetNSE(), getDofs );
+
+        TableCreator<PointIndex> creator(ncolors);
+        for ( ; !creator.Done(); creator++)
+            ParallelForRange( Range(colors), [&](auto myrange)
+                    {
+                    for(auto i : myrange)
+                    creator.Add(colors[i], PointIndex(i+PointIndex::BASE));
+                    });
+
+        color_table = creator.MoveTable();
+      }
     
-    /*
-    NgArray<int, PointIndex::BASE> nelementsonpoint(mesh.GetNP());
-    nelementsonpoint = 0;
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & el = mesh[seia[i]];
-	for (int j = 0; j < el.GetNP(); j++)
-	  nelementsonpoint[el[j]]++;
-      }
-
-    TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonpoint(nelementsonpoint);
-
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & el = mesh[seia[i]];
-	for (int j = 0; j < el.GetNP(); j++)
-	  elementsonpoint.Add (el[j], seia[i]);
-      }
-    */
-    
-
-
-    ld.loch = mp.maxh;
-    ld.locmetricweight = metricweight;
-    ld.meshthis = this;
-
-
-
-    Opti2SurfaceMinFunction surfminf(mesh, ld);
-    Opti2EdgeMinFunction edgeminf(mesh, ld);
-    Opti2SurfaceMinFunctionJacobian surfminfj(mesh, ld);
-
-    OptiParameters par;
-    par.maxit_linsearch = 8;
-    par.maxit_bfgs = 5;
-
     /*
     int i, j, k;
     Vector xedge(1);
@@ -891,27 +867,6 @@ namespace netgen
     */
 
 
-    bool printeddot = 0;
-    char plotchar = '.';
-    int modplot = 1;
-    if (mesh.GetNP() > 1000)
-      {
-	plotchar = '+';
-	modplot = 100;
-      }
-    if (mesh.GetNP() > 10000)
-      {
-	plotchar = 'o';
-	modplot = 1000;
-      }
-    if (mesh.GetNP() > 100000)
-      {
-	plotchar = 'O';
-	modplot = 10000;
-      }
-    int cnt = 0;
-
-
     // NgProfiler::StopTimer (timer1);
 
     /*
@@ -921,28 +876,39 @@ namespace netgen
 
     static Timer tloop("MeshSmooting 2D - loop");
     tloop.Start();
-    for (int hi = 0; hi < icompress.Size(); hi++)
+    for (auto icolor : Range(color_table))
-      {
+    {
-	PointIndex pi = icompress[hi];
+     if (multithread.terminate)
+         break;
+     ParallelForRange( Range(color_table[icolor].Size()), [&](auto myrange)
+     {
+      Opti2dLocalData ld;
+      ld.uselocalh = mp.uselocalh;
+      ld.loch = mp.maxh;
+      ld.locmetricweight = metricweight;
+      ld.meshthis = this;
+
+      Opti2SurfaceMinFunction surfminf(mesh, ld);
+      Opti2SurfaceMinFunctionJacobian surfminfj(mesh, ld);
+
+      MinFunction & minfunc = mixed ? static_cast<MinFunction&>(surfminfj) : surfminf;
+
+      OptiParameters par;
+      par.maxit_linsearch = 8;
+      par.maxit_bfgs = 5;
+      for (auto i : myrange)
+        {
+	PointIndex pi = color_table[icolor][i];
 	if (mesh[pi].Type() == SURFACEPOINT)
 	  {
 	    if (multithread.terminate)
-	      throw NgException ("Meshing stopped");
+	      return;
 	    
-	    cnt++;
+	    if (elementsonpoint[pi].Size() == 0) continue;
-	    if (cnt % modplot == 0 && writestatus)
-	      {
-		printeddot = 1;
-		PrintDot (plotchar);
-	      }
-	    
-	    // if (elementsonpoint[pi].Size() == 0) continue;
-	    if (elementsonpoint[hi].Size() == 0) continue;
             
 	    ld.sp1 = mesh[pi];
 	    
-	    // Element2d & hel = mesh[elementsonpoint[pi][0]];
+	    Element2d & hel = mesh[elementsonpoint[pi][0]];
-	    Element2d & hel = mesh[elementsonpoint[hi][0]];
 	    
 	    int hpi = 0;
 	    for (int j = 1; j <= hel.GetNP(); j++)
@@ -960,9 +926,9 @@ namespace netgen
             ld.loc_pnts2.SetSize (0);
             ld.loc_pnts3.SetSize (0);
 
-	    for (int j = 0; j < elementsonpoint[hi].Size(); j++)
+	    for (int j = 0; j < elementsonpoint[pi].Size(); j++)
 	      {
-		SurfaceElementIndex sei = elementsonpoint[hi][j];
+		SurfaceElementIndex sei = elementsonpoint[pi][j];
 		const Element2d & bel = mesh[sei];
 		ld.surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();
 		
@@ -989,38 +955,35 @@ namespace netgen
 	  ld.t1 = ld.normal.GetNormal ();
 	  ld.t2 = Cross (ld.normal, ld.t1);
 	  
-	  // save points, and project to tangential plane
+          if(mixed)
-	  for (int j = 0; j < ld.locelements.Size(); j++)
+            {
-	    {
+              // save points, and project to tangential plane (only for optimization with Opti2SurfaceMinFunctionJacobian in mixed element meshes)
-	      const Element2d & el = mesh[ld.locelements[j]];
+              for (int j = 0; j < ld.locelements.Size(); j++)
-	      for (int k = 0; k < el.GetNP(); k++)
+                {
-		savepoints[el[k]] = mesh[el[k]];
+                  const Element2d & el = mesh[ld.locelements[j]];
-	    }
+                  for (int k = 0; k < el.GetNP(); k++)
+                    savepoints[el[k]] = mesh[el[k]];
+                }
 
-	  for (int j = 0; j < ld.locelements.Size(); j++)
+              for (int j = 0; j < ld.locelements.Size(); j++)
-	    {
+                {
-	      const Element2d & el = mesh[ld.locelements[j]];
+                  const Element2d & el = mesh[ld.locelements[j]];
-	      for (int k = 0; k < el.GetNP(); k++)
+                  for (int k = 0; k < el.GetNP(); k++)
-		{
+                    {
-		  PointIndex hhpi = el[k];
+                      PointIndex hhpi = el[k];
-		  double lam = ld.normal * (mesh[hhpi] - ld.sp1);
+                      double lam = ld.normal * (mesh[hhpi] - ld.sp1);
-		  mesh[hhpi] -= lam * ld.normal;
+                      mesh[hhpi] -= lam * ld.normal;
-		}
+                    }
-	    }
+                }
+            }
 	  
+          Vector x(2);
 	  x = 0;
 	  par.typx = 0.3*ld.lochs[0];
 
           // NgProfiler::StartTimer (timer2);
 
-	  if (mixed)
+          BFGS (x, minfunc, par, 1e-6);
-	    {
-	      BFGS (x, surfminfj, par, 1e-6);
-	    }
-	  else
-	    {
-	      BFGS (x, surfminf, par, 1e-6);
-	    }
 
           // NgProfiler::StopTimer (timer2);
 
@@ -1029,16 +992,19 @@ namespace netgen
 	  double fact = 1;
 	  int moveisok = 0;
 
-	  // restore other points
+          if(mixed)
-	  for (int j = 0; j < ld.locelements.Size(); j++)
+            {
-	    {
+              // restore other points
-	      const Element2d & el = mesh[ld.locelements[j]];
+              for (int j = 0; j < ld.locelements.Size(); j++)
-	      for (int k = 0; k < el.GetNP(); k++)
+                {
-		{
+                  const Element2d & el = mesh[ld.locelements[j]];
-		  PointIndex hhpi = el[k];
+                  for (int k = 0; k < el.GetNP(); k++)
-		  if (hhpi != pi) mesh[hhpi] = savepoints[hhpi];
+                    {
-		}
+                      PointIndex hhpi = el[k];
-	    }
+                      if (hhpi != pi) mesh[hhpi] = savepoints[hhpi];
+                    }
+                }
+            }
 
 	  
 	  //optimizer loop (if whole distance is not possible, move only a bit!!!!)
@@ -1077,13 +1043,12 @@ namespace netgen
 		}
 	    
 	    }
-	}
+          }
-      }
+        }
+      }, mixed ? 1 : ngcore::TasksPerThread(4)); // mixed element smoothing not parallel yet
+    }
 
     tloop.Stop();
-    if (printeddot)
-      PrintDot ('\n');
-
     CheckMeshApproximation (mesh);
     mesh.SetNextTimeStamp();
   }