fix parallel surface optimization with occ

libsrc/meshing/basegeom.cpp

@@ -982,11 +982,10 @@ namespace netgen
     static Timer timer_opt2d("Optimization 2D");
     RegionTimer reg(timer_opt2d);
     auto meshopt = MeshOptimize2d(mesh);
+    meshopt.SetFaceIndex(0);
     for(auto i : Range(mparam.optsteps2d))
-    for(auto k : Range(mesh.GetNFD()))
       {
         PrintMessage(3, "Optimization step ", i);
-        meshopt.SetFaceIndex(k+1);
         int innerstep = 0;
         for(auto optstep : mparam.optimize2d)
           {

libsrc/meshing/improve2.cpp

@@ -18,6 +18,27 @@ namespace netgen
     { tnr = atnr; sidenr = asidenr; }
   };
 
+  // check if element is quad with at least one surface point -> relevant for optimization
+  // (quads with 4 edge points are not optimized and can be ignored)
+  bool checkMixedElement(const Mesh & mesh, FlatArray<SurfaceElementIndex> seia)
+  {
+      bool mixed = false;
+      ParallelForRange( Range(seia), [&] (auto myrange) NETGEN_LAMBDA_INLINE
+         {
+            for (SurfaceElementIndex i : myrange)
+            {
+              const auto & sel = mesh[i];
+
+              if(sel.GetNP() == 3)
+                  continue;
+
+              for(auto i : Range(sel.GetNP()))
+                  if(mesh[sel[i]].Type() == SURFACEPOINT)
+                      mixed = true;
+            }
+         });
+      return mixed;
+  }
 
   bool MeshOptimize2d :: EdgeSwapping (const int usemetric,
     Array<Neighbour> &neighbors,
@@ -181,34 +202,13 @@ namespace netgen
     timerstart.Start();
 
     Array<SurfaceElementIndex> seia;
-    bool mixed = false;
-
-    if(faceindex==0)
-      {
-        seia.SetSize(mesh.GetNSE());
-        ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
-            {
-              SurfaceElementIndex sei(i);
-              seia[i] = sei;
-              if (mesh[sei].GetNP() != 3)
-              {
-                  const auto & sel = mesh[sei];
-                  for(auto i : Range(sel.GetNP()))
-                    if(mesh[sel[i]].Type() == INNERPOINT)
-                      mixed = true;
-              }
-            });
-      }
-    else
-      {
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
-        for (SurfaceElementIndex sei : seia)
-            if (mesh[sei].GetNP() != 3)
-                mixed = true;
-      }
 
-    if(mixed)
+    if(checkMixedElement(mesh, seia))
+    {
+        timerstart.Stop();
         return GenericImprove();
+    }
 
     Array<Neighbour> neighbors(mesh.GetNSE());
     auto elements_on_node = mesh.CreatePoint2SurfaceElementTable(faceindex);
@@ -595,25 +595,14 @@ namespace netgen
 
     
     Array<SurfaceElementIndex> seia;
-
-    if(faceindex)
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
-    else
-      {
-        seia.SetSize(mesh.GetNSE());
-        ParallelFor( IntRange(mesh.GetNSE()), [&seia] (auto i) NETGEN_LAMBDA_INLINE
-                { seia[i] = i; });
-      }
 
-    bool mixed = false;
+    if(checkMixedElement(mesh, seia))
-    ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
     {
-                if (mesh[seia[i]].GetNP() != 3)
+        timerstart1.Stop();
-                    mixed = true;
+        timerstart.Stop();
-            });
-
-    if(mixed)
         return;
+    }
 
     int np = mesh.GetNP();
 
@@ -625,18 +614,8 @@ namespace netgen
     BuildEdgeList( mesh, elementsonnode, edges );
 
     Array<bool,PointIndex> fixed(np);
-    ParallelFor( fixed.Range(), [&fixed] (auto i) NETGEN_LAMBDA_INLINE
+    ParallelFor( fixed.Range(), [&] (auto i) NETGEN_LAMBDA_INLINE
-            { fixed[i] = false; });
+            { fixed[i] = mesh[i].Type() != SURFACEPOINT; });
-
-    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();
 

libsrc/meshing/meshclass.cpp

@@ -6412,25 +6412,17 @@ namespace netgen
     static int timer = NgProfiler::CreateTimer ("GetSurfaceElementsOfFace");
     NgProfiler::RegionTimer reg (timer);
 
-     /*
+    if(facenr==0)
-     sei.SetSize (0);
-     for (SurfaceElementIndex i = 0; i < GetNSE(); i++)
     {
-        if ( (*this)[i].GetIndex () == facenr && (*this)[i][0] >= PointIndex::BASE &&
+        sei.SetSize(GetNSE());
-           !(*this)[i].IsDeleted() )
+        ParallelForRange( IntRange(GetNSE()), [&sei] (auto myrange)
             {
-           sei.Append (i);
+                for(auto i : myrange)
+                    sei[i] = i;
+            });
+        return;
     }
-     }
-     */
 
-     /* Philippose - 01/10/2009
-     Commented out the following lines, and activated the originally 
-     commented out lines above because of a bug which causes corruption 
-     of the variable "facedecoding" when a mesh is converted to second order
-     */
-
-     //      int size1 = sei.Size();
      sei.SetSize(0);
 
      SurfaceElementIndex si = facedecoding[facenr-1].firstelement;
@@ -6444,16 +6436,6 @@ namespace netgen
 
         si = (*this)[si].next;
      }
-     
-     /*
-     // *testout << "with list = " << endl << sei << endl;
-
-     if (size1 != sei.Size()) 
-     {
-        cout << "size mismatch" << endl;
-        exit(1);
-     }
-     */
   }