Merge branch 'stl_manifold' into 'master'

stl manifold meshing

See merge request ngsolve/netgen!516

libsrc/stlgeom/meshstlsurface.cpp

@@ -106,7 +106,7 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
 		   << ", trig2 = " << trig2
 		   << ", trig2b = " << trig2b << endl;
 
-	  if (trig1 <= 0 || trig2 <= 0 || trig1b <= 0 || trig2b <= 0)
+	  if (trig1 <= 0 || trig2 < 0 || trig1b <= 0 || trig2b < 0)
 	    {
 	      cout << "negative trigs, "
 		   << ", trig1 = " << trig1
@@ -177,10 +177,13 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
 	  mesh.AddSegment (seg);
 
 
+          if(trig2 != 0)
+            {
 	  Segment seg2;
 	  seg2[0] = p2 + PointIndex::BASE-1;;
 	  seg2[1] = p1 + PointIndex::BASE-1;;
 	  seg2.si = geom.GetTriangle(trig2).GetFaceNum();
+
 	  seg2.edgenr = i;
 
 	  seg2.epgeominfo[0].edgenr = i;
@@ -219,8 +222,8 @@ static void STLFindEdges (STLGeometry & geom, Mesh & mesh,
 	      (*testout) << "Get GeomInfo PROBLEM" << endl;
 	    }
 	  */	  
-
 	  mesh.AddSegment (seg2);
+            }
 	}
     }
 

libsrc/stlgeom/python_stl.cpp

@@ -9,7 +9,7 @@
 using namespace netgen;
 namespace netgen
 {
-  //extern shared_ptr<Mesh> mesh;
+  extern shared_ptr<Mesh> mesh;
   extern shared_ptr<NetgenGeometry> ng_geometry;
 }
 
@@ -125,11 +125,12 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
 {
   py::class_<STLGeometry,shared_ptr<STLGeometry>, NetgenGeometry> (m,"STLGeometry")
     .def(py::init<>())
-    .def(py::init<>([](const string& filename)
+    .def(py::init<>([](const string& filename, bool surface)
                     {
                       ifstream ist(filename);
-                      return shared_ptr<STLGeometry>(STLGeometry::Load(ist));
+                      return shared_ptr<STLGeometry>(STLGeometry::Load(ist,
-                    }), py::arg("filename"),
+                                                                       surface));
+                    }), py::arg("filename"), py::arg("surface")=false,
       py::call_guard<py::gil_scoped_release>())
     .def(NGSPickle<STLGeometry>())
     .def("_visualizationData", [](shared_ptr<STLGeometry> stl_geo)
@@ -203,7 +204,10 @@ NGCORE_API_EXPORT void ExportSTL(py::module & m)
                            SetGlobalMesh(mesh);
                            auto result = STLMeshingDummy(geo.get(), mesh, mp, stlparam);
                            if(result != 0)
-                             throw Exception("Meshing failed!");
+                             {
+                               netgen::mesh = mesh;
+                               throw Exception("Meshing failed!");
+                             }
 
                            return mesh;
                          }, py::arg("mp") = nullptr,

libsrc/stlgeom/stlgeom.cpp

@@ -2578,7 +2578,7 @@ void STLGeometry :: CalcEdgeDataAngles()
   for (int i = 1; i <= GetNTE(); i++)
     {
       STLTopEdge & edge = GetTopEdge (i);
-      double cosang = 
+      double cosang = edge.TrigNum(2) == 0 ? 1. :
 	GetTriangle(edge.TrigNum(1)).Normal() *
 	GetTriangle(edge.TrigNum(2)).Normal();
       edge.SetCosAngle (cosang);
@@ -2611,6 +2611,8 @@ void STLGeometry :: FindEdgesFromAngles(const STLParameters& stlparam)
   for (int i = 1; i <= edgedata->Size(); i++)
     {
       STLTopEdge & sed = edgedata->Elem(i);
+      if(sed.TrigNum(2) == 0)
+        sed.SetStatus(ED_CONFIRMED);
       if (sed.GetStatus() == ED_CANDIDATE || 
 	  sed.GetStatus() == ED_UNDEFINED)
 	{
@@ -3187,7 +3189,7 @@ void STLGeometry :: BuildSmoothEdges ()
       ng1 = trig.GeomNormal(points);
       ng1 /= (ng1.Length() + 1e-24);
 
-      for (int j = 1; j <= 3; j++)
+      for (int j = 1; j <= NONeighbourTrigs(i); j++)
 	{ 
 	  int nbt = NeighbourTrig (i, j);
 
@@ -3261,7 +3263,7 @@ void STLGeometry :: AddConeAndSpiralEdges(const STLParameters& stlparam)
 	  STLTrigId t = chart.GetChartTrig1(j); 
 	  const STLTriangle& tt = GetTriangle(t);
 
-	  for (int k = 1; k <= 3; k++)
+	  for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	    {
 	      STLTrigId nt = NeighbourTrig(t,k); 
 	      if (GetChartNr(nt) != i && !TrigIsInOC(nt,i))

libsrc/stlgeom/stlgeomchart.cpp

@@ -238,7 +238,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
 			      */
 			      
 			      //find overlapping charts exacter (fast, too): 
-			      for (int k = 1; k <= 3; k++) 
+			      for (int k = 1; k <= NONeighbourTrigs(nt); k++)
 				{ 
 				  int nnt = NeighbourTrig(nt,k);
 				  if (GetMarker(nnt) != chartnum)
@@ -387,7 +387,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh, const MeshingParameters& mparam, cons
                       // NgProfiler::StartTimer (timer4a);
                       
 		      if (spiralcheckon && !isdirtytrig)
-			for (int k = 1; k <= 3; k++) 
+			for (int k = 1; k <= NONeighbourTrigs(nt); k++)
 			  {
                             // NgProfiler::StartTimer (timer4b);                            
 			    STLTrigId nnt = NeighbourTrig(nt,k);
@@ -695,7 +695,7 @@ void STLGeometry :: GetInnerChartLimes(NgArray<twoint>& limes, ChartId chartnum)
     {
       STLTrigId t = chart.GetChartTrig1(j); 
       const STLTriangle& tt = GetTriangle(t);
-      for (int k = 1; k <= 3; k++)
+      for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	{
 	  STLTrigId nt = NeighbourTrig(t,k); 
 	  if (GetChartNr(nt) != chartnum)
@@ -756,7 +756,7 @@ void STLGeometry :: GetDirtyChartTrigs(int chartnum, STLChart& chart,
       STLTrigId t = chart.GetChartTrig1(j); 
       const STLTriangle& tt = GetTriangle(t);
       
-      for (int k = 1; k <= 3; k++)
+      for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	{
 	  STLTrigId nt = NeighbourTrig(t,k); 
 	  if (GetChartNr(nt) != chartnum && outercharttrigs[nt] != chartnum)

libsrc/stlgeom/stlgeommesh.cpp

@@ -1169,7 +1169,7 @@ void STLGeometry :: RestrictHChartDistOneChart(ChartId chartnum, NgArray<int>& a
     {
       int t = chart.GetChartTrig1(j); 
       tt = GetTriangle(t);
-      for (int k = 1; k <= 3; k++)
+      for (int k = 1; k <= NONeighbourTrigs(t); k++)
 	{
 	  int nt = NeighbourTrig(t,k); 
 	  if (GetChartNr(nt) != chartnum)
@@ -1495,6 +1495,9 @@ int STLMeshingDummy (STLGeometry* stlgeometry, shared_ptr<Mesh> & mesh, const Me
       if (multithread.terminate)
 	return 0;
 
+      if(stlgeometry->IsSurfaceSTL())
+        return 0;
+
       if (mparam.perfstepsstart <= MESHCONST_MESHVOLUME && 
 	  mparam.perfstepsend >= MESHCONST_MESHVOLUME)
 	{

libsrc/stlgeom/stltopology.cpp

@@ -338,7 +338,7 @@ void STLTopology :: Save (const filesystem::path & filename) const
 }
 
 
-STLGeometry *  STLTopology ::Load (istream & ist)
+STLGeometry *  STLTopology ::Load (istream & ist, bool surface)
 {
   // Check if the file starts with "solid". If not, the file is binary
   {
@@ -457,6 +457,7 @@ STLGeometry *  STLTopology ::Load (istream & ist)
       PrintWarning("File has normal vectors which differ extremely from geometry->correct with stldoctor!!!");
     }
 
+  geom->surface = surface;
   geom->InitSTLGeometry(readtrigs);
   return geom;
 }
@@ -650,6 +651,7 @@ void STLTopology :: FindNeighbourTrigs()
 	  }
       }
 
+  if(!surface)
   for (int i = 1; i <= ne; i++)
     {
       const STLTopEdge & edge = GetTopEdge (i);
@@ -668,9 +670,12 @@ void STLTopology :: FindNeighbourTrigs()
 	  const STLTriangle & t = GetTriangle (i);
 	  for (int j = 1; j <= 3; j++)
 	    {
-	      const STLTriangle & nbt = GetTriangle (t.NBTrigNum(j));
+              if(t.NBTrigNum(j) != 0)
-	      if (!t.IsNeighbourFrom (nbt))
+                {
-		orientation_ok = 0;
+                  const STLTriangle & nbt = GetTriangle (t.NBTrigNum(j));
+                  if (!t.IsNeighbourFrom (nbt))
+                    orientation_ok = 0;
+                }
 	    }
 	}
     }
@@ -801,7 +806,8 @@ void STLTopology :: FindNeighbourTrigs()
       neighbourtrigs.SetSize(GetNT());
       for (int i = 1; i <= GetNT(); i++)
 	for (int k = 1; k <= 3; k++)
-	  AddNeighbourTrig (i, GetTriangle(i).NBTrigNum(k));
+          if(GetTriangle(i).NBTrigNum(k) != 0)
+            AddNeighbourTrig (i, GetTriangle(i).NBTrigNum(k));
     }
   else
     {

libsrc/stlgeom/stltopology.hpp

@@ -120,7 +120,12 @@ public:
 
 
   STLTriangle (const STLPointId * apts);
-  STLTriangle () {pts[0]=0;pts[1]=0;pts[2]=0;}
+  STLTriangle ()
+  {
+    pts[0]=0;pts[1]=0;pts[2]=0;
+    nbtrigs[0][0] = nbtrigs[0][1] = nbtrigs[0][2] = 0.;
+    nbtrigs[1][0] = nbtrigs[1][1] = nbtrigs[1][2] = 0.;
+  }
 
   void DoArchive(Archive& ar)
   {
@@ -282,6 +287,7 @@ protected:
   Array<STLTriangle, STLTrigId> trias;
   NgArray<STLTopEdge> topedges;
   Array<Point<3>, STLPointId> points;
+  bool surface = false;
 
   // mapping of sorted pair of points to topedge
   INDEX_2_HASHTABLE<int> * ht_topedges;
@@ -313,13 +319,15 @@ public:
   virtual ~STLTopology();
 
   static STLGeometry * LoadNaomi (istream & ist);
-  static STLGeometry * Load (istream & ist);
+  static STLGeometry * Load (istream & ist, bool surface=false);
   static STLGeometry * LoadBinary (istream & ist);
 
   void Save (const filesystem::path & filename) const;
   void SaveBinary (const filesystem::path & filename, const char* aname) const;
   void SaveSTLE (const filesystem::path & filename) const; // stores trigs and edges
 
+  bool IsSurfaceSTL() const { return surface; }
+
   virtual void DoArchive(Archive& ar)
   {
     ar & trias & points & boundingbox & pointtol;