diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp
index b00ecc34..3bc1c440 100644
--- a/libsrc/csg/genmesh.cpp
+++ b/libsrc/csg/genmesh.cpp
@@ -504,6 +504,8 @@ namespace netgen
 	for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++)
 	  mesh[sei].SetIndex (k);
 
+        auto n_illegal_trigs = mesh.FindIllegalTrigs();
+        PrintMessage (3, n_illegal_trigs, " illegal triangles");
 
 	//      mesh.CalcSurfacesOfNode();
 
diff --git a/libsrc/meshing/improve2.cpp b/libsrc/meshing/improve2.cpp
index 261ada54..2fc9f5f4 100644
--- a/libsrc/meshing/improve2.cpp
+++ b/libsrc/meshing/improve2.cpp
@@ -430,6 +430,7 @@ namespace netgen
   {
     if (!faceindex)
       {
+        SplitImprove(mesh);
 	PrintMessage (3, "Combine improve");
 
 	for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
@@ -694,7 +695,11 @@ namespace netgen
 		  if ( (normals[el[l]] * nv) < 0.5)
 		    bad2 += 1e10;
 
-		illegal2 += 1-mesh.LegalTrig(el);
+                Element2d el1 = el;
+                for (auto i : Range(3))
+                    if(el1[i]==pi2)
+                        el1[i] = pi1;
+                illegal2 += 1-mesh.LegalTrig(el1);
 	      }
 	    bad2 /= hasonepi.Size();
 
@@ -789,6 +794,142 @@ namespace netgen
     mesh.SetNextTimeStamp();
   }
 
+  void MeshOptimize2d :: SplitImprove (Mesh & mesh)
+  {
+    if (!faceindex)
+      {
+        PrintMessage (3, "Split improve");
+
+        mesh.CalcSurfacesOfNode(); // TODO: needed?
+        for (faceindex = 1; faceindex <= mesh.GetNFD(); faceindex++)
+          {
+            SplitImprove (mesh);
+
+            if (multithread.terminate)
+                throw NgException ("Meshing stopped");
+          }
+
+        faceindex = 0;
+        mesh.Compress(); // TODO: needed?
+        return;
+      }
+
+    Array<SurfaceElementIndex> elements;
+    mesh.GetSurfaceElementsOfFace (faceindex, elements);
+
+    // return if we have quads in this surface
+    for (auto & ei : elements)
+        if (mesh[ei].GetNP() != 3)
+            return;
+
+    // maps from edges to adjacent trigs
+    INDEX_2_HASHTABLE<tuple<SurfaceElementIndex, SurfaceElementIndex>> els_on_edge(2*elements.Size() + 2);
+
+    // build els_on_edge table
+    for (SurfaceElementIndex sei : elements)
+      {
+        const Element2d & sel = mesh[sei];
+
+        for (int j = 0; j < 3; j++)
+          {
+            PointIndex pi1 = sel.PNumMod(j+2);
+            PointIndex pi2 = sel.PNumMod(j+3);
+
+            if (mesh.IsSegment (pi1, pi2)) continue;
+
+            INDEX_2 ii2 (pi1, pi2);
+            ii2.Sort();
+            if (els_on_edge.Used (ii2))
+              {
+                auto els = els_on_edge.Get(ii2);
+                get<1>(els) = sei;
+                els_on_edge.Set(ii2, els);
+              }
+            else
+              {
+                els_on_edge.Set (ii2, make_tuple(sei, sei));
+              }
+          }
+      }
+
+    // split edges of illegal trigs
+    for (SurfaceElementIndex sei : elements)
+      {
+        Element2d & sel = mesh[sei];
+
+        if (sel.IsDeleted()) continue;
+
+        // TODO: split also bad trigs, nut just illegal ones
+        if (mesh.LegalTrig(sel)) continue;
+
+        for (int j = 0; j < 3; j++)
+          {
+            PointIndex pi1 = sel.PNumMod(j+2);
+            PointIndex pi2 = sel.PNumMod(j+3);
+            PointIndex pi3 = sel.PNumMod(j+1);
+            PointIndex pi4;
+            PointGeomInfo gi1 = sel.GeomInfoPiMod(j+2);
+            PointGeomInfo gi2 = sel.GeomInfoPiMod(j+3);
+            PointGeomInfo gi3 = sel.GeomInfoPiMod(j+1);
+            PointGeomInfo gi4;
+
+            if (mesh.IsSegment (pi1, pi2)) continue;
+
+            // get neighbor element
+            INDEX_2 ii2 (pi1, pi2);
+            ii2.Sort();
+            auto els = els_on_edge.Get(ii2);
+            SurfaceElementIndex other_i = get<0>(els);
+            if(other_i==sei) other_i = get<1>(els);
+            auto & other = mesh[other_i];
+
+            // find opposite point of neighbor element
+            for (int j = 0; j < 3; j++)
+                if(other[j]!=pi1 && other[j]!=pi2)
+                  {
+                    pi4 = other[j];
+                    gi4 = other.GeomInfoPi(j);
+                    break;
+                  }
+
+            // split edge pi1,pi2
+            Point<3> p5;
+            PointIndex pi5;
+            PointGeomInfo gi5;
+
+            mesh.GetGeometry()->GetRefinement().PointBetween  (mesh[pi1], mesh[pi2], 0.5,
+                    faceindex,
+                    gi1, gi2, p5, gi5);
+
+            pi5 = mesh.AddPoint(p5);
+
+            Element2d e1(3);
+            e1.SetIndex(faceindex);
+            e1={ {pi1,gi1}, {pi5,gi5}, {pi3,gi3} };
+            mesh.AddSurfaceElement( e1 );
+
+            Element2d e2(3);
+            e2.SetIndex(faceindex);
+            e2 ={ {pi5,gi5}, {pi2,gi2}, {pi3,gi3} };
+            mesh.AddSurfaceElement( e2 );
+
+            Element2d e3(3);
+            e3.SetIndex(faceindex);
+            e3 ={ {pi1,gi1}, {pi4,gi4}, {pi5,gi5} };
+            mesh.AddSurfaceElement( e3 );
+
+            Element2d e4(3);
+            e4.SetIndex(faceindex);
+            e4 ={ {pi4,gi4}, {pi2,gi2}, {pi5,gi5} };
+            mesh.AddSurfaceElement( e4 );
+
+            sel.Delete();
+            other.Delete();
+          }
+      }
+
+    mesh.SetNextTimeStamp();
+  }
 
   void MeshOptimize2d :: CheckMeshApproximation (Mesh & mesh)
   {
diff --git a/libsrc/meshing/improve2.hpp b/libsrc/meshing/improve2.hpp
index e0e0c70a..c9fbae4b 100644
--- a/libsrc/meshing/improve2.hpp
+++ b/libsrc/meshing/improve2.hpp
@@ -24,6 +24,7 @@ public:
 
   void EdgeSwapping (Mesh & mesh, int usemetric);
   void CombineImprove (Mesh & mesh);
+  void SplitImprove (Mesh & mesh);
 
   void GenericImprove (Mesh & mesh);
 
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 0249ba4a..e521c590 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -3849,9 +3849,53 @@ namespace netgen
     return 0;
   }
 
+  int Mesh :: FindIllegalTrigs ()
+  {
+    // Temporary table to store the vertex numbers of all triangles
+    INDEX_3_CLOSED_HASHTABLE<int> temp_tab(3*GetNSE() + 1);
+    size_t cnt = 0;
+    for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
+      {
+        const Element2d & sel = surfelements[sei];
+        if (sel.IsDeleted()) continue;
+
+        INDEX_3 i3(sel[0], sel[1], sel[2]);
+        i3.Sort();
+        if(temp_tab.Used(i3))
+          {
+            temp_tab.Set (i3, -1);
+            cnt++;
+          }
+        else
+          {
+            temp_tab.Set (i3, sei);
+          }
+      }
+
+    illegal_trigs = make_unique<INDEX_3_CLOSED_HASHTABLE<int>> (2*cnt+1);
+    for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++)
+      {
+        const Element2d & sel = surfelements[sei];
+        if (sel.IsDeleted()) continue;
+
+        INDEX_3 i3(sel[0], sel[1], sel[2]);
+        i3.Sort();
+        if(temp_tab.Get(i3)==-1)
+            illegal_trigs -> Set (i3, 1);
+      }
+    return cnt;
+  }
 
   bool Mesh :: LegalTrig (const Element2d & el) const
   {
+    // Search for surface trigs with same vertices ( may happen for instance with close surfaces in stl geometies )
+    if(!illegal_trigs)
+        throw Exception("In Mesh::LegalTrig() - illegal_trigs table not built");
+    INDEX_3 i3 (el[0], el[1], el[2]);
+    i3.Sort();
+    if(illegal_trigs->Used(i3))
+        return false;
+
     return 1;
     if ( /* hp */ 1)  // needed for old, simple hp-refinement
       { 
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index e873efa9..8fc56c73 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -55,6 +55,7 @@ namespace netgen
     unique_ptr<INDEX_2_CLOSED_HASHTABLE<int>> segmentht;
     ///
     unique_ptr<INDEX_3_CLOSED_HASHTABLE<int>> surfelementht;
+    unique_ptr<INDEX_3_CLOSED_HASHTABLE<int>> illegal_trigs;
 
     /// faces of rest-solid
     NgArray<Element2d> openelements;
@@ -570,6 +571,10 @@ namespace netgen
 
 
     ///
+    // Find trigs with same vertices
+    // return: number of illegal trigs
+    int FindIllegalTrigs ();
+
     bool LegalTrig (const Element2d & el) const;
     /**
        if values non-null, return values in 4-double array:
diff --git a/libsrc/occ/occgenmesh.cpp b/libsrc/occ/occgenmesh.cpp
index 414b52f6..78c72cb2 100644
--- a/libsrc/occ/occgenmesh.cpp
+++ b/libsrc/occ/occgenmesh.cpp
@@ -862,6 +862,9 @@ namespace netgen
 
         for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++)
           mesh[sei].SetIndex (k);
+
+        auto n_illegal_trigs = mesh.FindIllegalTrigs();
+        PrintMessage (3, n_illegal_trigs, " illegal triangles");
       }
 
     //      ofstream problemfile("occmesh.rep");
diff --git a/libsrc/stlgeom/meshstlsurface.cpp b/libsrc/stlgeom/meshstlsurface.cpp
index eaf5a8af..5a3d9ff5 100644
--- a/libsrc/stlgeom/meshstlsurface.cpp
+++ b/libsrc/stlgeom/meshstlsurface.cpp
@@ -279,6 +279,9 @@ int STLSurfaceMeshing (STLGeometry & geom, class Mesh & mesh, const MeshingParam
 	      mesh.FindOpenSegments();
 	      nopen = mesh.GetNOpenSegments();
 
+              auto n_illegal_trigs = mesh.FindIllegalTrigs();
+              PrintMessage (3, n_illegal_trigs, " illegal triangles");
+
 	      if (nopen)
 		{
 		  geom.ClearMarkedSegs();