modernize Delaunay

libsrc/meshing/delaunay.cpp

@@ -6,7 +6,6 @@
 namespace netgen
 {
 
-
   static const int deltetfaces[][3] = 
     { { 1, 2, 3 },
       { 2, 0, 3 },
@@ -14,10 +13,6 @@ namespace netgen
       { 1, 0, 2 } };
 
 
-
-
-
-
   class DelaunayTet
   {
     PointIndex pnums[4];
@@ -41,9 +36,6 @@ namespace netgen
     PointIndex & operator[] (int i) { return pnums[i]; }
     PointIndex operator[] (int i) const { return pnums[i]; }
 
-    // int & NB1(int i) { return nb[i-1]; }
-    // int NB1(int i) const { return nb[i-1]; }
-
     int & NB(int i) { return nb[i]; }
     int NB(int i) const { return nb[i]; }
 
@@ -58,29 +50,6 @@ namespace netgen
       return 3;
     }
 
-    /*
-    void GetFace1 (int i, INDEX_3 & face) const
-    {
-      face.I(1) = pnums[deltetfaces[i-1][0]];
-      face.I(2) = pnums[deltetfaces[i-1][1]];
-      face.I(3) = pnums[deltetfaces[i-1][2]];
-    }
-    */
-    void GetFace (int i, INDEX_3 & face) const
-    {
-      face.I(1) = pnums[deltetfaces[i][0]];
-      face.I(2) = pnums[deltetfaces[i][1]];
-      face.I(3) = pnums[deltetfaces[i][2]];
-    }
-
-    /*
-    INDEX_3 GetFace1 (int i) const
-    {
-      return INDEX_3 (pnums[deltetfaces[i-1][0]],
-		      pnums[deltetfaces[i-1][1]],
-		      pnums[deltetfaces[i-1][2]]);
-    }
-    */
     INDEX_3 GetFace (int i) const
     {
       return INDEX_3 (pnums[deltetfaces[i][0]],
@@ -88,15 +57,6 @@ namespace netgen
 		      pnums[deltetfaces[i][2]]);
     }
 
-    /*
-    void GetFace1 (int i, Element2d & face) const
-    {
-      // face.SetType(TRIG);
-      face[0] = pnums[deltetfaces[i-1][0]];
-      face[1] = pnums[deltetfaces[i-1][1]];
-      face[2] = pnums[deltetfaces[i-1][2]];
-    }
-    */
     void GetFace (int i, Element2d & face) const
     {
       // face.SetType(TRIG);
@@ -104,7 +64,6 @@ namespace netgen
       face[1] = pnums[deltetfaces[i][1]];
       face[2] = pnums[deltetfaces[i][2]];
     }
-
   };
 
 
@@ -113,8 +72,6 @@ namespace netgen
 
 
 
-
-
   /*
     Table to maintain neighbour elements
   */
@@ -147,7 +104,7 @@ namespace netgen
     // get neighbour of element elnr in direction fnr 
     int GetNB (int elnr, int fnr)
     { 
-      return tets.Get(elnr).NB(fnr-1); 
+      return tets.Get(elnr).NB(fnr); 
     }
 
     //
@@ -193,7 +150,6 @@ namespace netgen
 
 
 
-
   /*
     connected lists of cosphereical elements
   */
@@ -265,24 +221,23 @@ namespace netgen
 			 Array<int> & freelist, SphereList & list,
 			 IndexSet & insphere, IndexSet & closesphere)
   {
+    // static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t);
     /*
       find any sphere, such that newp is contained in
     */
-  
     DelaunayTet el;
     int cfelind = -1;
 
     const Point<3> * pp[4];
     Point<3> pc;
-    double r2;
     Point3d tpmin, tpmax;
 
     tettree.GetIntersecting (newp, newp, treesearch);
+
     double quot,minquot(1e20);
 
-    for (int j = 0; j < treesearch.Size(); j++)
+    for (auto jjj : treesearch)
       {
-	int jjj = treesearch[j];
 	quot = Dist2 (centers.Get(jjj), newp) / radi2.Get(jjj);
 	
 	if((cfelind == -1 || quot < 0.99*minquot) && quot < 1)
@@ -296,46 +251,12 @@ namespace netgen
       }
 
 
-    /*
-      int i, j, k, l;
-      if (!felind)
-      {
-      cerr << "not in any sphere, 1" << endl;
-      // old, non tree search
-
-      double mindist = 1e10;
-      for (j = 1; j <= tempels.Size(); j++)
-      {
-      if (tempels.Get(j).PNum(1))
-      {
-      double toofar = 
-      Dist2 (centers.Get(j), newp) - radi2.Get(j);
-      if (toofar < mindist || toofar < 1e-7) 
-      {
-      mindist = toofar;
-      cout << " dist2 = " << Dist2 (centers.Get(j), newp)
-      << " radi2 = " << radi2.Get(j) << endl;
-      }
-      if (toofar < 0)
-      {
-      el = tempels.Get(j);
-      felind = j;
-      cout << "sphere found !" << endl;
-      break; 
-      }
-      }
-      }
-      cout << "point is too far from sheres: " << mindist << endl;
-      }
-    */      
-
     if (cfelind == -1)
       {
 	PrintWarning ("Delaunay, point not in any sphere");
 	return;
       }
 	
-
     /*
       insphere:     point is in sphere -> delete element
       closesphere:  point is close to sphere -> considered for same center
@@ -388,7 +309,7 @@ namespace netgen
 	
 	// check neighbour-tets
 	for (int j = starti; j < nstarti; j++)
-	  for (int k = 1; k <= 4; k++)
+	  for (int k = 0; k < 4; k++)
 	    {
 	      int helind = insphere.GetArray().Get(j);
 	      int nbind = meshnb.GetNB (helind, k);
@@ -416,61 +337,48 @@ namespace netgen
 		    }
 		  else
 		    {
-		      /*
+		      INDEX_3 i3 = tempels.Get(helind).GetFace (k);
-		      Element2d face;
-		      tempels.Get(helind).GetFace (k, face);
 
-		      const Point3d & p1 = mesh.Point (face.PNum(1));
+		      const Point<3> & p1 = mesh.Point ( PointIndex (i3.I1()) );
-		      const Point3d & p2 = mesh.Point (face[1]);
+		      const Point<3> & p2 = mesh.Point ( PointIndex (i3.I2()) );
-		      const Point3d & p3 = mesh.Point (face[2]);
+		      const Point<3> & p3 = mesh.Point ( PointIndex (i3.I3()) );
-		      */
 
-		      INDEX_3 i3 = tempels.Get(helind).GetFace (k-1);
+		      Vec<3> v1 = p2-p1;
-
+		      Vec<3> v2 = p3-p1;
-		      const Point3d & p1 = mesh.Point ( PointIndex (i3.I1()) );
+		      Vec<3> n = Cross (v1, v2);
-		      const Point3d & p2 = mesh.Point ( PointIndex (i3.I2()) );
-		      const Point3d & p3 = mesh.Point ( PointIndex (i3.I3()) );
-
-
-		      Vec3d v1(p1, p2);
-		      Vec3d v2(p1, p3);
-		      Vec3d n = Cross (v1, v2);
 		      n /= n.Length();
 
-		      if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k-1])) > 0)
+		      if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k])) > 0)
 			n *= -1;
 			
 		      double dist = n * Vec3d (p1, newp);
 
-
 		      if (dist > -1e-10)  // 1e-10
 			{
 			  insphere.Add (nbind);
 			  changed = 1;
 			}
-
-
 		    }
 		}
 	    }
       } // while (changed)
-    //      (*testout) << "newels: " << endl;
+
-    Array<Element> newels;
+    // Array<Element> newels;
+    Array<DelaunayTet> newels;
 
     Element2d face(TRIG);
 
-    for (int j = 1; j <= insphere.GetArray().Size(); j++)
+    for (int celind : insphere.GetArray())
-      for (int k = 1; k <= 4; k++)
+      for (int k = 0; k < 4; k++)
 	{
-	  //	    int elind = insphere.GetArray().Get(j);
-	  int celind = insphere.GetArray().Get(j);
 	  int nbind = meshnb.GetNB (celind, k);
 
 	  if (!nbind || !insphere.IsIn (nbind))
 	    {
-	      tempels.Get (celind).GetFace (k-1, face);
+	      tempels.Get (celind).GetFace (k, face);
 		
-	      Element newel(TET);
+	      // Element newel(TET);
+              DelaunayTet newel;
 	      for (int l = 0; l < 3; l++)
                 newel[l] = face[l];
               newel[3] = newpi;
@@ -483,7 +391,7 @@ namespace netgen
 
 	      n.Normalize();
 	      if (n * Vec3d(mesh.Point (face[0]), 
-			    mesh.Point (tempels.Get(insphere.GetArray().Get(j))[k-1]))
+			    mesh.Point (tempels.Get(celind)[k]))
 		  > 0)
 		n *= -1;
 
@@ -504,38 +412,36 @@ namespace netgen
 
     meshnb.ResetFaceHT (10*insphere.GetArray().Size()+1);
 
-    for (int j = 1; j <= insphere.GetArray().Size(); j++)
+    for (auto celind : insphere.GetArray())
       {
-	//	  int elind = 
-	int celind = insphere.GetArray().Get(j);
-
 	meshnb.Delete (celind); 
 	list.DeleteElement (celind);
 	  
 	for (int k = 0; k < 4; k++)
 	  tempels.Elem(celind)[k] = -1;
         
-	// ((ADTree6&)tettree.Tree()).DeleteElement (celind);
         tettree.DeleteElement (celind);
 	freelist.Append (celind);
       }
 
-    int hasclose = 0;
+
-    for (int j = 1; j <= closesphere.GetArray().Size(); j++)
+    
+    bool hasclose = false;
+    for (int ind : closesphere.GetArray())
       {
-	int ind = closesphere.GetArray().Get(j);
 	if (!insphere.IsIn(ind) &&
 	    fabs (Dist2 (centers.Get (ind), newp) - radi2.Get(ind)) < 1e-8 )
-	  hasclose = 1;
+	  hasclose = true;
       }
 
     for (int j = 1; j <= newels.Size(); j++)
       {
+        const auto & newel = newels.Get(j);
 	int nelind;
 
 	if (!freelist.Size())
 	  {
-	    tempels.Append (newels.Get(j));
+	    tempels.Append (newel);
 	    nelind = tempels.Size();
 	  }
 	else
@@ -543,29 +449,29 @@ namespace netgen
 	    nelind = freelist.Last();
 	    freelist.DeleteLast();
 
-	    tempels.Elem(nelind) = newels.Get(j);
+	    tempels.Elem(nelind) = newel;
 	  }
 
 	meshnb.Add (nelind);
 	list.AddElement (nelind);
 
 	for (int k = 0; k < 4; k++)
-	  pp[k] = &mesh.Point (newels.Get(j)[k]);
+	  pp[k] = &mesh.Point (newel[k]);
 
 	if (CalcSphereCenter (&pp[0], pc) )
 	  {
 	    PrintSysError ("Delaunay: New tet is flat");
 
 	    (*testout) << "new tet is flat" << endl;
-	    for (int k = 1; k <= 4; k++)
+	    for (int k = 0; k < 4; k++)
-	      (*testout) << newels.Get(j).PNum(k) << " ";
+	      (*testout) << newel[k] << " ";
 	    (*testout) << endl;
-	    for (int k = 1; k <= 4; k++)
+	    for (int k = 0; k < 4; k++)
 	      (*testout) << *pp[k-1] << " ";
 	    (*testout) << endl;
 	  }
 
-	r2 = Dist2 (*pp[0], pc);
+	double r2 = Dist2 (*pp[0], pc);
 	if (hasclose)
 	  for (int k = 1; k <= closesphere.GetArray().Size(); k++)
 	    {
@@ -615,33 +521,11 @@ namespace netgen
 		  int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax)
   {
     static Timer t("Meshing3::Delaunay1"); RegionTimer reg(t);
+    static Timer tloop("Meshing3::Delaunay1 loop");
     
     Array<Point<3>> centers;
     Array<double> radi2;
   
-    Point3d tpmin, tpmax;
-
-
-    // new: local box
-    /*
-    mesh.GetBox (pmax, pmin);   // lower bound for pmax, upper for pmin
-    for (int i = 1; i <= adfront->GetNF(); i++)
-      {
-	const MiniElement2d & face = adfront->GetFace(i);
-        for (PointIndex pi : face.PNums())
-	  {
-	    pmin.SetToMin  (mesh.Point (pi));
-	    pmax.SetToMax  (mesh.Point (pi));
-	  }
-      }
-  
-    for (PointIndex pi : mesh.LockedPoints())
-      {
-	pmin.SetToMin (mesh.Point (pi));
-	pmax.SetToMax (mesh.Point (pi));
-      }
-    cout << "pmin = " << pmin << ", pmax = " << pmax << endl;
-    */
     Box<3> bbox(Box<3>::EMPTY_BOX);
 
     for (auto & face : adfront->Faces())
@@ -653,26 +537,24 @@ namespace netgen
     
     pmin = bbox.PMin();
     pmax = bbox.PMax();
-    // cout << "bbox = " << bbox << endl;
 
-    Vec3d vdiag(pmin, pmax);
+
+    Vec<3> vdiag = pmax-pmin;
     // double r1 = vdiag.Length();
-    double r1 = sqrt (3.0) * max3(vdiag.X(), vdiag.Y(), vdiag.Z());
+    double r1 = sqrt (3.0) * max3(vdiag(0), vdiag(1), vdiag(2));
-    vdiag = Vec3d (r1, r1, r1);
+    vdiag = Vec<3> (r1, r1, r1);
     //double r2;
 
-    Point3d pmin2 = pmin - 8 * vdiag;
+    Point<3> pmin2 = pmin - 8 * vdiag;
-    Point3d pmax2 = pmax + 8 * vdiag;
+    Point<3> pmax2 = pmax + 8 * vdiag;
 
-    Point3d cp1(pmin2), cp2(pmax2), cp3(pmax2), cp4(pmax2);
+    Point<3> cp1(pmin2), cp2(pmax2), cp3(pmax2), cp4(pmax2);
-    cp2.X() = pmin2.X();
+    cp2(0) = pmin2(0);
-    cp3.Y() = pmin2.Y();
+    cp3(1) = pmin2(1);
-    cp4.Z() = pmin2.Z();
+    cp4(2) = pmin2(2);
 
 
-
+    size_t np = mesh.GetNP();
-
-    int np = mesh.GetNP();
 
     startel[0] = mesh.AddPoint (cp1);
     startel[1] = mesh.AddPoint (cp2);
@@ -682,24 +564,21 @@ namespace netgen
     // flag points to use for Delaunay:
     BitArrayChar<PointIndex::BASE> usep(np);
     usep.Clear();
-    for (int i = 1; i <= adfront->GetNF(); i++)
-      {
-	const MiniElement2d & face = adfront->GetFace(i);
-	for (int j = 0; j < face.GetNP(); j++)
-	  usep.Set (face[j]);
-      }
 
-    for (int i = oldnp + PointIndex::BASE; 
+    for (auto & face : adfront->Faces())
+      for (PointIndex pi : face.Face().PNums())      
+        usep.Set (pi);
+    
+    for (size_t i = oldnp + PointIndex::BASE; 
 	 i < np + PointIndex::BASE; i++)
       usep.Set (i);
 
-    for (int i = 0; i < mesh.LockedPoints().Size(); i++)
+    for (PointIndex pi : mesh.LockedPoints())
-      usep.Set (mesh.LockedPoints()[i]);
+      usep.Set (pi);
     
 
     Array<int> freelist;
 
-
     int cntp = 0;
 
     MeshNB meshnb (tempels, mesh.GetNP() + 5);
@@ -716,13 +595,12 @@ namespace netgen
     list.AddElement (1);
     Array<int> connected, treesearch;
 
+    Box<3> tbox(Box<3>::EMPTY_BOX);
+    for (size_t k = 0; k < 4; k++)
+      tbox.Add (mesh.Point(startel[k]));
+    Point<3> tpmin = tbox.PMin();
+    Point<3> tpmax = tbox.PMax();
     
-    tpmin = tpmax = mesh.Point(startel[0]);
-    for (int k = 1; k < 4; k++)
-      {
-	tpmin.SetToMin (mesh.Point (startel[k]));
-	tpmax.SetToMax (mesh.Point (startel[k]));
-      }
     tpmax = tpmax + 0.01 * (tpmax - tpmin);
     tettree.Insert (tpmin, tpmax, 1);
 
@@ -754,6 +632,7 @@ namespace netgen
     for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++)
       mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE );
 
+    tloop.Start();
     for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++)
       {
 	if (pi % 1000 == 0)
@@ -782,6 +661,7 @@ namespace netgen
 			  connected, treesearch, freelist, list, insphere, closesphere);
 
       }
+    tloop.Stop();
     
     for (int i = tempels.Size(); i >= 1; i--)
       if (tempels.Get(i)[0] <= 0)
@@ -839,26 +719,27 @@ namespace netgen
       // improve delaunay - mesh by swapping !!!!
 
       Mesh tempmesh;
-      for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
-	tempmesh.AddPoint (mesh[pi]);
 
-      for (int i = 1; i <= tempels.Size(); i++)
+      for (auto & meshpoint : mesh.Points())
+        tempmesh.AddPoint (meshpoint);
+      
+      for (auto & tempel : tempels)
 	{   
 	  Element el(4);
 	  for (int j = 0; j < 4; j++)
-	    el[j] = tempels.Elem(i)[j];
+            el[j] = tempel[j];
 
 	  el.SetIndex (1);
 
-	  const Point3d & lp1 = mesh.Point (el[0]);
+	  const Point<3> & lp1 = mesh.Point (el[0]);
-	  const Point3d & lp2 = mesh.Point (el[1]);
+	  const Point<3> & lp2 = mesh.Point (el[1]);
-	  const Point3d & lp3 = mesh.Point (el[2]);
+	  const Point<3> & lp3 = mesh.Point (el[2]);
-	  const Point3d & lp4 = mesh.Point (el[3]);
+	  const Point<3> & lp4 = mesh.Point (el[3]);
-	  Vec3d v1(lp1, lp2);
+	  Vec<3> v1 = lp2-lp1;
-	  Vec3d v2(lp1, lp3);
+	  Vec<3> v2 = lp3-lp1;
-	  Vec3d v3(lp1, lp4);
+	  Vec<3> v3 = lp4-lp1;
 
-	  Vec3d n = Cross (v1, v2);
+	  Vec<3> n = Cross (v1, v2);
 	  double vol = n * v3;
 	  if (vol > 0) swap (el[2], el[3]);
 
@@ -1261,6 +1142,7 @@ namespace netgen
 
   
     INDEX_3_HASHTABLE<int> boundaryfaces(mesh.GetNOpenElements()/3+1);
+    /*
     for (int i = 1; i <= mesh.GetNOpenElements(); i++)
       {
 	const Element2d & tri = mesh.OpenElement(i);
@@ -1268,6 +1150,13 @@ namespace netgen
 	i3.Sort();
 	boundaryfaces.PrepareSet (i3);
       }
+    */
+    for (const Element2d & tri : mesh.OpenElements())
+      {
+	INDEX_3 i3 (tri[0], tri[1], tri[2]);
+	i3.Sort();
+	boundaryfaces.PrepareSet (i3);
+      }
     boundaryfaces.AllocateElements();
     for (int i = 1; i <= mesh.GetNOpenElements(); i++)
       {
@@ -1277,14 +1166,23 @@ namespace netgen
 	boundaryfaces.Set (i3, 1);
       }
 
+    /*
     for (int i = 0; i < tempels.Size(); i++)
       for (int j = 0; j < 4; j++)
 	tempels[i].NB(j) = 0;
+    */
+    for (auto & el : tempels)
+      for (int j = 0; j < 4; j++)
+	el.NB(j) = 0;
       
     TABLE<int,PointIndex::BASE> elsonpoint(mesh.GetNP());
+    /*
     for (int i = 0; i < tempels.Size(); i++)
       {
 	const DelaunayTet & el = tempels[i];
+    */
+    for (const DelaunayTet & el : tempels)
+      {
 	INDEX_4 i4(el[0], el[1], el[2], el[3]);
 	i4.Sort();
 	elsonpoint.IncSizePrepare (i4.I1());
@@ -1308,7 +1206,8 @@ namespace netgen
     INDEX_3_CLOSED_HASHTABLE<INDEX_2> faceht(100);   
   
     Element2d hel(TRIG);
-    for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
+    // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++)
+    for (PointIndex pi : mesh.Points().Range())
       {
 	faceht.SetSize (4 * elsonpoint[pi].Size());
 	for (int ii = 0; ii < elsonpoint[pi].Size(); ii++)

libsrc/meshing/meshclass.hpp

@@ -454,6 +454,7 @@ namespace netgen
     const Element2d & OpenElement(int i) const
     { return openelements.Get(i); }
 
+    auto & OpenElements() const { return openelements; }
     
     /// are also quads open elements
     bool HasOpenQuads () const;