New function Mesh::FindIllegalTrigs()

Find illegal trigs after surface meshing and use this information in
surface mesh optimization

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();
 

libsrc/meshing/meshclass.cpp

@@ -3689,31 +3689,51 @@ 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");
-        auto & tab = const_cast<decltype(illegal_trigs)&>(illegal_trigs);
-        tab = make_unique<INDEX_3_CLOSED_HASHTABLE<int>> (3*GetNSE() + 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(tab->Used(i3) && tab->Get(i3)!=sei)
-                tab -> Set (i3, -1);
-            else
-                tab -> Set (i3, sei);
-          }
-
-      }
     INDEX_3 i3 (el[0], el[1], el[2]);
     i3.Sort();
-    if(illegal_trigs->Used(i3) && illegal_trigs->Get(i3)==-1)
+    if(illegal_trigs->Used(i3))
         return false;
 
     return 1;

libsrc/meshing/meshclass.hpp

@@ -571,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:

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");

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();