Merge branch 'parallel_surf_meshing' into 'master'

Parallel surf meshing

See merge request jschoeberl/netgen!395

libsrc/csg/genmesh.cpp

@@ -517,7 +517,7 @@ namespace netgen
 	  }
 
 	if (multithread.terminate) return;
-      
+        
 	for (SurfaceElementIndex sei = oldnf; sei < mesh.GetNSE(); sei++)
 	  mesh[sei].SetIndex (k);
 

libsrc/linalg/densemat.cpp

@@ -1380,6 +1380,6 @@ namespace netgen
     return ost;
   }
 
-
+  
 
 }

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

@@ -53,15 +53,21 @@ namespace netgen
       if (!globalrules.Size())
         {
           LoadRules (NULL, mp.quad);
-          for (auto * rule : rules)
+          for (auto & rule : rules)
-            globalrules.Append (unique_ptr<netrule>(rule));
+            globalrules.Append (make_unique<netrule>(*rule));
+          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");
 
@@ -458,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";
 	  }
@@ -1475,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;
@@ -1491,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

@@ -60,10 +60,13 @@ void netrule :: LoadRule (istream & ist)
   ist.get (buf, sizeof(buf), '"');
   ist.get (ch);
 
-  // if(name != NULL) 
+  // 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);