Merge branch 'parallel_edgeswapping' into 'master'

Parallel edgeswapping

See merge request jschoeberl/netgen!274

libsrc/general/parthreads.hpp

@@ -96,12 +96,6 @@ void ParallelFor( int first, int next, const TFunc & f )
 
 
   
-  template<typename T>
-  inline atomic<T> & AsAtomic (T & d)
-  {
-    return reinterpret_cast<atomic<T>&> (d);
-  }
-
   typedef void (*TaskManager)(std::function<void(int,int)>);
   typedef void (*Tracer)(string, bool);  // false .. start, true .. stop
 

libsrc/meshing/improve2.cpp

@@ -7,24 +7,6 @@
 namespace netgen
 {
 
-  class Neighbour
-  {
-    int nr[3];
-    int orient[3];
-
-  public:
-    Neighbour () { ; } 
-
-    void SetNr (int side, int anr) { nr[side] = anr; }
-    int GetNr (int side) { return nr[side]; }
-
-    void SetOrientation (int side, int aorient) { orient[side] = aorient; }
-    int GetOrientation (int side) { return orient[side]; }
-  };
-
-
-
-
   class trionedge
   {
   public:
@@ -37,204 +19,293 @@ namespace netgen
   };
 
 
+  bool MeshOptimize2d :: EdgeSwapping (Mesh & mesh, const int usemetric,
+    Array<Neighbour> &neighbors,
+    Array<bool> &swapped,
+    const SurfaceElementIndex t1, const int o1,
+    const int t,
+    Array<int,PointIndex> &pdef,
+    const bool check_only)
+  {
+    bool should;
+    bool do_swap = false;
+
+    SurfaceElementIndex t2 = neighbors[t1].GetNr (o1);
+    int o2 = neighbors[t1].GetOrientation (o1);
+
+    if (t2 == -1) return false;
+    if (swapped[t1] || swapped[t2]) return false;
+
+    const int faceindex = mesh[t1].GetIndex();
+    const int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
+
+    PointIndex pi1 = mesh[t1].PNumMod(o1+1+1);
+    PointIndex pi2 = mesh[t1].PNumMod(o1+1+2);
+    PointIndex pi3 = mesh[t1].PNumMod(o1+1);
+    PointIndex pi4 = mesh[t2].PNumMod(o2+1);
+
+    PointGeomInfo gi1 = mesh[t1].GeomInfoPiMod(o1+1+1);
+    PointGeomInfo gi2 = mesh[t1].GeomInfoPiMod(o1+1+2);
+    PointGeomInfo gi3 = mesh[t1].GeomInfoPiMod(o1+1);
+    PointGeomInfo gi4 = mesh[t2].GeomInfoPiMod(o2+1);
+
+    bool allowswap = true;
+
+    Vec<3> auxvec1 = mesh[pi3]-mesh[pi4];
+    Vec<3> auxvec2 = mesh[pi1]-mesh[pi4];
+
+    allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4;
+
+    if(!allowswap)
+        return false;
+
+    // normal of new
+    Vec<3> nv1 = Cross (auxvec1, auxvec2);
+
+    auxvec1 = mesh.Point(pi4)-mesh.Point(pi3);
+    auxvec2 = mesh.Point(pi2)-mesh.Point(pi3);
+    allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4;
+
+
+    if(!allowswap)
+        return false;
+
+    Vec<3> nv2 = Cross (auxvec1, auxvec2);
+
+
+    // normals of original
+    Vec<3> nv3 = Cross (mesh[pi1]-mesh[pi4], mesh[pi2]-mesh[pi4]);
+    Vec<3> nv4 = Cross (mesh[pi2]-mesh[pi3], mesh[pi1]-mesh[pi3]);
+
+    nv3 *= -1;
+    nv4 *= -1;
+    nv3.Normalize();
+    nv4.Normalize();
+
+    nv1.Normalize();
+    nv2.Normalize();
+
+    Vec<3> nvp3, nvp4;
+    GetNormalVector (surfnr, mesh.Point(pi3), gi3, nvp3);
+
+    nvp3.Normalize();
+
+    GetNormalVector (surfnr, mesh.Point(pi4), gi4, nvp4);
+
+    nvp4.Normalize();
+
+
+
+    double critval = cos (M_PI / 6);  // 30 degree
+    allowswap = allowswap &&
+        (nv1 * nvp3 > critval) &&
+        (nv1 * nvp4 > critval) &&
+        (nv2 * nvp3 > critval) &&
+        (nv2 * nvp4 > critval) &&
+        (nvp3 * nv3 > critval) &&
+        (nvp4 * nv4 > critval);
+
+
+    double horder = Dist (mesh[pi1], mesh[pi2]);
+
+    if ( // nv1 * nv2 >= 0 &&
+            nv1.Length() > 1e-3 * horder * horder &&
+            nv2.Length() > 1e-3 * horder * horder &&
+            allowswap )
+      {
+        if (!usemetric)
+          {
+            int e = pdef[pi1] + pdef[pi2] - pdef[pi3] - pdef[pi4];
+            double d =
+                Dist2 (mesh[pi1], mesh[pi2]) -
+                Dist2 (mesh[pi3], mesh[pi4]);
+
+            should = e >= t && (e > 2 || d > 0);
+          }
+        else
+          {
+            double loch = mesh.GetH(mesh[pi1]);
+            should =
+                CalcTriangleBadness (mesh[pi4], mesh[pi3], mesh[pi1], metricweight, loch) +
+                CalcTriangleBadness (mesh[pi3], mesh[pi4], mesh[pi2], metricweight, loch) <
+                CalcTriangleBadness (mesh[pi1], mesh[pi2], mesh[pi3], metricweight, loch) +
+                CalcTriangleBadness (mesh[pi2], mesh[pi1], mesh[pi4], metricweight, loch);
+          }
+
+        if (allowswap)
+          {
+            Element2d sw1 (pi4, pi3, pi1);
+            Element2d sw2 (pi3, pi4, pi2);
+
+            int legal1 =
+                mesh.LegalTrig (mesh[t1]) +
+                mesh.LegalTrig (mesh[t2]);
+            int legal2 =
+                mesh.LegalTrig (sw1) + mesh.LegalTrig (sw2);
+
+            if (legal1 < legal2) should = true;
+            if (legal2 < legal1) should = false;
+          }
+
+        do_swap = should;
+        if (should && !check_only)
+          {
+            // do swapping !
+
+            mesh[t1] = { { pi1, gi1 }, { pi4, gi4 }, { pi3, gi3 } };
+            mesh[t2] = { { pi2, gi2 }, { pi3, gi3 }, { pi4, gi4 } };
+
+            pdef[pi1]--;
+            pdef[pi2]--;
+            pdef[pi3]++;
+            pdef[pi4]++;
+
+            swapped[t1] = true;
+            swapped[t2] = true;
+          }
+      }
+    return do_swap;
+  }
 
  
   void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
   {
     static Timer timer("EdgeSwapping (2D)"); RegionTimer reg(timer);
-    if (!faceindex)
+    static Timer timer_nb("EdgeSwapping-Find neighbors");
-      {
+    if (usemetric)
-	if (usemetric)
+      PrintMessage (3, "Edgeswapping, metric");
-	  PrintMessage (3, "Edgeswapping, metric");
+    else
-	else
+      PrintMessage (3, "Edgeswapping, topological");
-	  PrintMessage (3, "Edgeswapping, topological");
-
-	for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
-	  {
-	    EdgeSwapping (mesh, usemetric);
-
-	    if (multithread.terminate)
-	      throw NgException ("Meshing stopped");
-	  }
-
-	faceindex = 0;
-	mesh.CalcSurfacesOfNode();
-	return;
-      }
-
 
     static int timerstart = NgProfiler::CreateTimer ("EdgeSwapping 2D start");
     NgProfiler::StartTimer (timerstart);
 
-
     Array<SurfaceElementIndex> seia;
-    mesh.GetSurfaceElementsOfFace (faceindex, seia);
+    bool mixed = false;
 
-    /*
+    if(faceindex==0)
-    for (int i = 0; i < seia.Size(); i++)
+      {
-      if (mesh[seia[i]].GetNP() != 3)
+        seia.SetSize(mesh.GetNSE());
-	{
+        ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
-	  GenericImprove (mesh);
+            {
-	  return;
+              SurfaceElementIndex sei(i);
-	}
+              seia[i] = sei;
-    */
+              if (mesh[sei].GetNP() != 3)
-    for (SurfaceElementIndex sei : seia)
+                  mixed = true;
-      if (mesh[sei].GetNP() != 3)
+            });
-	{
+      }
-	  GenericImprove (mesh);
+    else
-	  return;
+      {
-	}
+        mesh.GetSurfaceElementsOfFace (faceindex, seia);
+        for (SurfaceElementIndex sei : seia)
+            if (mesh[sei].GetNP() != 3)
+                mixed = true;
+      }
+
+    if(mixed)
+        return GenericImprove(mesh);
       
-    int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
-
     Array<Neighbour> neighbors(mesh.GetNSE());
-    INDEX_2_HASHTABLE<trionedge> other(2*seia.Size() + 2);
+    auto elements_on_node = mesh.CreatePoint2SurfaceElementTable(faceindex);
-
 
     Array<bool> swapped(mesh.GetNSE());
-    NgArray<int,PointIndex::BASE> pdef(mesh.GetNP());
+    Array<int,PointIndex> pdef(mesh.GetNP());
-    NgArray<double,PointIndex::BASE> pangle(mesh.GetNP());
+    Array<double,PointIndex> pangle(mesh.GetNP());
 
-
-    // int e;
-    // double d;
-    // Vec3d nv1, nv2;
-
-    // double loch(-1);
     static const double minangle[] = { 0, 1.481, 2.565, 3.627, 4.683, 5.736, 7, 9 };
 
 
-    for (int i = 0; i < seia.Size(); i++)
+    if(faceindex == 0)
       {
-	const Element2d & sel = mesh[seia[i]];
+        ParallelFor( Range(pangle), [&] (auto i) NETGEN_LAMBDA_INLINE
-	for (int j = 0; j < 3; j++)
+            {
-	  pangle[sel[j]] = 0.0;
+              pangle[i] = 0.0;
+            });
       }
-    // pangle = 0;
+    else
-
-    for (int i = 0; i < seia.Size(); i++)
       {
-	const Element2d & sel = mesh[seia[i]];
+        ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
-	for (int j = 0; j < 3; j++)
+            {
-	  {
+              const Element2d & sel = mesh[seia[i]];
-	    POINTTYPE typ = mesh[sel[j]].Type();
+              for (int j = 0; j < 3; j++)
-	    if (typ == FIXEDPOINT || typ == EDGEPOINT)
+                  pangle[sel[j]] = 0.0;
-	      {
+            });
-		pangle[sel[j]] +=
-		  Angle (mesh[sel[(j+1)%3]] - mesh[sel[j]],
-			 mesh[sel[(j+2)%3]] - mesh[sel[j]]);
-	      }
-	  }
       }
 
-    // for (PointIndex pi = PointIndex::BASE; 
+    ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
-    // pi < mesh.GetNP()+PointIndex::BASE; pi++)
-    
-    // pdef = 0;
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & sel = mesh[seia[i]];
-	for (int j = 0; j < 3; j++)
-	  {
-	    PointIndex pi = sel[j];
-	    if (mesh[pi].Type() == INNERPOINT || mesh[pi].Type() == SURFACEPOINT)
-	      pdef[pi] = -6;
-	    else
-	      for (int j = 0; j < 8; j++)
-		if (pangle[pi] >= minangle[j])
-		  pdef[pi] = -1-j;
-	  }
-      }
-
-    /*
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & sel = mesh[seia[i]];
-	for (int j = 0; j < 3; j++)
-	  pdef[sel[j]]++;
-      }
-    */
-    for (SurfaceElementIndex sei : seia)
-      for (PointIndex pi : mesh[sei].PNums<3>())
-        pdef[pi]++;
-    
-    // for (int i = 0; i < seia.Size(); i++)
-    for (SurfaceElementIndex sei : seia)
-      for (int j = 0; j < 3; j++)
         {
-          neighbors[sei].SetNr (j, -1);
+          const Element2d & sel = mesh[seia[i]];
-          neighbors[sei].SetOrientation (j, 0);
+          for (int j = 0; j < 3; j++)
-        }
+            {
+              POINTTYPE typ = mesh[sel[j]].Type();
+              if (typ == FIXEDPOINT || typ == EDGEPOINT)
+                {
+                  AtomicAdd(pangle[sel[j]],
+                    Angle (mesh[sel[(j+1)%3]] - mesh[sel[j]],
+                           mesh[sel[(j+2)%3]] - mesh[sel[j]]));
+                }
+            }
+       });
 
-    /*
+    ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
-      NgArray<Vec3d> normals(mesh.GetNP());
+        {
-      for (i = 1; i <= mesh.GetNSE(); i++)
+          const Element2d & sel = mesh[seia[i]];
-      {
+          for (int j = 0; j < 3; j++)
-      Element2d & hel = mesh.SurfaceElement(i);
+            {
-      if (hel.GetIndex() == faceindex)
+              PointIndex pi = sel[j];
-      for (k = 1; k <= 3; k++)
+              if (mesh[pi].Type() == INNERPOINT || mesh[pi].Type() == SURFACEPOINT)
-      {
+                pdef[pi] = -6;
-      int pi = hel.PNum(k);
+              else
-      SelectSurfaceOfPoint (mesh.Point(pi), hel.GeomInfoPi(k));
+                for (int j = 0; j < 8; j++)
-      int surfi = mesh.GetFaceDescriptor(faceindex).SurfNr();
+                  if (pangle[pi] >= minangle[j])
-      GetNormalVector (surfi, mesh.Point(pi), normals.Elem(pi));
+                    pdef[pi] = -1-j;
-      normals.Elem(pi) /= normals.Elem(pi).Length();
+            }
-      }
+       });
-      }
-    */	    
 
-    /*
+    ParallelFor( Range(seia), [&pdef, &neighbors, &mesh, &seia, &elements_on_node] (auto i) NETGEN_LAMBDA_INLINE
-    for (int i = 0; i < seia.Size(); i++)
+        {
-      {
+          auto sei = seia[i];
-	const Element2d & sel = mesh[seia[i]];
+          for (PointIndex pi : mesh[sei].template PNums<3>())
-    */
+              AsAtomic(pdef[pi])++;
+          for (int j = 0; j < 3; j++)
+            {
+              neighbors[sei].SetNr (j, -1);
+              neighbors[sei].SetOrientation (j, 0);
+            }
 
-    for (SurfaceElementIndex sei : seia)
+          const auto sel = mesh[sei];
-      {
+          for (int j = 0; j < 3; j++)
-	const Element2d & sel = mesh[sei];
+            {
-        
+              PointIndex pi1 = sel.PNumMod(j+2);
-	for (int j = 0; j < 3; j++)
+              PointIndex pi2 = sel.PNumMod(j+3);
-	  {
+
-	    PointIndex pi1 = sel.PNumMod(j+2);
+              for (auto sei_other : elements_on_node[pi1])
-	    PointIndex pi2 = sel.PNumMod(j+3);
+                {
-	  
+                  if(sei_other==sei) continue;
-	    //	    double loch = mesh.GetH(mesh[pi1]);
+                  const auto & other = mesh[sei_other];
-	    
+                  int pi1_other = -1;
-	    // INDEX_2 edge(pi1, pi2);
+                  int pi2_other = -1;
-	    // edge.Sort();
+                  bool common_edge = false;
-	  
+                  for (int k = 0; k < 3; k++)
-	    if (mesh.IsSegment (pi1, pi2))
+                    {
-	      continue;
+                      if(other[k] == pi1)
-	  
+                          pi1_other = k;
-	    /*
+                      if(other[k] == pi2)
-	      if (segments.Used (edge))
+                        {
-	      continue;
+                          pi2_other = k;
-	    */
+                          common_edge = true;
-	    INDEX_2 ii2 (pi1, pi2);
+                        }
-	    if (other.Used (ii2))
+                    }
-	      {
+
-		// INDEX_2 i2s(ii2);
+                  if(common_edge)
-		// i2s.Sort();
+                    {
+                      neighbors[sei].SetNr (j, sei_other);
+                      neighbors[sei].SetOrientation (j, 3-pi1_other-pi2_other);
+                    }
+                }
+            }
+        });
 
-                /*
-		int i2 = other.Get(ii2).tnr;
-		int j2 = other.Get(ii2).sidenr;
-		*/
-                auto othertrig = other.Get(ii2);
-		SurfaceElementIndex i2 = othertrig.tnr;
-		int j2 = othertrig.sidenr;
-                
-		neighbors[sei].SetNr (j, i2);
-		neighbors[sei].SetOrientation (j, j2);
-		neighbors[i2].SetNr (j2, sei);
-		neighbors[i2].SetOrientation (j2, j);
-	      }
-	    else
-	      {
-		other.Set (INDEX_2 (pi2, pi1), trionedge (sei, j));
-	      }
-	  }
-      }
 
     for (SurfaceElementIndex sei : seia)
       swapped[sei] = false;
@@ -243,174 +314,37 @@ namespace netgen
   
 
 
+    Array<std::pair<SurfaceElementIndex,int>> improvement_candidates(3*seia.Size());
+    atomic<int> cnt(0);
+
     int t = 4;
     bool done = false;
     while (!done && t >= 2)
       {
-	for (int i = 0; i < seia.Size(); i++)
+        cnt = 0;
-	  {
+        ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
-	    SurfaceElementIndex t1 = seia[i];
+          {
+            SurfaceElementIndex t1 = seia[i];
 
-	    if (mesh[t1].IsDeleted())
+            if (mesh[t1].IsDeleted())
-	      continue;
+              return;
 
-	    if (mesh[t1].GetIndex() != faceindex)
+            if (mesh[t1].GetIndex() != faceindex)
-	      continue;
+              return;
 
-	    if (multithread.terminate)
+            if (multithread.terminate)
-	      throw NgException ("Meshing stopped");
+              throw NgException ("Meshing stopped");
 
-	    for (int o1 = 0; o1 < 3; o1++)
+            for (int o1 = 0; o1 < 3; o1++)
-	      {
+                if(EdgeSwapping(mesh, usemetric, neighbors, swapped, t1, o1, t, pdef, true))
-		bool should;
+                    improvement_candidates[cnt++]= std::make_pair(t1,o1);
+          });
 
-		SurfaceElementIndex t2 = neighbors[t1].GetNr (o1);
+        auto elements_with_improvement = improvement_candidates.Range(cnt.load());
-		int o2 = neighbors[t1].GetOrientation (o1);
+        QuickSort(elements_with_improvement);
 
-		if (t2 == -1) continue;
+        for (auto [t1,o1] : elements_with_improvement)
-		if (swapped[t1] || swapped[t2]) continue;
+            done |= EdgeSwapping(mesh, usemetric, neighbors, swapped, t1, o1, t, pdef, false);
-	      
-
-		PointIndex pi1 = mesh[t1].PNumMod(o1+1+1);
-		PointIndex pi2 = mesh[t1].PNumMod(o1+1+2);
-		PointIndex pi3 = mesh[t1].PNumMod(o1+1);
-		PointIndex pi4 = mesh[t2].PNumMod(o2+1);
-	      
-		PointGeomInfo gi1 = mesh[t1].GeomInfoPiMod(o1+1+1);
-		PointGeomInfo gi2 = mesh[t1].GeomInfoPiMod(o1+1+2);
-		PointGeomInfo gi3 = mesh[t1].GeomInfoPiMod(o1+1);
-		PointGeomInfo gi4 = mesh[t2].GeomInfoPiMod(o2+1);
-	    
-		bool allowswap = true;
-
-		Vec<3> auxvec1 = mesh[pi3]-mesh[pi4];
-		Vec<3> auxvec2 = mesh[pi1]-mesh[pi4];
-
-		allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4;
-
-		if(!allowswap)
-		  continue;
-
-		// normal of new
-		Vec<3> nv1 = Cross (auxvec1, auxvec2);
-
-		auxvec1 = mesh.Point(pi4)-mesh.Point(pi3);
-		auxvec2 = mesh.Point(pi2)-mesh.Point(pi3);
-		allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4;
-
-
-		if(!allowswap)
-		  continue;
-
-		Vec<3> nv2 = Cross (auxvec1, auxvec2);
-
-	      
-		// normals of original
-		Vec<3> nv3 = Cross (mesh[pi1]-mesh[pi4], mesh[pi2]-mesh[pi4]);
-		Vec<3> nv4 = Cross (mesh[pi2]-mesh[pi3], mesh[pi1]-mesh[pi3]);
-	      
-		nv3 *= -1;
-		nv4 *= -1;
-		nv3.Normalize();
-		nv4.Normalize();
-
-		nv1.Normalize();
-		nv2.Normalize();
-	    
-		Vec<3> nvp3, nvp4;
-		GetNormalVector (surfnr, mesh.Point(pi3), gi3, nvp3);
-
-		nvp3.Normalize();
-
-		GetNormalVector (surfnr, mesh.Point(pi4), gi4, nvp4);
-	    
-		nvp4.Normalize();
-	      
-	      
-	      
-		double critval = cos (M_PI / 6);  // 30 degree
-		allowswap = allowswap &&
-		  (nv1 * nvp3 > critval) && 
-		  (nv1 * nvp4 > critval) && 
-		  (nv2 * nvp3 > critval) && 
-		  (nv2 * nvp4 > critval) &&
-		  (nvp3 * nv3 > critval) && 
-		  (nvp4 * nv4 > critval);
-	      
-
-		double horder = Dist (mesh[pi1], mesh[pi2]);
-
-		if ( // nv1 * nv2 >= 0 &&
-		    nv1.Length() > 1e-3 * horder * horder &&
-		    nv2.Length() > 1e-3 * horder * horder &&
-		    allowswap )
-		  {
-		    if (!usemetric)
-		      {
-			int e = pdef[pi1] + pdef[pi2] - pdef[pi3] - pdef[pi4];
-			double d = 
-			  Dist2 (mesh[pi1], mesh[pi2]) - 
-			  Dist2 (mesh[pi3], mesh[pi4]);
-		      
-			should = e >= t && (e > 2 || d > 0);
-		      }
-		    else
-		      {
-			double loch = mesh.GetH(mesh[pi1]);
-			should = 
-			  CalcTriangleBadness (mesh[pi4], mesh[pi3], mesh[pi1], metricweight, loch) +
-			  CalcTriangleBadness (mesh[pi3], mesh[pi4], mesh[pi2], metricweight, loch) <
-			  CalcTriangleBadness (mesh[pi1], mesh[pi2], mesh[pi3], metricweight, loch) +
-			  CalcTriangleBadness (mesh[pi2], mesh[pi1], mesh[pi4], metricweight, loch);
-		      }
-		  
-		    if (allowswap)
-		      {
-			Element2d sw1 (pi4, pi3, pi1);
-			Element2d sw2 (pi3, pi4, pi2);
-
-			int legal1 = 
-			  mesh.LegalTrig (mesh[t1]) + 
-			  mesh.LegalTrig (mesh[t2]);
-			int legal2 = 
-			  mesh.LegalTrig (sw1) + mesh.LegalTrig (sw2);
-
-			if (legal1 < legal2) should = true;
-			if (legal2 < legal1) should = false;
-		      }
-		  
-		    if (should)
-		      {
-			// do swapping !
-		      
-			done = true;
-
-                        /*
-                        mesh[t1] = { pi1, pi4, pi3 };
-                        mesh[t2] = { pi2, pi3, pi4 };
-                        
-			mesh[t1].GeomInfoPi(1) = gi1;
-			mesh[t1].GeomInfoPi(2) = gi4;
-			mesh[t1].GeomInfoPi(3) = gi3;
-		      
-			mesh[t2].GeomInfoPi(1) = gi2;
-			mesh[t2].GeomInfoPi(2) = gi3;
-			mesh[t2].GeomInfoPi(3) = gi4;
-                        */
-                        mesh[t1] = { { pi1, gi1 }, { pi4, gi4 }, { pi3, gi3 } };
-                        mesh[t2] = { { pi2, gi2 }, { pi3, gi3 }, { pi4, gi4 } };
-                        
-			pdef[pi1]--;
-			pdef[pi2]--;
-			pdef[pi3]++;
-			pdef[pi4]++;
-		      
-			swapped[t1] = true;
-			swapped[t2] = true;
-		      }
-		  }
-	      }
-	  }
 	t--;
       }
 

libsrc/meshing/improve2.hpp

@@ -2,6 +2,20 @@
 #define FILE_IMPROVE2
 
 
+class Neighbour
+{
+  int nr[3];
+  int orient[3];
+
+public:
+  Neighbour () { ; }
+
+  void SetNr (int side, int anr) { nr[side] = anr; }
+  int GetNr (int side) { return nr[side]; }
+
+  void SetOrientation (int side, int aorient) { orient[side] = aorient; }
+  int GetOrientation (int side) { return orient[side]; }
+};
 
 ///
 class MeshOptimize2d
@@ -22,6 +36,8 @@ public:
   void ProjectBoundaryPoints(NgArray<int> & surfaceindex, 
 			     const NgArray<Point<3>* > & from, NgArray<Point<3>* > & dest);
 
+  bool EdgeSwapping (Mesh & mesh, const int usemetric, Array<Neighbour> &neighbors, Array<bool> &swapped,
+    const SurfaceElementIndex t1, const int edge, const int t, Array<int,PointIndex> &pdef, const bool check_only=false);
   void EdgeSwapping (Mesh & mesh, int usemetric);
   void CombineImprove (Mesh & mesh);
   void SplitImprove (Mesh & mesh);

libsrc/meshing/improve2gen.cpp

@@ -21,6 +21,7 @@ namespace netgen
 
   void MeshOptimize2d :: GenericImprove (Mesh & mesh)
   {
+    static Timer timer("MeshOptimize2d::GenericImprove"); RegionTimer reg(timer);
     if (!faceindex)
       {
 	if (writestatus)

libsrc/meshing/meshclass.cpp

@@ -6205,22 +6205,41 @@ namespace netgen
        {
          for (PointIndex pi : myrange)
            QuickSort(elementsonnode[pi]);
-       });
+       }, ngcore::TasksPerThread(4));
 
     return move(elementsonnode);
   }
 
-  Table<SurfaceElementIndex, PointIndex> Mesh :: CreatePoint2SurfaceElementTable() const
+  Table<SurfaceElementIndex, PointIndex> Mesh :: CreatePoint2SurfaceElementTable( int faceindex ) const
   {
+    static Timer timer("Mesh::CreatePoint2SurfaceElementTable"); RegionTimer rt(timer);
+
     TableCreator<SurfaceElementIndex, PointIndex> creator(GetNP());
-    for ( ; !creator.Done(); creator++)
+
-      ngcore::ParallelForRange
+    if(faceindex==0)
-        (Range(surfelements), [&] (auto myrange)
+      {
-         {
+        for ( ; !creator.Done(); creator++)
-           for (SurfaceElementIndex ei : myrange)
+          ngcore::ParallelForRange
-             for (PointIndex pi : (*this)[ei].PNums())
+            (Range(surfelements), [&] (auto myrange)
-               creator.Add (pi, ei);
+             {
-         });
+               for (SurfaceElementIndex ei : myrange)
+                 for (PointIndex pi : (*this)[ei].PNums())
+                   creator.Add (pi, ei);
+             }, ngcore::TasksPerThread(4));
+      }
+    else
+      {
+        Array<SurfaceElementIndex> face_els;
+        GetSurfaceElementsOfFace(faceindex, face_els);
+        for ( ; !creator.Done(); creator++)
+          ngcore::ParallelForRange
+            (Range(face_els), [&] (auto myrange)
+             {
+               for (auto i : myrange)
+                 for (PointIndex pi : (*this)[face_els[i]].PNums())
+                   creator.Add (pi, face_els[i]);
+             }, ngcore::TasksPerThread(4));
+      }
     
     auto elementsonnode = creator.MoveTable();
     ngcore::ParallelForRange

libsrc/meshing/meshclass.hpp

@@ -762,7 +762,7 @@ namespace netgen
   
 
     Table<ElementIndex, PointIndex> CreatePoint2ElementTable() const;
-    Table<SurfaceElementIndex, PointIndex> CreatePoint2SurfaceElementTable() const;
+    Table<SurfaceElementIndex, PointIndex> CreatePoint2SurfaceElementTable( int faceindex=0 ) const;
 
     DLL_HEADER bool PureTrigMesh (int faceindex = 0) const;
     DLL_HEADER bool PureTetMesh () const;

libsrc/stlgeom/meshstlsurface.cpp

@@ -306,7 +306,6 @@ int STLSurfaceMeshing (STLGeometry & geom, class Mesh & mesh, const MeshingParam
 		  
 		  mesh.CalcSurfacesOfNode();
 		  optmesh.EdgeSwapping (mesh, 0);
-		  mesh.CalcSurfacesOfNode();
 		  optmesh.ImproveMesh (mesh, mparam);
 		}
 

tests/pytest/test_tutorials.py

@@ -58,6 +58,8 @@ def getMeshingparameters(filename):
         return standard[:3] # this gets too big for finer meshsizes
     if filename == "screw.step":
         return standard[3:] # coarser meshes don't work here
+    if filename == "part1.stl":
+        return standard[0:1] + standard[2:] # very coarse does not work
     return standard
 
 _geofiles = [f for f in getFiles(".geo")] + [f for f in getFiles(".stl")]