Parallel 2d MeshImprove

libsrc/meshing/smoothing2.cpp

@@ -709,35 +709,23 @@ 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++)
-	  {
-	    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;
+    auto elementsonpoint = mesh.CreatePoint2SurfaceElementTable( faceindex );
+    NgArray<MeshPoint, PointIndex::BASE> savepoints(mesh.GetNP());
 
+    Table<PointIndex> color_table;
+    if(faceindex)
+      {
         Array<SurfaceElementIndex> seia;
         mesh.GetSurfaceElementsOfFace (faceindex, seia);
-    bool mixed = 0;
         for (auto sei : seia)
           if (mesh[sei].GetNP() != 3)
             {
@@ -745,13 +733,7 @@ namespace netgen
               break;
             }
 
-    Vector x(2);
+        Array<int, PointIndex> compress(mesh.GetNP());
-
-    NgArray<MeshPoint, PointIndex::BASE> savepoints(mesh.GetNP());
-
-    ld.uselocalh = mp.uselocalh;
-
-    NgArray<int, PointIndex::BASE> compress(mesh.GetNP());
         NgArray<PointIndex> icompress;
         for (int i = 0; i < seia.Size(); i++)
           {
@@ -769,58 +751,52 @@ namespace netgen
                   icompress.Append(el[j]);
                 }
           }
-    NgArray<int> cnta(icompress.Size());
+
-    cnta = 0;
+        const auto & getDofs = [&] (int i)
-    for (int i = 0; i < seia.Size(); i++)
           {
-	const Element2d & el = mesh[seia[i]];
+            return elementsonpoint[icompress[i]];
-	for (int j = 0; j < el.GetNP(); j++)
+          };
-	  cnta[compress[el[j]]]++;
+
-      }
+        colors.SetSize(icompress.Size());
-    TABLE<SurfaceElementIndex> elementsonpoint(cnta);
+
-    for (int i = 0; i < seia.Size(); i++)
+        ncolors = ngcore::ComputeColoring( colors, mesh.GetNSE(), getDofs );
+
+        TableCreator<PointIndex> creator(ncolors);
+        for ( ; !creator.Done(); creator++)
+            ParallelForRange( Range(colors), [&](auto myrange)
                     {
-	const Element2d & el = mesh[seia[i]];
+                    for(auto i : myrange)
-	for (int j = 0; j < el.GetNP(); j++)
+                    creator.Add(colors[i], icompress[i]);
-	  elementsonpoint.Add (compress[el[j]], seia[i]);
+                    });
+
+        color_table = creator.MoveTable();
       }
-
+    else
-    
-    /*
-    NgArray<int, PointIndex::BASE> nelementsonpoint(mesh.GetNP());
-    nelementsonpoint = 0;
-    for (int i = 0; i < seia.Size(); i++)
       {
-	const Element2d & el = mesh[seia[i]];
+        for (auto & se : mesh.SurfaceElements())
-	for (int j = 0; j < el.GetNP(); j++)
+          if (se.GetNP() != 3)
-	  nelementsonpoint[el[j]]++;
-      }
-
-    TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonpoint(nelementsonpoint);
-
-    for (int i = 0; i < seia.Size(); i++)
             {
-	const Element2d & el = mesh[seia[i]];
+              mixed = true;
-	for (int j = 0; j < el.GetNP(); j++)
+              break;
-	  elementsonpoint.Add (el[j], seia[i]);
             }
-    */
+        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));
+                    });
 
-    ld.loch = mp.maxh;
+        color_table = creator.MoveTable();
-    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;
@@ -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,7 +955,9 @@ namespace netgen
 	  ld.t1 = ld.normal.GetNormal ();
 	  ld.t2 = Cross (ld.normal, ld.t1);
 	  
-	  // save points, and project to tangential plane
+          if(mixed)
+            {
+              // save points, and project to tangential plane (only for optimization with Opti2SurfaceMinFunctionJacobian in mixed element meshes)
               for (int j = 0; j < ld.locelements.Size(); j++)
                 {
                   const Element2d & el = mesh[ld.locelements[j]];
@@ -1007,20 +975,15 @@ namespace netgen
                       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,6 +992,8 @@ namespace netgen
 	  double fact = 1;
 	  int moveisok = 0;
 
+          if(mixed)
+            {
               // restore other points
               for (int j = 0; j < ld.locelements.Size(); j++)
                 {
@@ -1039,6 +1004,7 @@ namespace netgen
                       if (hhpi != pi) mesh[hhpi] = savepoints[hhpi];
                     }
                 }
+            }
 
 	  
 	  //optimizer loop (if whole distance is not possible, move only a bit!!!!)
@@ -1079,11 +1045,10 @@ namespace netgen
 	    }
           }
         }
+      }, mixed ? 1 : ngcore::TasksPerThread(4)); // mixed element smoothing not parallel yet
+    }
 
     tloop.Stop();
-    if (printeddot)
-      PrintDot ('\n');
-
     CheckMeshApproximation (mesh);
     mesh.SetNextTimeStamp();
   }