Merge branch 'import_stl_as_mesh' into 'master'

Import stl as mesh

See merge request jschoeberl/netgen!240

libsrc/interface/readuser.cpp

@@ -7,6 +7,7 @@
 #include <myadt.hpp>
 #include <linalg.hpp>
 #include <csg.hpp>
+#include <stlgeom.hpp>
 #include <meshing.hpp>
 
 #include "writeuser.hpp"
@@ -648,6 +649,32 @@ namespace netgen
         ReadFNFFormat (mesh, filename);
       }
 
+    if ( ( (strlen (filename) > 4) && strcmp (&filename[strlen (filename)-4], ".stl") == 0 ) ||
+         ( (strlen (filename) > 5) && strcmp (&filename[strlen (filename)-5], ".stlb") == 0 ) )
+      {
+        ifstream ist{string{filename}};
+        auto geom = shared_ptr<STLGeometry>(STLGeometry::Load(ist));
+
+        mesh.SetDimension (3);
+
+        auto & points = geom->GetPoints();
+
+        for (auto & p : points)
+          mesh.AddPoint(MeshPoint(p));
+
+        mesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 1));
+
+        for (auto ti : IntRange(geom->GetNT()))
+        {
+          Element2d el(TRIG);
+          for (auto i : IntRange(3))
+            el[i] = (*geom)[STLTrigIndex(ti)][i];
+
+          el.SetIndex(1);
+
+          mesh.AddSurfaceElement(el);
+        }
+      }
   }
   
 }

libsrc/meshing/improve3.cpp

@@ -893,29 +893,127 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 
 
       
-bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
+  
-				    const NgBitArray * working_elements,
+
-  TABLE<ElementIndex,PointIndex::BASE> & elementsonnode,
+void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
-  INDEX_3_HASHTABLE<int> & faces,
+				    const NgBitArray * working_elements)
-  PointIndex pi1, PointIndex pi2, bool check_only
-                                    )
 {
+  static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t);
+  static Timer tloop("MeshOptimize3d::SwapImprove loop");
+  
   PointIndex pi3(PointIndex::INVALID), pi4(PointIndex::INVALID),
     pi5(PointIndex::INVALID), pi6(PointIndex::INVALID);
-
+  int cnt = 0;
-  double bad1, bad2, bad3;
 
   Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET);
   Element el1(TET), el2(TET), el3(TET), el4(TET);
   Element el1b(TET), el2b(TET), el3b(TET), el4b(TET);
-  ArrayMem<ElementIndex, 20> hasbothpoints;
   
-  bool do_swap = false;
+  double bad1, bad2, bad3;
+
+  int np = mesh.GetNP();
+  int ne = mesh.GetNE();
+  
+  // contains at least all elements at node
+  TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np);
+
+  NgArray<ElementIndex> hasbothpoints;
+
+  PrintMessage (3, "SwapImprove ");
+  (*testout) << "\n" << "Start SwapImprove" << endl;
+
+  const char * savetask = multithread.task;
+  multithread.task = "Swap Improve";
+  
+  //  mesh.CalcSurfacesOfNode ();
+    
+  INDEX_3_HASHTABLE<int> faces(mesh.GetNOpenElements()/3 + 2);
+  if (goal == OPT_CONFORM)
+    {
+      for (int i = 1; i <= mesh.GetNOpenElements(); i++)
+	{
+	  const Element2d & hel = mesh.OpenElement(i);
+	  INDEX_3 face(hel[0], hel[1], hel[2]);
+	  face.Sort();
+	  faces.Set (face, 1);
+	}
+    }
+  
+  // Calculate total badness
+  if (goal == OPT_QUALITY)
+    {
+      bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
+      (*testout) << "Total badness = " << bad1 << endl;
+    }
+  
+  // find elements on node
+  for (ElementIndex ei = 0; ei < ne; ei++)
+    for (PointIndex pi : mesh[ei].PNums())
+      elementsonnode.Add (pi, ei);
+    /*
+    for (int j = 0; j < mesh[ei].GetNP(); j++)
+      elementsonnode.Add (mesh[ei][j], ei);
+    */
+
+
+  // INDEX_2_HASHTABLE<int> edgeused(2 * ne + 5);
+  INDEX_2_CLOSED_HASHTABLE<int> edgeused(12 * ne + 5);
+
+  tloop.Start();
+  for (ElementIndex ei = 0; ei < ne; ei++)
+    {
+      if (multithread.terminate)
+	break;
+      
+      if (mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
+      	continue;
+      
+      multithread.percent = 100.0 * (ei+1) / ne;
+
+      if ((mesh.ElementType(ei)) == FIXEDELEMENT)
+	continue;
+
+      if(working_elements && 
+	 ei < working_elements->Size() &&
+	 !working_elements->Test(ei))
+	continue;
+
+      if (mesh[ei].IsDeleted())
+	continue;
+
+      if ((goal == OPT_LEGAL) && 
+	  mesh.LegalTet (mesh[ei]) &&
+	  CalcBad (mesh.Points(), mesh[ei], 0) < 1e3)
+	continue;
+
+      //      int onlybedges = 1;
+
+      for (int j = 0; j < 6; j++)
+	{
+	  // loop over edges
+
+	  const Element & elemi = mesh[ei];
+	  if (elemi.IsDeleted()) continue;
+
+	  int mattyp = elemi.GetIndex();
+	  
+	  static const int tetedges[6][2] =
+	    { { 0, 1 }, { 0, 2 }, { 0, 3 },
+	      { 1, 2 }, { 1, 3 }, { 2, 3 } };
+
+	  PointIndex pi1 = elemi[tetedges[j][0]];
+	  PointIndex pi2 = elemi[tetedges[j][1]];
+
 	  if (pi2 < pi1) Swap (pi1, pi2);
 	 	  
-	  if (mesh.BoundaryEdge (pi1, pi2)) return false;
+	  if (mesh.BoundaryEdge (pi1, pi2)) continue;
 
 
+	  INDEX_2 i2 (pi1, pi2);
+	  i2.Sort();
+	  if (edgeused.Used(i2)) continue;
+	  edgeused.Set (i2, 1);
+	  
 	  hasbothpoints.SetSize (0);
 	  for (int k = 0; k < elementsonnode[pi1].Size(); k++)
 	    {
@@ -923,7 +1021,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 	      ElementIndex elnr = elementsonnode[pi1][k];
 	      const Element & elem = mesh[elnr];
 	      
-	      if (elem.IsDeleted()) return false;
+	      if (elem.IsDeleted()) continue;
 	      
 	      for (int l = 0; l < elem.GetNP(); l++)
 		{
@@ -943,34 +1041,14 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 		}
 	    }
 	  
+	  bool puretet = true;
 	  for (ElementIndex ei : hasbothpoints)
-          {
 	    if (mesh[ei].GetType () != TET)
-              return false;
+	      puretet = false;
           
-            if (mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
+	  if (!puretet) continue;
-              return false;
-
-
-            if ((mesh.ElementType(ei)) == FIXEDELEMENT)
-              return false;
-
-            if(working_elements &&
-                ei < working_elements->Size() &&
-                !working_elements->Test(ei))
-              return false;
-
-            if (mesh[ei].IsDeleted())
-              return false;
-
-            if ((goal == OPT_LEGAL) &&
-                mesh.LegalTet (mesh[ei]) &&
-                CalcBad (mesh.Points(), mesh[ei], 0) < 1e3)
-              return false;
-          }
 
 	  int nsuround = hasbothpoints.Size();
-          int mattyp = mesh[hasbothpoints[0]].GetIndex();
 
 	  if ( nsuround == 3 )
 	    {
@@ -1027,6 +1105,9 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 	      el33[3] = pi3;  
 	      el33.SetIndex (mattyp);
 	      
+	      elementsonnode.Add (pi4, hasbothpoints[1]);
+	      elementsonnode.Add (pi3, hasbothpoints[2]);
+	      
 	      bad1 = CalcBad (mesh.Points(), el31, 0) + 
 		CalcBad (mesh.Points(), el32, 0) +
 		CalcBad (mesh.Points(), el33, 0);
@@ -1096,8 +1177,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 		{
 		  //		  (*mycout) << "3->2 " << flush;
 		  //		  (*testout) << "3->2 conversion" << endl;
-		  do_swap = true;
+		  cnt++;
-                  if(check_only) return do_swap;
 		  
 
 		  /*
@@ -1118,9 +1198,6 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 		    mesh[hasbothpoints[2]][l].Invalidate();
 		  mesh[hasbothpoints[2]].Delete();
 
-                  elementsonnode.Add (pi4, hasbothpoints[1]);
-                  elementsonnode.Add (pi3, hasbothpoints[2]);
-
 		  for (int k = 0; k < 2; k++)
 		    for (int l = 0; l < 4; l++)
 		      elementsonnode.Add (mesh[hasbothpoints[k]][l], hasbothpoints[k]);
@@ -1328,8 +1405,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 	      if (swap2 || swap3)
 		{
 		  // (*mycout) << "4->4 " << flush;
-		  do_swap = true;
+		  cnt++;
-                  if(check_only) return do_swap;
 		  //		  (*testout) << "4->4 conversion" << "\n";
 		  /*
 		    (*testout) << "bad1 = " << bad1 
@@ -1591,8 +1667,7 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 	      if (bestl != -1)
 		{
 		  // (*mycout) << nsuround << "->" << 2 * (nsuround-2) << " " << flush;
-		  do_swap = true;
+		  cnt++;
-                  if(check_only) return do_swap;
 		  
 		  for (int k = bestl+1; k <= nsuround + bestl - 2; k++)
 		    {
@@ -1642,105 +1717,53 @@ bool MeshOptimize3d :: SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal,
 		    }
 		}
 	    }
-          return do_swap;
 	}
 
-void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
+      /*
-				    const NgBitArray * working_elements)
+	if (onlybedges)
 	{
-  static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t);
+	(*testout) << "bad tet: " 
-  static Timer tloop("MeshOptimize3d::SwapImprove loop");
+	<< volelements.Get(i)[0] 
+	<< volelements.Get(i)[1] 
+	<< volelements.Get(i)[2] 
+	<< volelements.Get(i)[3] << "\n";
 
-  int cnt = 0;
+	if (!mesh.LegalTet (volelements.Get(i)))
-
+	cerr << "Illegal tet" << "\n";
-  int np = mesh.GetNP();
-  int ne = mesh.GetNE();
-
-  mesh.BoundaryEdge (1,2); // ensure the boundary-elements table is built
-
-  // contains at least all elements at node
-  TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np);
-
-  NgArray<ElementIndex> hasbothpoints;
-
-  PrintMessage (3, "SwapImprove ");
-  (*testout) << "\n" << "Start SwapImprove" << endl;
-
-  const char * savetask = multithread.task;
-  multithread.task = "Swap Improve";
-
-  INDEX_3_HASHTABLE<int> faces(mesh.GetNOpenElements()/3 + 2);
-  if (goal == OPT_CONFORM)
-    {
-      for (int i = 1; i <= mesh.GetNOpenElements(); i++)
-	{
-	  const Element2d & hel = mesh.OpenElement(i);
-	  INDEX_3 face(hel[0], hel[1], hel[2]);
-	  face.Sort();
-	  faces.Set (face, 1);
 	}
+      */
     }
+  //  (*mycout) << endl;
+  tloop.Stop();
+  /*  
+      cout << "edgeused: ";
+      edgeused.PrintMemInfo(cout);
+  */
+  PrintMessage (5, cnt, " swaps performed");
 
-  // Calculate total badness
+
+
+
+
+  mesh.Compress ();
+
+  /*
   if (goal == OPT_QUALITY)
     {
-      double bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
+      bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
-      (*testout) << "Total badness = " << bad1 << endl;
+      //      (*testout) << "Total badness = " << bad1 << endl;
-    }
-
-  // find elements on node
-  for (ElementIndex ei = 0; ei < ne; ei++)
-    for (PointIndex pi : mesh[ei].PNums())
-      elementsonnode.Add (pi, ei);
-
-  Array<std::tuple<PointIndex,PointIndex>> edges;
-  BuildEdgeList(mesh, elementsonnode, edges);
-
-  Array<int> candidate_edges(edges.Size());
-  std::atomic<int> improvement_counter(0);
-
-  tloop.Start();
-
-  ParallelForRange(Range(edges), [&] (auto myrange)
-  {
-    for(auto i : myrange)
-    {
-      if (multithread.terminate)
-        break;
-
-      auto [pi0, pi1] = edges[i];
-      if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true))
-        candidate_edges[improvement_counter++] = i;
-    }
-  });
-  // Sequential version:
-  /*
-  for(auto i : edges.Range())
-  {
-    if (multithread.terminate)
-      break;
-
-    auto [pi0, pi1] = edges[i];
-    if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, true))
-      candidate_edges[improvement_counter++] = i;
     }
   */
 
-  auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
+  /*
-  QuickSort(edges_with_improvement);
+    for (i = 1; i <= GetNE(); i++)
-
+    if (volelements.Get(i)[0])
-  for(auto ei : edges_with_improvement)
+    if (!mesh.LegalTet (volelements.Get(i)))
     {
-      auto [pi0,pi1] = edges[ei];
+    cout << "detected illegal tet, 2" << endl;
-        if(SwapImproveEdge (mesh, goal, working_elements, elementsonnode, faces, pi0, pi1, false))
+    (*testout) << "detected illegal tet1: " << i << endl;
-          cnt++;
     }
-
+  */
-  tloop.Stop();
-
-  PrintMessage (5, cnt, " swaps performed");
-
-  mesh.Compress ();
 
   multithread.task = savetask;
 }

libsrc/meshing/improve3.hpp

@@ -26,8 +26,6 @@ public:
   void CombineImproveSequential (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
 
   void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
-
-  bool SwapImproveEdge (Mesh & mesh, OPTIMIZEGOAL goal, const NgBitArray * working_elements, TABLE<ElementIndex,PointIndex::BASE> & elementsonnode, INDEX_3_HASHTABLE<int> & faces, PointIndex pi1, PointIndex pi2, bool check_only=false);
   void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
 		    const NgBitArray * working_elements = NULL);
   void SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,

libsrc/stlgeom/stltopology.cpp

@@ -338,6 +338,18 @@ void STLTopology :: Save (const char* filename) const
 
 STLGeometry *  STLTopology ::Load (istream & ist)
 {
+  // Check if the file starts with "solid". If not, the file is binary
+  {
+    char buf[5+1];
+    FIOReadStringE(ist,buf,5);
+    if (strcmp(buf, "solid") != 0)
+    {
+      for (auto i : Range(5))
+        ist.unget();
+      return LoadBinary(ist);
+    }
+  }
+
   STLGeometry * geom = new STLGeometry();
 
   NgArray<STLReadTriangle> readtrigs;

ng/dialog.tcl

@@ -2228,7 +2228,7 @@ proc stloptionsdialog { } {
 }
 
 proc stldoctordialog { } {
-    Ng_STLDoctor 0 0
+    Ng_STLDoctor
     set wd .stldoctor_dlg
 
     if {[winfo exists .stldoctor_dlg] == 1} {
@@ -2642,17 +2642,17 @@ proc stldoctordialog { } {
 
     
     #tixControl $f.gtol -label "load-geometry tolerance factor" -integer false \
-	-variable stldoctor.geom_tol_fact \
+	# -variable stldoctor.geom_tol_fact \
-	-options {
+	# -options {
 	#    entry.width 8
 	#    label.width 30
 	#    label.anchor e
 	#}	
-    ttk::spinbox $f.gtol -from 1 -to 20 -textvariable stldoctor.geom_tol_fact -width 8
+    ttk::label $f.gtol_lbl -text "LoadSTL tolerance factor"
-    pack $f.gtol
+    ttk::spinbox $f.gtol -from 1e-15 -to 0.001 -textvariable stldoctor.geom_tol_fact -width 8
+    pack $f.gtol_lbl $f.gtol
 
     ttk::button $f.adap -text "Apply" -command {
-	.stldoctor_dlg.nb.advanced.gtol invoke
 	Ng_STLDoctor; 
     }
     pack $f.adap -expand yes

ng/menustat.tcl

@@ -231,6 +231,7 @@ loadmeshinifile;
 	    {"Universal format"  {.unv} }
 	    {"Olaf format"  {.emt} }
 	    {"TET format" {.tet} }
+	    {"STL format" {.stl .stlb} }
 	    {"Pro/ENGINEER neutral format" {.fnf} }
 	          }
 	set file [tk_getOpenFile -filetypes $types ]

ng/ngpkg.cpp

@@ -405,6 +405,7 @@ namespace netgen
     PrintMessage (2, mesh->GetNP(), " Points, ",
 		  mesh->GetNE(), " Elements.");
 
+    SetGlobalMesh (mesh);
     mesh->SetGlobalH (mparam.maxh);
     mesh->CalcLocalH(mparam.grading);
 
@@ -1372,6 +1373,12 @@ namespace netgen
 		  return 0;
 		}
 	    }
+          else if (mesh)
+            {
+              MeshVolume(mparam, *mesh);
+              OptimizeVolume(mparam, *mesh);
+              return 0;
+            }
           else // no ng_geometry
             {
               multithread.task = savetask;

ng/onetcl.cpp

@@ -861,6 +861,7 @@ const char * ngscript[] = {""
 ,"{\"Universal format\"  {.unv} }\n"
 ,"{\"Olaf format\"  {.emt} }\n"
 ,"{\"TET format\" {.tet} }\n"
+,"{\"STL format\" {.stl .stlb} }\n"
 ,"{\"Pro/ENGINEER neutral format\" {.fnf} }\n"
 ,"}\n"
 ,"set file [tk_getOpenFile -filetypes $types ]\n"
@@ -2803,7 +2804,7 @@ const char * ngscript[] = {""
 ,"}\n"
 ,"}\n"
 ,"proc stldoctordialog { } {\n"
-,"Ng_STLDoctor 0 0\n"
+,"Ng_STLDoctor\n"
 ,"set wd .stldoctor_dlg\n"
 ,"if {[winfo exists .stldoctor_dlg] == 1} {\n"
 ,"wm withdraw $wd\n"
@@ -3054,10 +3055,10 @@ const char * ngscript[] = {""
 ,"-variable stldoctor.conecheck \\\n"
 ,"-command {Ng_STLDoctor;}\n"
 ,"pack $f.sc.bu $f.sc.bu2\n"
-,"ttk::spinbox $f.gtol -from 1 -to 20 -textvariable stldoctor.geom_tol_fact -width 8\n"
+,"ttk::label $f.gtol_lbl -text \"LoadSTL tolerance factor\"\n"
-,"pack $f.gtol\n"
+,"ttk::spinbox $f.gtol -from 1e-15 -to 0.001 -textvariable stldoctor.geom_tol_fact -width 8\n"
+,"pack $f.gtol_lbl $f.gtol\n"
 ,"ttk::button $f.adap -text \"Apply\" -command {\n"
-,".stldoctor_dlg.nb.advanced.gtol invoke\n"
 ,"Ng_STLDoctor;\n"
 ,"}\n"
 ,"pack $f.adap -expand yes\n"