Merge remote-tracking branch 'origin/master' into parallel_meshing

libsrc/core/profiler.hpp

@@ -188,6 +188,8 @@ namespace ngcore
 
     Timer( const std::string & name, TTracing, TTiming ) : timernr(Init(name)) { }
 
+    [[deprecated ("Use Timer(name, NoTracing/NoTiming) instead")]] Timer( const std::string & name, int ) : timernr(Init(name)) {}
+
     void SetName (const std::string & name)
     {
       NgProfiler::SetName (timernr, name);
@@ -276,6 +278,25 @@ namespace ngcore
     void operator=(RegionTimer &&) = delete;
   };
 
+  class [[deprecated("Use RegionTimer instead (now thread safe)")]] ThreadRegionTimer
+  {
+    size_t nr;
+    size_t tid;
+  public:
+    /// start timer
+    ThreadRegionTimer (size_t _nr, size_t _tid) : nr(_nr), tid(_tid)
+    { NgProfiler::StartThreadTimer(nr, tid); }
+    /// stop timer
+    ~ThreadRegionTimer ()
+    { NgProfiler::StopThreadTimer(nr, tid); }
+
+    ThreadRegionTimer() = delete;
+    ThreadRegionTimer(ThreadRegionTimer &&) = delete;
+    ThreadRegionTimer(const ThreadRegionTimer &) = delete;
+    void operator=(const ThreadRegionTimer &) = delete;
+    void operator=(ThreadRegionTimer &&) = delete;
+  };
+
   class RegionTracer
     {
       int nr;

libsrc/csg/genmesh.cpp

@@ -10,15 +10,19 @@
 
 namespace netgen
 {
-  NgArray<SpecialPoint> specpoints;
+
-  static NgArray<MeshPoint> spoints;
+  NgArray<SpecialPoint> global_specpoints;  // for visualization
+  //static NgArray<MeshPoint> spoints;
   
 #define TCL_OK 0
 #define TCL_ERROR 1
 
 
 
-  static void FindPoints (CSGeometry & geom, Mesh & mesh)
+  static void FindPoints (CSGeometry & geom,
+                          NgArray<SpecialPoint> &  specpoints,
+                          NgArray<MeshPoint> & spoints,
+                          Mesh & mesh)
   {
     PrintMessage (1, "Start Findpoints");
 
@@ -49,6 +53,12 @@ namespace netgen
     PrintMessage (2, "Analyze spec points");
     spc.AnalyzeSpecialPoints (geom, spoints, specpoints);
 
+    {
+      static mutex mut;
+      lock_guard<mutex> guard(mut);
+      global_specpoints = specpoints;
+    }
+    
     PrintMessage (5, "done");
 
     (*testout) << specpoints.Size() << " special points:" << endl;
@@ -67,7 +77,10 @@ namespace netgen
 
 
 
-  static void FindEdges (CSGeometry & geom, Mesh & mesh, MeshingParameters & mparam,
+  static void FindEdges (CSGeometry & geom, Mesh & mesh,
+                         NgArray<SpecialPoint> &  specpoints,
+                         NgArray<MeshPoint> & spoints,
+                         MeshingParameters & mparam,
                          const bool setmeshsize = false)
   {
     EdgeCalculation ec (geom, specpoints, mparam);
@@ -669,6 +682,10 @@ namespace netgen
   int CSGGenerateMesh (CSGeometry & geom, 
 		       shared_ptr<Mesh> & mesh, MeshingParameters & mparam)
   {
+    NgArray<SpecialPoint> specpoints;
+    NgArray<MeshPoint> spoints;
+
+    
     if (mesh && mesh->GetNSE() &&
 	!geom.GetNSolids())
       {
@@ -705,7 +722,7 @@ namespace netgen
 	  }
 
 	spoints.SetSize(0);
-	FindPoints (geom, *mesh);
+	FindPoints (geom, specpoints, spoints, *mesh);
       
 	PrintMessage (5, "find points done");
 
@@ -723,7 +740,7 @@ namespace netgen
 
     if (mparam.perfstepsstart <= MESHCONST_MESHEDGES)
       {
-	FindEdges (geom, *mesh, mparam, true);
+	FindEdges (geom, *mesh, specpoints, spoints, mparam, true);
 	if (multithread.terminate) return TCL_OK;
 #ifdef LOG_STREAM      
 	(*logout) << "Edges meshed" << endl
@@ -740,16 +757,16 @@ namespace netgen
 	    mesh->CalcLocalH(mparam.grading);
 	    mesh->DeleteMesh();
 	    
-	    FindPoints (geom, *mesh);
+	    FindPoints (geom, specpoints, spoints, *mesh);
 	    if (multithread.terminate) return TCL_OK;
-	    FindEdges (geom, *mesh, mparam, true);
+	    FindEdges (geom, *mesh, specpoints, spoints, mparam, true);
 	    if (multithread.terminate) return TCL_OK;
 	    
 	    mesh->DeleteMesh();
 	  
-	    FindPoints (geom, *mesh);
+	    FindPoints (geom, specpoints, spoints, *mesh);
 	    if (multithread.terminate) return TCL_OK;
-	    FindEdges (geom, *mesh, mparam);
+	    FindEdges (geom, *mesh, specpoints, spoints, mparam);
 	    if (multithread.terminate) return TCL_OK;
 	  }
       }

libsrc/csg/specpoin.cpp

@@ -21,7 +21,7 @@
 
 namespace netgen
 {
-  NgArray<Box<3> > boxes;
+  NgArray<Box<3> > boxes; // for visualizaton
 
 
   void ProjectToEdge (const Surface * f1, const Surface * f2, Point<3> & hp);
@@ -64,7 +64,7 @@ namespace netgen
   }
 
 
-  static NgArray<int> numprim_hist;
+  // static NgArray<int> numprim_hist;
 
   SpecialPointCalculation :: SpecialPointCalculation ()
   {
@@ -75,8 +75,8 @@ namespace netgen
   CalcSpecialPoints (const CSGeometry & ageometry, 
 		     NgArray<MeshPoint> & apoints)
   {
-    static int timer = NgProfiler::CreateTimer ("CSG: find special points");
+    // static int timer = NgProfiler::CreateTimer ("CSG: find special points");
-    NgProfiler::RegionTimer reg (timer);
+    // NgProfiler::RegionTimer reg (timer);
 
 
     geometry = &ageometry;
@@ -100,8 +100,8 @@ namespace netgen
     box.CalcDiamCenter();
     PrintMessage (3, "main-solids: ", geometry->GetNTopLevelObjects());
 
-    numprim_hist.SetSize (geometry->GetNSurf()+1);
+    // numprim_hist.SetSize (geometry->GetNSurf()+1);
-    numprim_hist = 0;
+    // numprim_hist = 0;
 
     for (int i = 0; i < geometry->GetNTopLevelObjects(); i++)
       {
@@ -161,10 +161,12 @@ namespace netgen
 
     PrintMessage (3, "Found points ", apoints.Size());
 
+    /*
     for (int i = 0; i < boxesinlevel.Size(); i++)
       (*testout) << "level " << i << " has " 
 		 << boxesinlevel[i] << " boxes" << endl;
     (*testout) << "numprim_histogramm = " << endl << numprim_hist << endl;
+    */
   }
   
 
@@ -184,8 +186,8 @@ namespace netgen
 
     if (multithread.terminate)
       {
-	*testout << "boxes = " << boxes << endl;
+	// *testout << "boxes = " << boxes << endl;
-	*testout << "boxesinlevel = " << boxesinlevel << endl;
+	// *testout << "boxesinlevel = " << boxesinlevel << endl;
 	throw NgException ("Meshing stopped");
       }
 
@@ -215,11 +217,12 @@ namespace netgen
 
     // static int cntbox = 0;
     // cntbox++;
-
+    /*
     if (level <= boxesinlevel.Size())
       boxesinlevel.Elem(level)++;
     else
       boxesinlevel.Append (1);
+    */
     
     /*
       numprim = sol -> NumPrimitives();
@@ -233,7 +236,7 @@ namespace netgen
     (*testout) << "numprim = " << numprim << endl;
 #endif
 
-    numprim_hist[numprim]++;
+    // numprim_hist[numprim]++;
 
     Point<3> p = box.Center();
 

libsrc/csg/vscsg.cpp

@@ -18,7 +18,9 @@ namespace netgen
   /* *********************** Draw Geometry **************** */
 
   extern shared_ptr<Mesh> mesh;
-  extern NgArray<SpecialPoint> specpoints;
+  extern NgArray<SpecialPoint> global_specpoints;
+  NgArray<SpecialPoint> & specpoints = global_specpoints;
+  
   extern NgArray<Box<3> > boxes;
 
 

libsrc/general/optmem.cpp

@@ -28,6 +28,7 @@ namespace netgen
 
   BlockAllocator :: ~BlockAllocator ()
   {
+    lock_guard<mutex> guard(block_allocator_mutex);     
     // cout << "****************** delete BlockAllocator " << endl;
     for (int i = 0; i < bablocks.Size(); i++)
       delete [] bablocks[i];

libsrc/gprim/adtree.cpp

@@ -400,23 +400,23 @@ namespace netgen
 				   const float * bmax,
 				   NgArray<int> & pis) const
   {
-    static NgArray<ADTreeNode3*> stack(1000);
+    ArrayMem<ADTreeNode3*, 1000> stack(1000);
-    static NgArray<int> stackdir(1000);
+    ArrayMem<int, 1000> stackdir(1000);
     ADTreeNode3 * node;
     int dir, stacks;
 
-    stack.SetSize (1000);
+    // stack.SetSize (1000);
-    stackdir.SetSize(1000);
+    // stackdir.SetSize(1000);
     pis.SetSize(0);
 
-    stack.Elem(1) = root;
+    stack[0] = root;
-    stackdir.Elem(1) = 0;
+    stackdir[0] = 0;
-    stacks = 1;
+    stacks = 0;
 
-    while (stacks)
+    while (stacks >= 0)
       {
-	node = stack.Get(stacks);
+	node = stack[stacks];
-	dir = stackdir.Get(stacks); 
+	dir = stackdir[stacks]; 
 	stacks--;
 
 	if (node->pi != -1)
@@ -436,14 +436,14 @@ namespace netgen
 	if (node->left && bmin[dir] <= node->sep)
 	  {
 	    stacks++;
-	    stack.Elem(stacks) = node->left;
+	    stack[stacks] = node->left;
-	    stackdir.Elem(stacks) = ndir;
+	    stackdir[stacks] = ndir;
 	  }
 	if (node->right && bmax[dir] >= node->sep)
 	  {
 	    stacks++;
-	    stack.Elem(stacks) = node->right;
+	    stack[stacks] = node->right;
-	    stackdir.Elem(stacks) = ndir;
+	    stackdir[stacks] = ndir;
 	  }
       }
   }

libsrc/interface/nginterface.cpp

@@ -146,27 +146,37 @@ void Ng_LoadMesh (const char * filename, ngcore::NgMPI_Comm comm)
 
   if( id == 0) {
 
-    string fn(filename);
-    if (fn.substr (fn.length()-3, 3) == ".gz")
-      infile = new igzstream (filename);
-    else
-      infile = new ifstream (filename);
     mesh.reset (new Mesh());
     mesh->SetCommunicator(comm);
-    mesh -> Load(*infile);
-    SetGlobalMesh (mesh);
     
-    // make string from rest of file (for geometry info!)
+    string fn(filename);
-    // (this might be empty, in which case we take the global ng_geometry)
+    if (fn.substr (fn.length()-8, 8) == ".vol.bin")
-    stringstream geom_part;
+      {
-    geom_part << infile->rdbuf();
+        mesh -> Load(fn);
-    string geom_part_string = geom_part.str();
+        SetGlobalMesh (mesh);        
-    strs = geom_part_string.size();
+      }
-    // buf = new char[strs];
+    else
-    buf.SetSize(strs);
+      {
-    memcpy(&buf[0], geom_part_string.c_str(), strs*sizeof(char));
+        if (fn.substr (fn.length()-3, 3) == ".gz")
+          infile = new igzstream (filename);
+        else
+          infile = new ifstream (filename);
+        mesh -> Load(*infile);
+        SetGlobalMesh (mesh);
+        
+        // make string from rest of file (for geometry info!)
+        // (this might be empty, in which case we take the global ng_geometry)
+        stringstream geom_part;
+        geom_part << infile->rdbuf();
+        string geom_part_string = geom_part.str();
+        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/linalg/densemat.hpp

@@ -171,8 +171,14 @@ public:
   { height = h; data = adata; ownmem = false; }
   ///
   MatrixFixWidth (const MatrixFixWidth & m2)
-    : height(m2.height), data(m2.data), ownmem(false)
+    : height(m2.height), ownmem(true)
-  { ; } 
+  {
+    data = new T[height*WIDTH];
+    for (int i = 0; i < WIDTH*height; i++)
+      data[i] = m2.data[i];
+  }
+  // : height(m2.height), data(m2.data), ownmem(false)
+  //{ ; } 
   ///
   ~MatrixFixWidth ()
   { if (ownmem) delete [] data; }
@@ -277,6 +283,15 @@ public:
   ///
   MatrixFixWidth (int h)
   { height = h; data = new double[WIDTH*height]; ownmem = true; }
+  
+  MatrixFixWidth (const MatrixFixWidth & m2)
+    : height(m2.height), ownmem(true)
+  {
+    data = new double[height*WIDTH];
+    for (int i = 0; i < WIDTH*height; i++)
+      data[i] = m2.data[i];
+  }
+
   /// 
   MatrixFixWidth (int h, double * adata)
   { height = h; data = adata; ownmem = false; }

libsrc/meshing/curvedelems.cpp

@@ -1936,7 +1936,7 @@ namespace netgen
 
   template <typename T>
   void CurvedElements :: 
-  CalcElementDShapes (SurfaceElementInfo & info, const Point<2,T> xi, MatrixFixWidth<2,T> dshapes) const
+  CalcElementDShapes (SurfaceElementInfo & info, const Point<2,T> xi, MatrixFixWidth<2,T> & dshapes) const
   {
     const Element2d & el = mesh[info.elnr];
     ELEMENT_TYPE type = el.GetType();
@@ -2981,7 +2981,7 @@ namespace netgen
 
   template <typename T>
   void CurvedElements :: 
-  CalcElementDShapes (ElementInfo & info, const Point<3,T> xi, MatrixFixWidth<3,T> dshapes) const
+  CalcElementDShapes (ElementInfo & info, const Point<3,T> xi, MatrixFixWidth<3,T> & dshapes) const
   {
     // static int timer = NgProfiler::CreateTimer ("calcelementdshapes");
     

libsrc/meshing/curvedelems.hpp

@@ -212,7 +212,7 @@ private:
   void CalcElementShapes (ElementInfo & info, Point<3,T> xi, TFlatVector<T> shapes) const;
   void GetCoefficients (ElementInfo & info, Vec<3> * coefs) const;
   template <typename T>  
-  void CalcElementDShapes (ElementInfo & info, const Point<3,T> xi, MatrixFixWidth<3,T> dshapes) const;
+  void CalcElementDShapes (ElementInfo & info, const Point<3,T> xi, MatrixFixWidth<3,T> & dshapes) const;
 
   template <typename T>
   bool EvaluateMapping (ElementInfo & info, const Point<3,T> xi, Point<3,T> & x, Mat<3,3,T> & jac) const;  
@@ -233,7 +233,7 @@ private:
   template <int DIM_SPACE>
   void GetCoefficients (SurfaceElementInfo & elinfo, NgArray<Vec<DIM_SPACE> > & coefs) const;
   template <typename T>
-  void CalcElementDShapes (SurfaceElementInfo & elinfo, const Point<2,T> xi, MatrixFixWidth<2,T> dshapes) const;
+  void CalcElementDShapes (SurfaceElementInfo & elinfo, const Point<2,T> xi, MatrixFixWidth<2,T> & dshapes) const;
 
   template <int DIM_SPACE, typename T>
   bool EvaluateMapping (SurfaceElementInfo & info, const Point<2,T> xi, Point<DIM_SPACE,T> & x, Mat<DIM_SPACE,2,T> & jac) const;  

libsrc/meshing/meshing2.cpp

@@ -46,17 +46,28 @@ namespace netgen
     static Timer t("Mesing2::Meshing2"); RegionTimer r(t);
 
     auto & globalrules = mp.quad ? global_quad_rules : global_trig_rules;
-    if (!globalrules.Size())
+
-      {
+    {
-        LoadRules (NULL, mp.quad);
+      static mutex mut;
-        for (auto * rule : rules)
+      lock_guard<mutex> guard(mut);
-          globalrules.Append (unique_ptr<netrule>(rule));
+      if (!globalrules.Size())
-      }
+        {
-    else
+          LoadRules (NULL, mp.quad);
-      {
+          for (auto & rule : rules)
-        for (auto i : globalrules.Range())
+            globalrules.Append (make_unique<netrule>(*rule));
-          rules.Append (globalrules[i].get());
+          rules.SetSize(0);
-      }
+        }
+      /*
+      else
+        {
+          for (auto i : globalrules.Range())
+            rules.Append (globalrules[i].get());
+        }
+      */
+    }
+    for (auto i : globalrules.Range())
+      rules.Append (make_unique<netrule>(*globalrules[i]));
+
     // LoadRules ("rules/quad.rls");
     // LoadRules ("rules/triangle.rls");
 
@@ -453,7 +464,7 @@ namespace netgen
 	      {
 		(*testout) << foundmap.Get(i) << "/" 
 			   << canuse.Get(i) << "/"
-			   << ruleused.Get(i) << " map/can/use rule " << rules.Get(i)->Name() << "\n";
+			   << ruleused.Get(i) << " map/can/use rule " << rules[i-1]->Name() << "\n";
 	      }
 	    (*testout) << "\n";
 	  }
@@ -1470,7 +1481,7 @@ namespace netgen
 	    if ( debugparam.haltsuccess || debugflag )
 	      {
 		// adfront.PrintOpenSegments (*testout);
-		cout << "success of rule" << rules.Get(rulenr)->Name() << endl;
+		cout << "success of rule" << rules[rulenr-1]->Name() << endl;
 		multithread.drawing = 1;
 		multithread.testmode = 1;
 		multithread.pause = 1;
@@ -1486,7 +1497,7 @@ namespace netgen
 		  }
 		*/
 
-		(*testout) << "success of rule" << rules.Get(rulenr)->Name() << endl;
+		(*testout) << "success of rule" << rules[rulenr-1]->Name() << endl;
 		(*testout) << "trials = " << trials << endl;
 
 		(*testout) << "locpoints " << endl;

libsrc/meshing/meshing2.hpp

@@ -31,7 +31,7 @@ class Meshing2
   /// the current advancing front
   AdFront2 adfront;
   /// rules for mesh generation
-  NgArray<netrule*> rules;
+  Array<unique_ptr<netrule>> rules;
   /// statistics
   NgArray<int> ruleused, canuse, foundmap;
   /// 

libsrc/meshing/netrule2.cpp

@@ -6,18 +6,20 @@ namespace netgen
 
 netrule :: netrule ()
 {
-  name = new char[1];
+  // name = new char[1];
-  name[0] = char(0);
+  // name[0] = char(0);
   quality = 0;
 }
 
 netrule ::  ~netrule()
 {
-  delete [] name;
+  // delete [] name;
+  /*
   for(int i = 0; i < oldutofreearea_i.Size(); i++)
     delete oldutofreearea_i[i];
   for(int i = 0; i < freezone_i.Size(); i++)
     delete freezone_i[i];
+  */
 }
 
 
@@ -37,9 +39,9 @@ void netrule :: SetFreeZoneTransformation (const Vector & devp, int tolclass)
 
   if (tolclass <= oldutofreearea_i.Size())
     {
-      oldutofreearea_i[tolclass-1] -> Mult (devp, devfree);
+      oldutofreearea_i[tolclass-1].Mult (devp, devfree);
 
-      auto& fzi = *freezone_i[tolclass-1];
+      auto& fzi = freezone_i[tolclass-1];
       for (int i = 0; i < fzs; i++)
 	{
 	  transfreezone[i][0] = fzi[i][0] + devfree[2*i];

libsrc/meshing/parser2.cpp

@@ -61,9 +61,12 @@ void netrule :: LoadRule (istream & ist)
   ist.get (ch);
 
   // if(name != NULL)
+  /*
   delete [] name;
   name = new char[strlen (buf) + 1];
   strcpy (name, buf);
+  */
+  name = string(buf);
   //(*testout) << "name " << name << endl;
   //  (*mycout) << "Rule " << name << " found." << endl;
 
@@ -474,14 +477,14 @@ void netrule :: LoadRule (istream & ist)
     {
       double lam1 = 1.0/(i+1);
 
-      oldutofreearea_i[i] = new DenseMatrix (oldutofreearea.Height(), oldutofreearea.Width());
+      oldutofreearea_i[i] = move(DenseMatrix (oldutofreearea.Height(), oldutofreearea.Width()));
-      DenseMatrix & mati = *oldutofreearea_i[i];
+      DenseMatrix & mati = oldutofreearea_i[i];
       for (j = 0; j < oldutofreearea.Height(); j++)
 	for (int k = 0; k < oldutofreearea.Width(); k++)
 	  mati(j,k) = lam1 * oldutofreearea(j,k) + (1 - lam1) * oldutofreearealimit(j,k);
 
-      freezone_i[i] = new NgArray<Point<2>> (freezone.Size());
+      freezone_i[i] = NgArray<Point<2>> (freezone.Size());
-      auto& fzi = *freezone_i[i];
+      auto& fzi = freezone_i[i];
       for (int j = 0; j < freezone.Size(); j++)
 	fzi[j] = freezonelimit[j] + lam1 * (freezone[j] - freezonelimit[j]);
     }
@@ -589,12 +592,12 @@ void Meshing2 :: LoadRules (const char * filename, bool quad)
       if (strcmp (buf, "rule") == 0)
 	{
 	  //(*testout) << "found rule" << endl;
-	  netrule * rule = new netrule;
+	  auto rule = make_unique<netrule>();
 	  //(*testout) << "fr1" << endl;
 	  rule -> LoadRule(*ist);
 	  //(*testout) << "fr2" << endl;
 	  
-	  rules.Append (rule);
+	  rules.Append (move(rule));
 	}
       //(*testout) << "loop" << endl;
     }

libsrc/meshing/ruler2.cpp

@@ -209,7 +209,7 @@ namespace netgen
     for (int ri = 1; ri <= rules.Size(); ri++)
       {
 	// NgProfiler::RegionTimer reg(timers[ri-1]);
-	netrule * rule = rules.Get(ri);
+	netrule * rule = rules[ri-1].get();
 
 #ifdef LOCDEBUG
 	if (loctestmode)

libsrc/meshing/ruler2.hpp

@@ -21,7 +21,7 @@ private:
   ///
   int quality;
   ///
-  char * name;
+  string name;
   ///
   NgArray<Point<2>> points;
   ///
@@ -29,7 +29,7 @@ private:
   ///
   NgArray<Point<2>> freezone, freezonelimit;
   ///
-  NgArray<NgArray<Point<2>>*> freezone_i;
+  NgArray<NgArray<Point<2>>> freezone_i;
   ///
   NgArray<Point<2>> transfreezone;
 
@@ -44,7 +44,7 @@ private:
   ///
   DenseMatrix oldutonewu, oldutofreearea, oldutofreearealimit;
   ///
-  NgArray<DenseMatrix*> oldutofreearea_i;
+  NgArray<DenseMatrix> oldutofreearea_i;
   ///
   MatrixFixWidth<3> freesetinequ;
 
@@ -154,7 +154,7 @@ public:
   ///
   const DenseMatrix & GetOldUToFreeArea () const { return oldutofreearea; }
   ///
-  const char * Name () const { return name; }
+  const string & Name () const { return name; }
 
   ///
   void LoadRule (istream & ist);