Separate function CombineImproveEdge()

libsrc/meshing/improve2.cpp

@@ -354,138 +354,31 @@ namespace netgen
 
 
 
-
+  bool CombineImproveEdge( Mesh & mesh,
-
+                           const Table<SurfaceElementIndex, PointIndex> & elementsonnode,
-
+                           Array<Vec<3>, PointIndex> & normals,
- 
+                           Array<bool, PointIndex> & fixed,
-  void MeshOptimize2d :: CombineImprove (Mesh & mesh)
+                           PointIndex pi1, PointIndex pi2,
+                           bool check_only = true)
   {
-    if (!faceindex)
-      {
-        SplitImprove(mesh);
-	PrintMessage (3, "Combine improve");
-
-	for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
-	  {
-	    CombineImprove (mesh);
-
-	    if (multithread.terminate)
-	      throw NgException ("Meshing stopped");
-	  }
-	faceindex = 0;
-	return;
-      }
-
-
-    static Timer timer ("Combineimprove 2D");
-    RegionTimer reg (timer);
-
-    static Timer timerstart ("Combineimprove 2D start");
-    timerstart.Start();
-
-
-    static Timer timerstart1 ("Combineimprove 2D start1");
-    timerstart1.Start();
-
-    
-    Array<SurfaceElementIndex> seia;
-    mesh.GetSurfaceElementsOfFace (faceindex, seia);
-
-
-    for (SurfaceElementIndex sei : seia)
-      if (mesh[sei].GetNP() != 3)
-	return;
-
-
-    int surfnr = 0;
-    if (faceindex)
-      surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
-
-
     Vec<3> nv;
+    ArrayMem<SurfaceElementIndex, 20> hasonepi, hasbothpi;
 
-    int np = mesh.GetNP();
-
-    auto elementsonnode = mesh.CreatePoint2SurfaceElementTable(faceindex);
-    Array<SurfaceElementIndex> hasonepi, hasbothpi;
-
-    int ntasks = ngcore::TaskManager::GetMaxThreads();
-    Array<std::tuple<PointIndex, PointIndex>> edges;
-
-    BuildEdgeList( mesh, elementsonnode, edges );
-
-    Array<bool,PointIndex> fixed(np);
-    ParallelFor( fixed.Range(), [&fixed] (auto i) NETGEN_LAMBDA_INLINE
-            { fixed[i] = false; });
-
-    ParallelFor( edges.Range(), [&] (auto i) NETGEN_LAMBDA_INLINE
-            {
-              auto [pi0, pi1] = edges[i];
-              if (mesh.IsSegment (pi0, pi1))
-                {
-                  fixed[pi0] = true;
-                  fixed[pi1] = true;
-                }
-            });
-
-    timerstart1.Stop();
-
-    ParallelFor( mesh.LockedPoints().Range(), [&] (auto i) NETGEN_LAMBDA_INLINE
-            {
-              fixed[mesh.LockedPoints()[i]] = true;
-            });
-
-
-    Array<Vec<3>,PointIndex> normals(np);
-
-    ParallelFor( mesh.Points().Range(), [&] (auto pi) NETGEN_LAMBDA_INLINE
-        {
-            if (elementsonnode[pi].Size())
-              {
-                Element2d & hel = mesh[elementsonnode[pi][0]];
-                for (int k = 0; k < 3; k++)
-                  if (hel[k] == pi)
-                    {
-                      GetNormalVector (surfnr, mesh[pi], hel.GeomInfoPi(k+1), normals[pi]);
-                      break;
-                    }
-              }
-        });
-
-    timerstart.Stop();
-
-    for (auto e : edges)
-      {
-        auto [pi1, pi2] = e;
-
-	  {
-
-            /*
-	    if (pi1 < PointIndex::BASE || 
-		pi2 < PointIndex::BASE)
-	      continue;
-            */
     if (!pi1.IsValid() || !pi2.IsValid())
-	      continue;              
+        return false;
-	    /*
-	      INDEX_2 i2(pi1, pi2);
-	      i2.Sort();
-	      if (segmentht.Used(i2))
-	      continue;
-	    */
 
     bool debugflag = 0;
 
     if (debugflag)
       {
-		(*testout) << "Combineimprove, face = " << faceindex 
+        (*testout) << "Combineimprove "
             << "pi1 = " << pi1 << " pi2 = " << pi2 << endl;
       }
 
     /*
     // save version:
     if (fixed.Get(pi1) || fixed.Get(pi2))
-	    continue;
+    return false;
     if (pi2 < pi1) swap (pi1, pi2);
     */
 
@@ -494,23 +387,10 @@ namespace netgen
         Swap (pi1, pi2);
 
     if (fixed[pi2])
-	      continue;
+        return false;
 
     double loch = mesh.GetH (mesh[pi1]);
 
-	    // INDEX_2 si2 (pi1, pi2);
-	    // si2.Sort();
-
-	    /*	  
-	      if (edgetested.Used (si2))
-	      continue;
-	      edgetested.Set (si2, 1);
-	    */
-
-	    hasonepi.SetSize(0);
-	    hasbothpi.SetSize(0);
-
-	    // for (int k = 0; k < elementsonnode[pi1].Size(); k++)
     for (SurfaceElementIndex sei2 : elementsonnode[pi1])
       {
         const Element2d & el2 = mesh[sei2];
@@ -530,18 +410,10 @@ namespace netgen
       }
 
     if(hasbothpi.Size()==0)
-                continue;
+        return false;
 
 
-	    Element2d & hel = mesh[hasbothpi[0]];
+    nv = normals[pi1];
-	    for (int k = 0; k < 3; k++)
-	      if (hel[k] == pi1)
-		{
-		  GetNormalVector (surfnr, mesh[pi1], hel.GeomInfoPi(k+1), nv);
-		  break;
-		}
-
-	    //	  nv = normals.Get(pi1);
 
 
     for (SurfaceElementIndex sei2 :  elementsonnode[pi2])
@@ -600,8 +472,10 @@ namespace netgen
             bad2 += 1e10;
 
         for (int l = 0; l < 3; l++)
+          {
             if ( (normals[el[l]] * nv) < 0.5)
                 bad2 += 1e10;
+          }
 
         illegal2 += 1-mesh.LegalTrig(el);
       }
@@ -620,6 +494,9 @@ namespace netgen
       }
 
 
+    if(check_only)
+        return should;
+
     if (should)
       {
         /*
@@ -661,7 +538,6 @@ namespace netgen
           }
 
 
-
         // (*testout) << "Connect point " << pi2 << " to " << pi1 << "\n";
         // for (int k = 0; k < elementsonnode[pi2].Size(); k++)
         for (SurfaceElementIndex sei2 : elementsonnode[pi2])
@@ -684,8 +560,108 @@ namespace netgen
 
         for (auto sei : hasbothpi)
             mesh[sei].Delete();
+
+      }
+    return should;
+  }
+
+  void MeshOptimize2d :: CombineImprove (Mesh & mesh)
+  {
+    if (!faceindex)
+      {
+        SplitImprove(mesh);
+        PrintMessage (3, "Combine improve");
+
+        for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
+          {
+            CombineImprove (mesh);
+
+            if (multithread.terminate)
+              throw NgException ("Meshing stopped");
+          }
+        faceindex = 0;
+        return;
+      }
+
+
+    static Timer timer ("Combineimprove 2D");
+    RegionTimer reg (timer);
+
+    static Timer timerstart ("Combineimprove 2D start");
+    timerstart.Start();
+
+
+    static Timer timerstart1 ("Combineimprove 2D start1");
+    timerstart1.Start();
+
+    
+    Array<SurfaceElementIndex> seia;
+    mesh.GetSurfaceElementsOfFace (faceindex, seia);
+
+
+    for (SurfaceElementIndex sei : seia)
+      if (mesh[sei].GetNP() != 3)
+	return;
+
+
+    int surfnr = 0;
+    if (faceindex)
+      surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
+
+
+    int np = mesh.GetNP();
+
+    auto elementsonnode = mesh.CreatePoint2SurfaceElementTable(faceindex);
+
+    int ntasks = ngcore::TaskManager::GetMaxThreads();
+    Array<std::tuple<PointIndex, PointIndex>> edges;
+
+    BuildEdgeList( mesh, elementsonnode, edges );
+
+    Array<bool,PointIndex> fixed(np);
+    ParallelFor( fixed.Range(), [&fixed] (auto i) NETGEN_LAMBDA_INLINE
+            { fixed[i] = false; });
+
+    ParallelFor( edges.Range(), [&] (auto i) NETGEN_LAMBDA_INLINE
+            {
+              auto [pi0, pi1] = edges[i];
+              if (mesh.IsSegment (pi0, pi1))
+                {
+                  fixed[pi0] = true;
+                  fixed[pi1] = true;
+                }
+            });
+
+    timerstart1.Stop();
+
+    ParallelFor( mesh.LockedPoints().Range(), [&] (auto i) NETGEN_LAMBDA_INLINE
+            {
+              fixed[mesh.LockedPoints()[i]] = true;
+            });
+
+
+    Array<Vec<3>,PointIndex> normals(np);
+
+    ParallelFor( mesh.Points().Range(), [&] (auto pi) NETGEN_LAMBDA_INLINE
+        {
+            if (elementsonnode[pi].Size())
+              {
+                Element2d & hel = mesh[elementsonnode[pi][0]];
+                for (int k = 0; k < 3; k++)
+                  if (hel[k] == pi)
+                    {
+                      GetNormalVector (surfnr, mesh[pi], hel.GeomInfoPi(k+1), normals[pi]);
+                      break;
                     }
               }
+        });
+
+    timerstart.Stop();
+
+    for (auto e : edges)
+      {
+        auto [pi1, pi2] = e;
+        CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi1, pi2, false);
       }
 
     //  mesh.Compress();