Merge remote-tracking branch 'origin/master' into test_refactoring_meshing_design

libsrc/core/exception.cpp

@@ -167,6 +167,11 @@ namespace ngcore
 
 static void ngcore_signal_handler(int sig)
 {
+  static bool first_call = true;
+  if(!first_call)
+      exit(1); // avoid endless recursions if signals are caused by this handler
+  first_call = false;
+
   switch(sig)
     {
       case SIGABRT:

libsrc/core/logging.cpp

@@ -15,7 +15,7 @@ namespace ngcore
 {
   std::ostream* testout = new std::ostream(nullptr); // NOLINT
 
-  level::level_enum Logger::global_level;
+  level::level_enum Logger::global_level = level::warn;
 
   void Logger::log(level::level_enum level, std::string && s)
   {

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/core/utils.cpp

@@ -15,10 +15,22 @@ namespace ngcore
   // windows does demangling in typeid(T).name()
   NGCORE_API std::string Demangle(const char* typeinfo) { return typeinfo; }
 #else
-  NGCORE_API std::string Demangle(const char* typeinfo) { int status; return abi::__cxa_demangle(typeinfo,
+  NGCORE_API std::string Demangle(const char* typeinfo)
-                                                                                      nullptr,
+  {
-                                                                                      nullptr,
+    int status=0;
-                                                                                      &status); }
+    try
+      {
+        char *s = abi::__cxa_demangle(typeinfo, nullptr, nullptr, &status);
+        std::string result{s};
+        free(s);
+        return result;
+      }
+    catch( const std::exception & e )
+      {
+        GetLogger("utils")->warn("{}:{} cannot demangle {}, status: {}, error:{}", __FILE__, __LINE__, typeinfo, status, e.what());
+      }
+    return typeinfo;
+  }
 #endif
 
   double seconds_per_tick = [] () noexcept

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/delaunay.cpp

@@ -1,9 +1,289 @@
 #include <mystdlib.h>
 #include "meshing.hpp"
 
-
 namespace netgen
 {
+  template<int dim, typename T=INDEX, typename TSCAL=double>
+  class DelaunayTree
+  {
+  public:
+    // Number of entries per leaf
+    static constexpr int N = 100;
+
+    struct Node;
+
+    struct Leaf
+    {
+      Point<2*dim, TSCAL> p[N];
+      T index[N];
+      int n_elements;
+      int nr;
+
+      Leaf() : n_elements(0)
+      { }
+
+
+      void Add( Array<Leaf*> &leaves, Array<T> &leaf_index, const Point<2*dim> &ap, T aindex )
+        {
+          p[n_elements] = ap;
+          index[n_elements] = aindex;
+          n_elements++;
+          if(leaf_index.Size()<aindex+1)
+              leaf_index.SetSize(aindex+1);
+          leaf_index[aindex] = nr;
+        }
+    };
+
+    struct Node
+    {
+      union
+      {
+        Node *children[2];
+        Leaf *leaf;
+      };
+      double sep;
+      int level;
+
+      Node()
+          : children{nullptr,nullptr}
+      { }
+
+      ~Node()
+       { }
+
+      Leaf *GetLeaf() const
+        {
+          return children[1] ? nullptr : leaf;
+        }
+    };
+
+  private:
+    Node root;
+
+    Array<Leaf*> leaves;
+    Array<T> leaf_index;
+
+    Point<dim> global_min, global_max;
+    double tol;
+    size_t n_leaves;
+    size_t n_nodes;
+    BlockAllocator ball_nodes;
+    BlockAllocator ball_leaves;
+
+  public:
+
+    DelaunayTree (const Point<dim> & pmin, const Point<dim> & pmax)
+        : global_min(pmin), global_max(pmax), n_leaves(1), n_nodes(1), ball_nodes(sizeof(Node)), ball_leaves(sizeof(Leaf))
+      {
+        root.leaf = (Leaf*) ball_leaves.Alloc(); new (root.leaf) Leaf();
+        root.leaf->nr = 0;
+        leaves.Append(root.leaf);
+        root.level = 0;
+        tol = 1e-7 * Dist(pmax, pmin);
+      }
+
+    DelaunayTree (const Box<dim> & box)
+        : DelaunayTree(box.PMin(), box.PMax())
+      { }
+
+    size_t GetNLeaves()
+      {
+        return n_leaves;
+      }
+
+    size_t GetNNodes()
+      {
+        return n_nodes;
+      }
+
+    template<typename TFunc>
+    void GetFirstIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
+            TFunc func=[](auto pi){return false;}) const
+      {
+        // static Timer timer("DelaunayTree::GetIntersecting"); RegionTimer rt(timer);
+        // static Timer timer1("DelaunayTree::GetIntersecting-LinearSearch");
+        ArrayMem<const Node*, 100> stack;
+        ArrayMem<int, 100> dir_stack;
+
+
+        Point<2*dim> tpmin, tpmax;
+
+        for (size_t i : IntRange(dim))
+          {
+            tpmin(i) = global_min(i);
+            tpmax(i) = pmax(i)+tol;
+
+            tpmin(i+dim) = pmin(i)-tol;
+            tpmax(i+dim) = global_max(i);
+          }
+
+        stack.SetSize(0);
+        stack.Append(&root);
+        dir_stack.SetSize(0);
+        dir_stack.Append(0);
+
+        while(stack.Size())
+          {
+            const Node *node = stack.Last();
+            stack.DeleteLast();
+
+            int dir = dir_stack.Last();
+            dir_stack.DeleteLast();
+
+            if(Leaf *leaf = node->GetLeaf())
+              {
+                //               RegionTimer rt1(timer1);
+                for (auto i : IntRange(leaf->n_elements))
+                  {
+                    bool intersect = true;
+                    const auto p = leaf->p[i];
+
+                    for (int d = 0; d < dim; d++)
+                        if (p[d] > tpmax[d])
+                            intersect = false;
+                    for (int d = dim; d < 2*dim; d++)
+                        if (p[d] < tpmin[d])
+                            intersect = false;
+                    if(intersect)
+                        if(func(leaf->index[i])) return;
+                  }
+              }
+            else
+              {
+                int newdir = dir+1;
+                if(newdir==2*dim) newdir = 0;
+                if (tpmin[dir] <= node->sep)
+                  {
+                    stack.Append(node->children[0]);
+                    dir_stack.Append(newdir);
+                  }
+                if (tpmax[dir] >= node->sep)
+                  {
+                    stack.Append(node->children[1]);
+                    dir_stack.Append(newdir);
+                  }
+              }
+          }
+      }
+
+    void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
+            NgArray<T> & pis) const
+      {
+        pis.SetSize(0);
+        GetFirstIntersecting(pmin, pmax, [&pis](auto pi) { pis.Append(pi); return false;});
+      }
+
+    void Insert (const Box<dim> & box, T pi)
+      {
+        Insert (box.PMin(), box.PMax(), pi);
+      }
+
+    void Insert (const Point<dim> & pmin, const Point<dim> & pmax, T pi)
+      {
+        // static Timer timer("DelaunayTree::Insert"); RegionTimer rt(timer);
+        int dir = 0;
+        Point<2*dim> p;
+        for (auto i : IntRange(dim))
+          {
+            p(i) = pmin[i];
+            p(i+dim) = pmax[i];
+          }
+
+        Node * node = &root;
+        Leaf * leaf = node->GetLeaf();
+
+        // search correct leaf to add point
+        while(!leaf)
+          {
+            node = p[dir] < node->sep ? node->children[0] : node->children[1];
+            dir++;
+            if(dir==2*dim) dir = 0;
+            leaf = node->GetLeaf();
+          }
+
+        // add point to leaf
+        if(leaf->n_elements < N)
+            leaf->Add(leaves, leaf_index, p,pi);
+        else // assume leaf->n_elements == N
+          {
+            // add two new nodes and one new leaf
+            int n_elements = leaf->n_elements;
+            ArrayMem<TSCAL, N> coords(n_elements);
+            ArrayMem<int, N> order(n_elements);
+
+            // separate points in two halves, first sort all coordinates in direction dir
+            for (auto i : IntRange(n_elements))
+              {
+                order[i] = i;
+                coords[i] = leaf->p[i][dir];
+              }
+
+            QuickSortI(coords, order);
+            int isplit = N/2;
+            Leaf *leaf1 = (Leaf*) ball_leaves.Alloc(); new (leaf1) Leaf();
+            Leaf *leaf2 = (Leaf*) ball_leaves.Alloc(); new (leaf2) Leaf();
+
+            leaf1->nr = leaf->nr;
+            leaf2->nr = leaves.Size();
+            leaves.Append(leaf2);
+            leaves[leaf1->nr] = leaf1;
+
+            for (auto i : order.Range(isplit))
+                leaf1->Add(leaves, leaf_index, leaf->p[i], leaf->index[i] );
+            for (auto i : order.Range(isplit, N))
+                leaf2->Add(leaves, leaf_index, leaf->p[i], leaf->index[i] );
+
+            Node *node1 = (Node*) ball_nodes.Alloc(); new (node1) Node();
+            node1->leaf = leaf1;
+            node1->level = node->level+1;
+
+            Node *node2 = (Node*) ball_nodes.Alloc(); new (node2) Node();
+            node2->leaf = leaf2;
+            node2->level = node->level+1;
+
+            node->children[0] = node1;
+            node->children[1] = node2;
+            node->sep = 0.5 * (leaf->p[order[isplit-1]][dir] + leaf->p[order[isplit]][dir]);
+
+            // add new point to one of the new leaves
+            if (p[dir] < node->sep)
+                leaf1->Add( leaves, leaf_index, p, pi );
+            else
+                leaf2->Add( leaves, leaf_index, p, pi );
+
+            ball_leaves.Free(leaf);
+            n_leaves++;
+            n_nodes+=2;
+          }
+      }
+
+    void DeleteElement (T pi)
+      {
+        // static Timer timer("DelaunayTree::DeleteElement"); RegionTimer rt(timer);
+        Leaf *leaf = leaves[leaf_index[pi]];
+        leaf_index[pi] = -1;
+        auto & n_elements = leaf->n_elements;
+        auto & index = leaf->index;
+        auto & p = leaf->p;
+
+        for (auto i : IntRange(n_elements))
+          {
+            if(index[i] == pi)
+              {
+                n_elements--;
+                if(i!=n_elements)
+                  {
+                    index[i] = index[n_elements];
+                    p[i] = p[n_elements];
+                  }
+                return;
+              }
+          }
+      }
+  };
+
+  // typedef BoxTree<3> DTREE;
+  typedef DelaunayTree<3> DTREE;
 
   static const int deltetfaces[][3] = 
     { { 1, 2, 3 },
@@ -227,14 +507,14 @@ namespace netgen
   void AddDelaunayPoint (PointIndex newpi, const Point3d & newp, 
 			 NgArray<DelaunayTet> & tempels, 
 			 Mesh & mesh,
-			 BoxTree<3> & tettree,
+			 DTREE & tettree,
 			 MeshNB & meshnb,
 			 NgArray<Point<3> > & centers, NgArray<double> & radi2,
 			 NgArray<int> & connected, NgArray<int> & treesearch, 
 			 NgArray<int> & freelist, SphereList & list,
 			 IndexSet & insphere, IndexSet & closesphere)
   {
-    static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t);
+    static Timer t("Meshing3::AddDelaunayPoint");// RegionTimer reg(t);
     // static Timer tsearch("addpoint, search");
     // static Timer tinsert("addpoint, insert");
 
@@ -635,7 +915,7 @@ namespace netgen
     pmin2 = pmin2 + 0.1 * (pmin2 - pmax2);
     pmax2 = pmax2 + 0.1 * (pmax2 - pmin2);
 
-    BoxTree<3> tettree(pmin2, pmax2);
+    DTREE tettree(pmin2, pmax2);
 
 
     tempels.Append (startel);
@@ -798,6 +1078,7 @@ namespace netgen
 	  tempmesh.AddVolumeElement (el);
 	}
 
+      tempels.DeleteAll();
 
       MeshQuality3d (tempmesh);
 
@@ -852,6 +1133,7 @@ namespace netgen
     
       MeshQuality3d (tempmesh);
     
+      tempels.SetSize(tempmesh.GetNE());
       tempels.SetSize(0);
       for (auto & el : tempmesh.VolumeElements())
         tempels.Append (el);

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,239 +19,25 @@ namespace netgen
   };
 
 
-
+  bool MeshOptimize2d :: EdgeSwapping (const int usemetric,
- 
+    Array<Neighbour> &neighbors,
-  void MeshOptimize2d :: EdgeSwapping (int usemetric)
+    Array<bool> &swapped,
-  {
+    const SurfaceElementIndex t1, const int o1,
-    static Timer timer("EdgeSwapping (2D)"); RegionTimer reg(timer);
+    const int t,
-    if (!faceindex)
+    Array<int,PointIndex> &pdef,
-      {
+    const bool check_only)
-	if (usemetric)
-	  PrintMessage (3, "Edgeswapping, metric");
-	else
-	  PrintMessage (3, "Edgeswapping, topological");
-
-	for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
-	  {
-	    EdgeSwapping (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);
-
-    /*
-    for (int i = 0; i < seia.Size(); i++)
-      if (mesh[seia[i]].GetNP() != 3)
-	{
-	  GenericImprove (mesh);
-	  return;
-	}
-    */
-    for (SurfaceElementIndex sei : seia)
-      if (mesh[sei].GetNP() != 3)
-	{
-	  GenericImprove();
-	  return;
-	}
-      
-    int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
-
-    Array<Neighbour> neighbors(mesh.GetNSE());
-    INDEX_2_HASHTABLE<trionedge> other(2*seia.Size() + 2);
-
-
-    Array<bool> swapped(mesh.GetNSE());
-    NgArray<int,PointIndex::BASE> pdef(mesh.GetNP());
-    NgArray<double,PointIndex::BASE> 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++)
-      {
-	const Element2d & sel = mesh[seia[i]];
-	for (int j = 0; j < 3; j++)
-	  pangle[sel[j]] = 0.0;
-      }
-    // pangle = 0;
-
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & sel = mesh[seia[i]];
-	for (int j = 0; j < 3; j++)
-	  {
-	    POINTTYPE typ = mesh[sel[j]].Type();
-	    if (typ == FIXEDPOINT || typ == EDGEPOINT)
-	      {
-		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; 
-    // 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);
-          neighbors[sei].SetOrientation (j, 0);
-        }
-
-    /*
-      NgArray<Vec3d> normals(mesh.GetNP());
-      for (i = 1; i <= mesh.GetNSE(); i++)
-      {
-      Element2d & hel = mesh.SurfaceElement(i);
-      if (hel.GetIndex() == faceindex)
-      for (k = 1; k <= 3; k++)
-      {
-      int pi = hel.PNum(k);
-      SelectSurfaceOfPoint (mesh.Point(pi), hel.GeomInfoPi(k));
-      int surfi = mesh.GetFaceDescriptor(faceindex).SurfNr();
-      GetNormalVector (surfi, mesh.Point(pi), normals.Elem(pi));
-      normals.Elem(pi) /= normals.Elem(pi).Length();
-      }
-      }
-    */	    
-
-    /*
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	const Element2d & sel = mesh[seia[i]];
-    */
-
-    for (SurfaceElementIndex sei : seia)
-      {
-	const Element2d & sel = mesh[sei];
-        
-	for (int j = 0; j < 3; j++)
-	  {
-	    PointIndex pi1 = sel.PNumMod(j+2);
-	    PointIndex pi2 = sel.PNumMod(j+3);
-	  
-	    //	    double loch = mesh.GetH(mesh[pi1]);
-	    
-	    // INDEX_2 edge(pi1, pi2);
-	    // edge.Sort();
-	  
-	    if (mesh.IsSegment (pi1, pi2))
-	      continue;
-	  
-	    /*
-	      if (segments.Used (edge))
-	      continue;
-	    */
-	    INDEX_2 ii2 (pi1, pi2);
-	    if (other.Used (ii2))
-	      {
-		// INDEX_2 i2s(ii2);
-		// i2s.Sort();
-
-                /*
-		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;
-    
-    NgProfiler::StopTimer (timerstart);
-  
-
-
-    int t = 4;
-    bool done = false;
-    while (!done && t >= 2)
-      {
-	for (int i = 0; i < seia.Size(); i++)
-	  {
-	    SurfaceElementIndex t1 = seia[i];
-
-	    if (mesh[t1].IsDeleted())
-	      continue;
-
-	    if (mesh[t1].GetIndex() != faceindex)
-	      continue;
-
-	    if (multithread.terminate)
-	      throw NgException ("Meshing stopped");
-
-	    for (int o1 = 0; o1 < 3; o1++)
   {
     bool should;
+    bool do_swap = false;
 
     SurfaceElementIndex t2 = neighbors[t1].GetNr (o1);
     int o2 = neighbors[t1].GetOrientation (o1);
 
-		if (t2 == -1) continue;
+    if (t2 == -1) return false;
-		if (swapped[t1] || swapped[t2]) continue;
+    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);
@@ -289,7 +57,7 @@ namespace netgen
     allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4;
 
     if(!allowswap)
-		  continue;
+        return false;
 
     // normal of new
     Vec<3> nv1 = Cross (auxvec1, auxvec2);
@@ -300,7 +68,7 @@ namespace netgen
 
 
     if(!allowswap)
-		  continue;
+        return false;
 
     Vec<3> nv2 = Cross (auxvec1, auxvec2);
 
@@ -317,12 +85,10 @@ namespace netgen
     nv1.Normalize();
     nv2.Normalize();
 
-		// SelectSurfaceOfPoint (mesh.Point(pi3), gi3);
     auto nvp3 = geo.GetNormal (surfnr, mesh.Point(pi3), gi3);
 
     nvp3.Normalize();
 
-		// SelectSurfaceOfPoint (mesh.Point(pi4), gi4);
     auto nvp4 = geo.GetNormal (surfnr, mesh.Point(pi4), gi4);
 
     nvp4.Normalize();
@@ -380,24 +146,11 @@ namespace netgen
             if (legal2 < legal1) should = false;
           }
 
-		    if (should)
+        do_swap = should;
+        if (should && !check_only)
           {
             // 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 } };
 
@@ -410,8 +163,187 @@ namespace netgen
             swapped[t2] = true;
           }
       }
+    return do_swap;
+  }
+
+ 
+  void MeshOptimize2d :: EdgeSwapping (int usemetric)
+  {
+    static Timer timer("EdgeSwapping (2D)"); RegionTimer reg(timer);
+    static Timer timer_nb("EdgeSwapping-Find neighbors");
+    if (usemetric)
+      PrintMessage (3, "Edgeswapping, metric");
+    else
+      PrintMessage (3, "Edgeswapping, topological");
+
+    static Timer timerstart("EdgeSwapping 2D start");
+    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)
+                  mixed = true;
+            });
+      }
+    else
+      {
+        mesh.GetSurfaceElementsOfFace (faceindex, seia);
+        for (SurfaceElementIndex sei : seia)
+            if (mesh[sei].GetNP() != 3)
+                mixed = true;
+      }
+
+    if(mixed)
+        return GenericImprove();
+      
+    Array<Neighbour> neighbors(mesh.GetNSE());
+    auto elements_on_node = mesh.CreatePoint2SurfaceElementTable(faceindex);
+
+    Array<bool> swapped(mesh.GetNSE());
+    Array<int,PointIndex> pdef(mesh.GetNP());
+    Array<double,PointIndex> pangle(mesh.GetNP());
+
+    static const double minangle[] = { 0, 1.481, 2.565, 3.627, 4.683, 5.736, 7, 9 };
+
+
+    if(faceindex == 0)
+      {
+        ParallelFor( Range(pangle), [&] (auto i) NETGEN_LAMBDA_INLINE
+            {
+              pangle[i] = 0.0;
+            });
+      }
+    else
+      {
+        ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
+            {
+              const Element2d & sel = mesh[seia[i]];
+              for (int j = 0; j < 3; j++)
+                  pangle[sel[j]] = 0.0;
+            });
+      }
+
+    ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
+        {
+          const Element2d & sel = mesh[seia[i]];
+          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
+        {
+          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;
+            }
+       });
+
+    ParallelFor( Range(seia), [this, &pdef, &neighbors, &seia, &elements_on_node] (auto i) NETGEN_LAMBDA_INLINE
+        {
+          auto sei = 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);
+            }
+
+          const auto sel = mesh[sei];
+          for (int j = 0; j < 3; j++)
+            {
+              PointIndex pi1 = sel.PNumMod(j+2);
+              PointIndex pi2 = sel.PNumMod(j+3);
+
+              for (auto sei_other : elements_on_node[pi1])
+                {
+                  if(sei_other==sei) continue;
+                  const auto & other = mesh[sei_other];
+                  int pi1_other = -1;
+                  int pi2_other = -1;
+                  bool common_edge = false;
+                  for (int k = 0; k < 3; k++)
+                    {
+                      if(other[k] == pi1)
+                          pi1_other = k;
+                      if(other[k] == pi2)
+                        {
+                          pi2_other = k;
+                          common_edge = true;
+                        }
+                    }
+
+                  if(common_edge)
+                    {
+                      neighbors[sei].SetNr (j, sei_other);
+                      neighbors[sei].SetOrientation (j, 3-pi1_other-pi2_other);
+                    }
+                }
+            }
+        });
+
+
+    for (SurfaceElementIndex sei : seia)
+      swapped[sei] = false;
+    
+    timerstart.Stop();
+  
+
+
+    Array<std::pair<SurfaceElementIndex,int>> improvement_candidates(3*seia.Size());
+    atomic<int> cnt(0);
+
+    int t = 4;
+    bool done = false;
+    while (!done && t >= 2)
+      {
+        cnt = 0;
+        ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
+          {
+            SurfaceElementIndex t1 = seia[i];
+
+            if (mesh[t1].IsDeleted())
+              return;
+
+            if (mesh[t1].GetIndex() != faceindex)
+              return;
+
+            if (multithread.terminate)
+              throw NgException ("Meshing stopped");
+
+            for (int o1 = 0; o1 < 3; o1++)
+                if(EdgeSwapping(usemetric, neighbors, swapped, t1, o1, t, pdef, true))
+                    improvement_candidates[cnt++]= std::make_pair(t1,o1);
+          });
+
+        auto elements_with_improvement = improvement_candidates.Range(cnt.load());
+        QuickSort(elements_with_improvement);
+
+        for (auto [t1,o1] : elements_with_improvement)
+            done |= EdgeSwapping(usemetric, neighbors, swapped, t1, o1, t, pdef, false);
 	t--;
       }
 
@@ -421,160 +353,31 @@ namespace netgen
 
 
 
-
+  double CombineImproveEdge( Mesh & mesh,
-
+                           const Table<SurfaceElementIndex, PointIndex> & elementsonnode,
-
+                           Array<Vec<3>, PointIndex> & normals,
- 
+                           Array<bool, PointIndex> & fixed,
-  void MeshOptimize2d :: CombineImprove()
+                           PointIndex pi1, PointIndex pi2,
+                           bool check_only = true)
   {
-    if (!faceindex)
-      {
-        SplitImprove();
-	PrintMessage (3, "Combine improve");
-
-	for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
-	  {
-	    CombineImprove();
-
-	    if (multithread.terminate)
-	      throw NgException ("Meshing stopped");
-	  }
-	faceindex = 0;
-	return;
-      }
-
-
-    static int timer = NgProfiler::CreateTimer ("Combineimprove 2D");
-    NgProfiler::RegionTimer reg (timer);
-
-    static int timerstart = NgProfiler::CreateTimer ("Combineimprove 2D start");
-    NgProfiler::StartTimer  (timerstart);
-
-
-    static int timerstart1 = NgProfiler::CreateTimer ("Combineimprove 2D start1");
-    NgProfiler::StartTimer  (timerstart1);
-
-    
-    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();
-
-    TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonnode(np); 
-    Array<SurfaceElementIndex> hasonepi, hasbothpi;
-
-    for (SurfaceElementIndex sei : seia)
-      for (PointIndex pi : mesh[sei].PNums<3>())
-        elementsonnode.Add (pi, sei);
-
-    Array<bool,PointIndex> fixed(np);
-    fixed = false;
-
-    NgProfiler::StopTimer  (timerstart1);
-
-    /*
-    for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
-      {
-	INDEX_2 i2(mesh[si][0], mesh[si][1]);
-	fixed[i2.I1()] = true;
-	fixed[i2.I2()] = true;
-      }
-    */
-
-    for (SurfaceElementIndex sei : seia)
-      {
-	Element2d & sel = mesh[sei];
-	for (int j = 0; j < sel.GetNP(); j++)
-	  {
-	    PointIndex pi1 = sel.PNumMod(j+2);
-	    PointIndex pi2 = sel.PNumMod(j+3);
-	    if (mesh.IsSegment (pi1, pi2))
-	      {	
-		fixed[pi1] = true;
-		fixed[pi2] = true;
-	      }
-	  }
-      }
-
-
-    /*
-    for(int i = 0; i < mesh.LockedPoints().Size(); i++)
-      fixed[mesh.LockedPoints()[i]] = true;
-    */
-    for (PointIndex pi : mesh.LockedPoints())
-      fixed[pi] = true;
-
-
-    Array<Vec<3>,PointIndex> normals(np);
-
-    // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
-    for (PointIndex pi : mesh.Points().Range())
-      {
-	if (elementsonnode[pi].Size())
-	  {
-	    Element2d & hel = mesh[elementsonnode[pi][0]];
-	    for (int k = 0; k < 3; k++)
-	      if (hel[k] == pi)
-		{
-		  // SelectSurfaceOfPoint (mesh[pi], hel.GeomInfoPi(k+1));
-		  normals[pi] = geo.GetNormal(surfnr, mesh[pi], hel.GeomInfoPi(k+1));
-		  break;
-		}
-	  }
-      }
-
-    NgProfiler::StopTimer  (timerstart);
-
-    for (int i = 0; i < seia.Size(); i++)
-      {
-	SurfaceElementIndex sei = seia[i];
-	Element2d & elem = mesh[sei];
-
-	for (int j = 0; j < 3; j++)
-	  {
-            if (elem.IsDeleted()) continue;
-	    PointIndex pi1 = elem[j];
-	    PointIndex pi2 = elem[(j+1) % 3];
-
-            /*
-	    if (pi1 < PointIndex::BASE || 
-		pi2 < PointIndex::BASE)
-	      continue;
-            */
     if (!pi1.IsValid() || !pi2.IsValid())
-	      continue;              
+        return 0.0;
-	    /*
-	      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 0.0;
     if (pi2 < pi1) swap (pi1, pi2);
     */
 
@@ -583,23 +386,10 @@ namespace netgen
         Swap (pi1, pi2);
 
     if (fixed[pi2])
-	      continue;
+        return 0.0;
 
     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];
@@ -618,18 +408,11 @@ namespace netgen
           }
       }
 
+    if(hasbothpi.Size()==0)
+        return 0.0;
 
-	    Element2d & hel = mesh[hasbothpi[0]];
-	    for (int k = 0; k < 3; k++)
-	      if (hel[k] == pi1)
-		{
-		  // SelectSurfaceOfPoint (mesh[pi1],
-                  //			hel.GeomInfoPi(k+1));
-		  nv = geo.GetNormal(surfnr, mesh[pi1], hel.GeomInfoPi(k+1));
-		  break;
-		}
 
-	    //	  nv = normals.Get(pi1);
+    nv = normals[pi1];
 
 
     for (SurfaceElementIndex sei2 :  elementsonnode[pi2])
@@ -667,17 +450,14 @@ namespace netgen
       }
     bad1 /= (hasonepi.Size()+hasbothpi.Size());
 
-	    MeshPoint p1 = mesh[pi1];
-	    MeshPoint p2 = mesh[pi2];
-
-	    MeshPoint pnew = p1;
-	    mesh[pi1] = pnew;
-	    mesh[pi2] = pnew;
-
     double bad2 = 0;
     for (int k = 0; k < hasonepi.Size(); k++)
       {
-		Element2d & el = mesh[hasonepi[k]];
+        Element2d el = mesh[hasonepi[k]];
+        for (auto i : Range(3))
+            if(el[i]==pi2)
+                el[i] = pi1;
+
         double err =
             CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]],
                     nv, -1, loch);
@@ -691,32 +471,31 @@ namespace netgen
             bad2 += 1e10;
 
         for (int l = 0; l < 3; l++)
+          {
             if ( (normals[el[l]] * nv) < 0.5)
                 bad2 += 1e10;
+          }
 
-                Element2d el1 = el;
+        illegal2 += 1-mesh.LegalTrig(el);
-                for (auto i : Range(3))
-                    if(el1[i]==pi2)
-                        el1[i] = pi1;
-                illegal2 += 1-mesh.LegalTrig(el1);
       }
     bad2 /= hasonepi.Size();
 
-	    mesh[pi1] = p1;
-	    mesh[pi2] = p2;
-       
     if (debugflag)
       {
         (*testout) << "bad1 = " << bad1 << ", bad2 = " << bad2 << endl;
       }
 
-	    bool should = (bad2 < bad1 && bad2 < 1e4);
+    bool should = (illegal2<=illegal1 && bad2 < bad1 && bad2 < 1e4);
-	    if (bad2 < 1e4)
+    if(illegal2 < illegal1)
       {
-		if (illegal1 > illegal2) should = true;
+        should = true;
-		if (illegal2 > illegal1) should = false;
+        bad1 += 1e4;
       }
 
+    double d_badness = should * (bad2-bad1);
+
+    if(check_only)
+        return d_badness;
 
     if (should)
       {
@@ -727,7 +506,6 @@ namespace netgen
            (*testout) << "loch = " << loch << endl;
            */
 
-		mesh[pi1] = pnew;
         PointGeomInfo gi;
         // bool gi_set(false);
 
@@ -760,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])
@@ -769,8 +546,6 @@ namespace netgen
             if (el.IsDeleted()) continue;
             if (el.PNums().Contains(pi1)) continue;
 
-		    elementsonnode.Add (pi1, sei2);
-                    
             for (auto l : Range(el.GetNP()))
               {
                 if (el[l] == pi2)
@@ -785,8 +560,121 @@ namespace netgen
 
         for (auto sei : hasbothpi)
             mesh[sei].Delete();
+
+      }
+    return d_badness;
+  }
+
+  void MeshOptimize2d :: CombineImprove ()
+  {
+    SplitImprove();
+    PrintMessage (3, "Combine improve");
+
+    if (multithread.terminate)
+        throw NgException ("Meshing stopped");
+
+    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;
+
+    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;
+    ParallelFor( Range(seia), [&] (auto i) NETGEN_LAMBDA_INLINE
+            {
+                if (mesh[seia[i]].GetNP() != 3)
+                    mixed = true;
+            });
+
+    if(mixed)
+        return;
+
+    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)
+                    {
+                      const int faceindex = hel.GetIndex();
+                      const int surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
+                      normals[pi] = geo.GetNormal (surfnr, mesh[pi], hel.GeomInfoPi(k+1));
+                      break;
                     }
               }
+        }, TasksPerThread(4));
+
+    timerstart.Stop();
+
+    // Find edges with improvement
+    Array<std::tuple<double, int>> candidate_edges(edges.Size());
+    std::atomic<int> improvement_counter(0);
+
+    ParallelFor( Range(edges), [&] (auto i) NETGEN_LAMBDA_INLINE
+      {
+        auto [pi1, pi2] = edges[i];
+        double d_badness = CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi1, pi2, true);
+        if(d_badness < 0.0)
+            candidate_edges[improvement_counter++] = make_tuple(d_badness, i);
+      }, TasksPerThread(4));
+
+    auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
+    QuickSort(edges_with_improvement);
+
+    for(auto [d_badness, ei] : edges_with_improvement)
+      {
+        auto [pi1, pi2] = edges[ei];
+        CombineImproveEdge(mesh, elementsonnode, normals, fixed, pi1, pi2, false);
       }
 
     //  mesh.Compress();

libsrc/meshing/improve2.hpp

@@ -1,7 +1,76 @@
 #ifndef FILE_IMPROVE2
 #define FILE_IMPROVE2
 
+template<typename TINDEX>
+void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges )
+{
+  static Timer tbuild_edges("Build edges"); RegionTimer reg(tbuild_edges);
 
+  static constexpr int tetedges[6][2] =
+    { { 0, 1 }, { 0, 2 }, { 0, 3 },
+        { 1, 2 }, { 1, 3 }, { 2, 3 } };
+
+  int ntasks = 2*ngcore::TaskManager::GetMaxThreads();
+  Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);
+
+  ParallelFor(IntRange(ntasks), [&] (int ti)
+    {
+      auto myrange = mesh.Points().Range().Split(ti, ntasks);
+      ArrayMem<std::tuple<PointIndex,PointIndex>, 100> local_edges;
+      for (auto pi : myrange)
+      {
+        local_edges.SetSize(0);
+
+        for(auto ei : elementsonnode[pi])
+        {
+            const auto & elem = mesh[ei];
+            if (elem.IsDeleted()) continue;
+
+            for (int j = 0; j < 6; j++)
+            {
+                PointIndex pi0 = elem[tetedges[j][0]];
+                PointIndex pi1 = elem[tetedges[j][1]];
+                if (pi1 < pi0) Swap(pi0, pi1);
+                if(pi0==pi)
+                    local_edges.Append(std::make_tuple(pi0, pi1));
+            }
+        }
+        QuickSort(local_edges);
+
+        auto edge_prev = std::make_tuple<PointIndex, PointIndex>(-1,-1);
+
+        for(auto edge : local_edges)
+            if(edge != edge_prev)
+            {
+                task_edges[ti].Append(edge);
+                edge_prev = edge;
+            }
+      }
+    }, ntasks);
+
+  int num_edges = 0;
+  for (auto & edg : task_edges)
+      num_edges += edg.Size();
+  edges.SetAllocSize(num_edges);
+  for (auto & edg : task_edges)
+      edges.Append(edg);
+}
+
+
+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
@@ -24,6 +93,8 @@ public:
   void ProjectBoundaryPoints(NgArray<int> & surfaceindex, 
 			     const NgArray<Point<3>* > & from, NgArray<Point<3>* > & dest);
 
+  bool EdgeSwapping (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 (int usemetric);
   void CombineImprove ();
   void SplitImprove ();

libsrc/meshing/improve2gen.cpp

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

libsrc/meshing/improve3.cpp

@@ -409,60 +409,6 @@ void MeshOptimize3d :: CombineImproveSequential (Mesh & mesh,
   multithread.task = savetask;
 } 
 
-void MeshOptimize3d :: BuildEdgeList( const Mesh & mesh, const Table<ElementIndex, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges )
-{
-  static Timer tbuild_edges("Build edges"); RegionTimer reg(tbuild_edges);
-
-  static constexpr int tetedges[6][2] =
-    { { 0, 1 }, { 0, 2 }, { 0, 3 },
-      { 1, 2 }, { 1, 3 }, { 2, 3 } };
-
-  int ntasks = 2*ngcore::TaskManager::GetMaxThreads();
-  Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);
-
-  ParallelFor(IntRange(ntasks), [&] (int ti)
-  {
-    auto myrange = mesh.Points().Range().Split(ti, ntasks);
-    ArrayMem<std::tuple<PointIndex,PointIndex>, 100> local_edges;
-    for (auto pi : myrange)
-    {
-      local_edges.SetSize(0);
-
-      for(auto ei : elementsonnode[pi])
-      {
-          const Element & elem = mesh[ei];
-          if (elem.IsDeleted()) continue;
-
-          for (int j = 0; j < 6; j++)
-          {
-              PointIndex pi0 = elem[tetedges[j][0]];
-              PointIndex pi1 = elem[tetedges[j][1]];
-              if (pi1 < pi0) Swap(pi0, pi1);
-              if(pi0==pi)
-                  local_edges.Append(std::make_tuple(pi0, pi1));
-          }
-      }
-      QuickSort(local_edges);
-
-      auto edge_prev = std::make_tuple<PointIndex, PointIndex>(-1,-1);
-
-      for(auto edge : local_edges)
-          if(edge != edge_prev)
-          {
-              task_edges[ti].Append(edge);
-              edge_prev = edge;
-          }
-    }
-  }, ntasks);
-
-  int num_edges = 0;
-  for (auto & edg : task_edges)
-      num_edges += edg.Size();
-  edges.SetAllocSize(num_edges);
-  for (auto & edg : task_edges)
-      edges.Append(edg);
-}
-
 void MeshOptimize3d :: CombineImprove (Mesh & mesh,
 				       OPTIMIZEGOAL goal)
 {

libsrc/meshing/improve3.hpp

@@ -12,8 +12,6 @@ class MeshOptimize3d
 {
   const MeshingParameters & mp;
 
-  void BuildEdgeList( const Mesh & mesh, const Table<ElementIndex, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges );
-
 public:
   MeshOptimize3d (const MeshingParameters & amp) : mp(amp) { ; }
 

libsrc/meshing/meshclass.cpp

@@ -6205,14 +6205,19 @@ 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());
+
+    if(faceindex==0)
+      {
         for ( ; !creator.Done(); creator++)
           ngcore::ParallelForRange
             (Range(surfelements), [&] (auto myrange)
@@ -6220,7 +6225,21 @@ namespace netgen
                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/meshing/msghandler.cpp

@@ -3,7 +3,7 @@
 
 namespace netgen
 {
-int printmessage_importance = 5;
+int printmessage_importance = 3;
 int printwarnings = 1;
 int printerrors = 1;
 int printdots = 1;

libsrc/meshing/smoothing2.cpp

@@ -690,35 +690,23 @@ namespace netgen
 
   void MeshOptimize2d :: ImproveMesh (const MeshingParameters & mp)
   {
-    if (!faceindex)
+    static Timer timer("MeshSmoothing 2D"); RegionTimer reg (timer);
-      {
+
     PrintMessage (3, "Smoothing");
 
-	for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
-	  {
-	    ImproveMesh (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)
             {
@@ -726,13 +714,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++)
           {
@@ -750,58 +732,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;
@@ -872,27 +848,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);
 
     /*
@@ -902,28 +857,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++)
@@ -942,9 +908,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();
 		
@@ -971,7 +937,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]];
@@ -989,20 +957,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);
 
@@ -1011,6 +974,8 @@ namespace netgen
 	  double fact = 1;
 	  int moveisok = 0;
 
+          if(mixed)
+            {
               // restore other points
               for (int j = 0; j < ld.locelements.Size(); j++)
                 {
@@ -1021,6 +986,7 @@ namespace netgen
                       if (hhpi != pi) mesh[hhpi] = savepoints[hhpi];
                     }
                 }
+            }
 
 	  
 	  //optimizer loop (if whole distance is not possible, move only a bit!!!!)
@@ -1061,11 +1027,10 @@ namespace netgen
 	    }
           }
         }
+      }, mixed ? 1 : ngcore::TasksPerThread(4)); // mixed element smoothing not parallel yet
+    }
 
     tloop.Stop();
-    if (printeddot)
-      PrintDot ('\n');
-
     CheckMeshApproximation (mesh);
     mesh.SetNextTimeStamp();
   }

libsrc/occ/occmeshsurf.cpp

@@ -371,6 +371,7 @@ namespace netgen
   void OCCSurface :: Project (Point<3> & ap, PointGeomInfo & gi)
   {
     static Timer t("OccSurface::Project"); RegionTimer reg(t);
+    static Timer t2("OccSurface::Project actural"); 
 
 
     // try Newton's method ...
@@ -475,7 +476,10 @@ namespace netgen
     Handle( ShapeAnalysis_Surface ) su = new ShapeAnalysis_Surface( occface );
     auto toltool =  BRep_Tool::Tolerance( topods_face );
 
-    gp_Pnt2d suval = su->ValueOfUV ( pnt, toltool);
+    // gp_Pnt2d suval = su->ValueOfUV ( pnt, toltool);
+    t2.Start();
+    gp_Pnt2d suval = su->NextValueOfUV (gp_Pnt2d(u,v), pnt, toltool);
+    t2.Stop();
     suval.Coord( u, v);
     pnt = occface->Value( u, v );
     

libsrc/stlgeom/meshstlsurface.cpp

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

libsrc/visualization/vsmesh.cpp

@@ -1029,6 +1029,9 @@ namespace netgen
 	else
 	  glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, matcol);
 	
+        static Point<3> xa[129];
+        static Vec<3> na[129];
+
         
         
 	for (int hi = 0; hi < seia.Size(); hi++)
@@ -1058,8 +1061,6 @@ namespace netgen
                   if (curv.IsHighOrder()) //  && curv.IsSurfaceElementCurved(sei))
 		    {
 		      if (hoplotn > 128) hoplotn = 128;
-		      Point<3> xa[129];
-		      Vec<3> na[129];
 
 		      for (int i = 0; i < hoplotn; i++)
 			{

libsrc/visualization/vssolution.cpp

@@ -29,6 +29,11 @@ namespace netgen
     // vssolution.AddUserVisualizationObject (vis);
     GetVSSolution().AddUserVisualizationObject (vis);
   }
+  void DeleteUserVisualizationObject (UserVisualizationObject * vis)
+  {
+    // vssolution.AddUserVisualizationObject (vis);
+    GetVSSolution().DeleteUserVisualizationObject (vis);
+  }
 
 
   VisualSceneSolution :: SolData :: SolData ()

libsrc/visualization/vssolution.hpp

@@ -233,7 +233,12 @@ public:
   {
     user_vis.Append (vis);
   }
-  
+  void DeleteUserVisualizationObject (UserVisualizationObject * vis)
+  {
+    int pos = user_vis.Pos(vis);
+    if (pos >= 0)
+      user_vis.Delete(pos);
+  }
 
 private:
   void GetClippingPlaneTrigs (NgArray<ClipPlaneTrig> & trigs, NgArray<ClipPlanePoint> & pts);

tests/pytest/test_tutorials.py

@@ -58,6 +58,10 @@ 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 == "cylsphere.geo":
+        return standard[0:2] + standard[3:] # coarse gives inconsistent reults (other mesh on MacOS)
+    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")]