Merge branch 'vol_bin_format' into 'master'

.vol.bin format using archiver

See merge request jschoeberl/netgen!385

libsrc/general/ngarray.hpp

@@ -242,7 +242,7 @@ namespace netgen
     using NgFlatArray<T,BASE,TIND>::data;
 
     /// physical size of array
-    size_t allocsize;
+    size_t allocsize = 0;
     /// memory is responsibility of container
     bool ownmem;
 

libsrc/include/nginterface_v2.hpp

@@ -358,8 +358,8 @@ namespace netgen
     int GetParentSElement (int ei) const;
 
     bool HasParentEdges() const;
-    tuple<int, std::array<int,3>> GetParentEdges (int enr) const;
+    std::tuple<int, std::array<int,3>> GetParentEdges (int enr) const;
-    tuple<int, std::array<int,4>> GetParentFaces (int fnr) const;
+    std::tuple<int, std::array<int,4>> GetParentFaces (int fnr) const;
     
     int GetNIdentifications() const;
     int GetIdentificationType(int idnr) const;

libsrc/meshing/delaunay.cpp

@@ -864,6 +864,8 @@ namespace netgen
 
 
     // remove degenerated
+    static Timer tdegenerated("Delaunay - remove degenerated");
+    tdegenerated.Start();
 
     NgBitArray badnode(mesh.GetNP());
     badnode.Clear();
@@ -914,6 +916,10 @@ namespace netgen
 
   
     PrintMessage (3, ndeg, " degenerated elements removed");
+    tdegenerated.Stop();
+
+    static Timer topenel("Delaunay - find openel");
+    topenel.Start();
 
     // find surface triangles which are no face of any tet
 
@@ -1073,6 +1079,10 @@ namespace netgen
       }
 
 
+    topenel.Stop();
+
+    static Timer trem_intersect("Delaunay - remove intersecting");
+    trem_intersect.Start();
 
 
     // find intersecting:
@@ -1190,6 +1200,10 @@ namespace netgen
       }
   
 
+    trem_intersect.Stop();
+
+    static Timer trem_outer("Delaunay - remove outer");
+    trem_outer.Start();
 
 
     PrintMessage (3, "Remove outer");
@@ -1400,15 +1414,17 @@ namespace netgen
       }
     */
 
+    int lowest_undefined_el = 1;
     while (1)
       {
 	int inside;
 	bool done = 1;
 
 	int i;
-	for (i = 1; i <= ne; i++)
+	for (i = lowest_undefined_el; i <= ne; i++)
 	  if (!inner.Test(i) && !outer.Test(i))
 	    {
+              lowest_undefined_el = i+1;
 	      done = 0;
 	      break;
 	    }
@@ -1643,6 +1659,8 @@ namespace netgen
 
     PrintMessage (5, "outer removed");
 
+    trem_outer.Stop();
+
     mesh.FindOpenElements(domainnr);
 
     mesh.Compress();

libsrc/meshing/improve2.hpp

@@ -10,7 +10,7 @@ void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & element
     { { 0, 1 }, { 0, 2 }, { 0, 3 },
         { 1, 2 }, { 1, 3 }, { 2, 3 } };
 
-  int ntasks = 2*ngcore::TaskManager::GetMaxThreads();
+  int ntasks = 4*ngcore::TaskManager::GetMaxThreads();
   Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);
 
   ParallelFor(IntRange(ntasks), [&] (int ti)

libsrc/meshing/improve3.cpp

@@ -2768,7 +2768,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
         candidate_edges[index] = make_tuple(d_badness, i);
       }
     }
-  });
+  }, TasksPerThread (4));
 
   auto edges_with_improvement = candidate_edges.Part(0, improvement_counter.load());
   QuickSort(edges_with_improvement);

libsrc/meshing/meshclass.cpp

@@ -413,6 +413,13 @@ namespace netgen
 
   void Mesh :: Save (const string & filename) const
   {
+    if (filename.find(".vol.bin") != string::npos)
+    {
+        BinaryOutArchive in(filename);
+        in & const_cast<Mesh&>(*this);
+        return;
+    }
+
     ostream * outfile;
     if (filename.find(".vol.gz")!=string::npos)
       outfile = new ogzstream(filename.c_str());
@@ -872,6 +879,14 @@ namespace netgen
   void Mesh :: Load (const string & filename)
   {
     cout << "filename = " << filename << endl;
+
+    if (filename.find(".vol.bin") != string::npos)
+    {
+        BinaryInArchive in(filename);
+        in & (*this);
+        return;
+    }
+
     istream * infile = NULL;
 
     if (filename.find(".vol.gz") != string::npos)
@@ -2100,38 +2115,28 @@ namespace netgen
     
     t_table.Start();
 
-    TableCreator<ElementIndex, PointIndex> creator(np);
+    auto elsonpoint = ngcore::CreateSortedTable<ElementIndex, PointIndex>( volelements.Range(),
+           [&](auto & table, ElementIndex ei)
+           {
+             const Element & el = (*this)[ei];
+             if (dom == 0 || dom == el.GetIndex())
+               {
+                 if (el.GetNP() == 4)
+                   {
+                     INDEX_4 i4(el[0], el[1], el[2], el[3]);
+                     i4.Sort();
+                     table.Add (PointIndex(i4.I1()), ei);
+                     table.Add (PointIndex(i4.I2()), ei);
+                   }
+                 else
+                   {
+                     for (PointIndex pi : el.PNums())
+                       table.Add(pi, ei);
+                   }
+               }
+           }, GetNP());
 
-    for ( ; !creator.Done(); creator++)
-      // for (ElementIndex ei : Range(VolumeElements()))
-      ParallelFor
-        (Range(VolumeElements()), [&] (ElementIndex ei)
-         {
-           const Element & el = (*this)[ei];
-           if (dom == 0 || dom == el.GetIndex())
-             {
-               if (el.GetNP() == 4)
-                 {
-                   INDEX_4 i4(el[0], el[1], el[2], el[3]);
-                   i4.Sort();
-                   creator.Add (PointIndex(i4.I1()), ei);
-                   creator.Add (PointIndex(i4.I2()), ei);
-                 }
-               else
-                 {
-                   for (PointIndex pi : el.PNums())
-                     creator.Add(pi, ei);
-                 }
-             }
-         });
-    
-    auto elsonpoint = creator.MoveTable();
 
-    ParallelFor (Range(elsonpoint), [&] (auto i)
-         {
-           QuickSort(elsonpoint[i]);
-         });
-           
     NgArray<int,PointIndex::BASE> numonpoint(np);
     /*
     numonpoint = 0;
@@ -2381,7 +2386,7 @@ namespace netgen
 
     */
 
-    size_t numtasks = ngcore::TaskManager::GetNumThreads();
+    size_t numtasks = 4*ngcore::TaskManager::GetNumThreads();
     Array<Array<Element2d>> thread_openelements(numtasks);
     ParallelJob
       ( [&](TaskInfo & ti)
@@ -2506,7 +2511,7 @@ namespace netgen
                         thread_openelements[ti.task_nr].Append (tri);
                       }
                   }
-            }});
+            }}, numtasks);
 
     for (auto & a : thread_openelements)
       for (auto & el : a)

libsrc/meshing/meshtype.hpp

@@ -446,7 +446,7 @@ namespace netgen
         pnum[cnt++] = val;
       return *this;
     }
-    Element2d & operator= (initializer_list<tuple<PointIndex,PointGeomInfo>> list)
+    Element2d & operator= (initializer_list<std::tuple<PointIndex,PointGeomInfo>> list)
     {
       size_t cnt = 0;
       for (auto val : list)

libsrc/meshing/python_mesh.cpp

@@ -685,29 +685,33 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
 		return;
 	      }
 
-	    istream * infile;
+	    istream * infile = nullptr;
 	    NgArray<char> buf; // for distributing geometry!
 	    int strs;
 
 	    if( id == 0) {
-
+	      if (filename.substr (filename.length()-8, 8) == ".vol.bin")
-	      if (filename.substr (filename.length()-3, 3) == ".gz")
+                mesh -> Load(filename);
+              else if (filename.substr (filename.length()-3, 3) == ".gz")
 		infile = new igzstream (filename.c_str());
 	      else
 		infile = new ifstream (filename.c_str());
-	      mesh -> Load(*infile);
 
-	      // make string from rest of file (for geometry info!)
+              if(infile)
-	      // (this might be empty, in which case we take the global ng_geometry)
+                {
-	      stringstream geom_part;
+                  mesh -> Load(*infile);
-	      geom_part << infile->rdbuf();
+                  // make string from rest of file (for geometry info!)
-	      string geom_part_string = geom_part.str();
+                  // (this might be empty, in which case we take the global ng_geometry)
-	      strs = geom_part_string.size();
+                  stringstream geom_part;
-	      // buf = new char[strs];
+                  geom_part << infile->rdbuf();
-	      buf.SetSize(strs);
+                  string geom_part_string = geom_part.str();
-	      memcpy(&buf[0], geom_part_string.c_str(), strs*sizeof(char));
+                  strs = geom_part_string.size();
+                  // buf = new char[strs];
+                  buf.SetSize(strs);
+                  memcpy(&buf[0], geom_part_string.c_str(), strs*sizeof(char));
+                  delete infile;
+                }
 
-	      delete infile;
 
 	      if (ntasks > 1)
 		{

libsrc/meshing/topology.hpp

@@ -195,10 +195,10 @@ public:
 
 
 private:
-  Array<tuple<int, std::array<int,3>>> parent_edges; 
+  Array<std::tuple<int, std::array<int,3>>> parent_edges;
   void BuildParentEdges ();
 
-  Array<tuple<int, std::array<int,4>>> parent_faces;
+  Array<std::tuple<int, std::array<int,4>>> parent_faces;
   void BuildParentFaces ();
 public:
   auto GetParentEdges (int enr) const { return parent_edges[enr]; }