Merge branch 'pointnd_to_point' into 'master'

change PointNd to Point<N> adfront2

See merge request jschoeberl/netgen!253

libsrc/gprim/geomobjects.hpp

@@ -162,8 +162,43 @@ namespace netgen
     Vec<D> GetNormal () const;
   };
 
+  template<int D>
+  inline Vec<D> operator-(const Point<D>& p1, const Point<D>& p2)
+  {
+    Vec<D> result;
+    for(auto i : Range(D))
+      result[i] = p1[i] - p2[i];
+    return result;
+  }
 
+  inline double Cross2(const Vec<2>& v1, const Vec<2>& v2)
+  {
+    return v1[0] * v2[1] - v1[1] * v2[0];
+  }
 
+  // are points clockwise?
+  inline bool CW(const Point<2>& p1, const Point<2>& p2,
+                  const Point<2>& p3)
+  {
+    return Cross2(p2-p1, p3-p2) < 0;
+  }
+
+  // are points counterclockwise?
+  inline bool CCW(const Point<2>& p1, const Point<2>& p2,
+                  const Point<2>& p3)
+  {
+    return Cross2(p2-p1, p3-p2) > 0;
+  }
+
+  // are strictly points counterclockwise?
+  inline bool CCW(const Point<2>& p1, const Point<2>& p2,
+                  const Point<2>& p3, double eps)
+  {
+    auto v1 = p2-p1;
+    auto v2 = p3-p2;
+    return Cross2(v1, v2) > eps*eps*max2(v1.Length2(),
+                                         v2.Length2());
+  }
 
 
   template <int H, int W=H, typename T = double>

libsrc/meshing/adfront2.cpp

@@ -274,7 +274,7 @@ namespace netgen
 
 
   int AdFront2 :: GetLocals (int baselineindex,
-			     NgArray<Point3d> & locpoints,
+			     NgArray<Point<3>> & locpoints,
 			     NgArray<MultiPointGeomInfo> & pgeominfo,
 			     NgArray<INDEX_2> & loclines,   // local index
 			     NgArray<INDEX> & pindex,

libsrc/meshing/adfront2.hpp

@@ -216,7 +216,7 @@ public:
 
   ///
   int GetLocals (int baseline, 
-		 NgArray<Point3d> & locpoints,
+		 NgArray<Point<3>> & locpoints,
 		 NgArray<MultiPointGeomInfo> & pgeominfo,
                  NgArray<INDEX_2> & loclines,   // local index
                  NgArray<int> & pindex,

libsrc/meshing/meshing2.cpp

@@ -246,11 +246,11 @@ namespace netgen
 
     // double h;
 
-    auto locpointsptr = make_shared<NgArray<Point3d>>();
+    auto locpointsptr = make_shared<NgArray<Point<3>>>();
     vssurfacemeshing.locpointsptr = locpointsptr;
     auto& locpoints = *locpointsptr;
     NgArray<int> legalpoints;
-    auto plainpointsptr = make_shared<NgArray<Point2d>>();
+    auto plainpointsptr = make_shared<NgArray<Point<2>>>();
     auto& plainpoints = *plainpointsptr;
     vssurfacemeshing.plainpointsptr = plainpointsptr;
     NgArray<int> plainzones;
@@ -698,7 +698,7 @@ namespace netgen
             
 	    double avy = 0;
 	    for (size_t i = 0; i < plainpoints.Size(); i++)
-	      avy += plainpoints[i].Y();
+	      avy += plainpoints[i][1];
 	    avy *= 1./plainpoints.Size();
 		
 
@@ -715,7 +715,7 @@ namespace netgen
 		  }
 		    
 
-		if (plainpoints[i].Y() < -1e-10*avy) // changed
+		if (plainpoints[i][1] < -1e-10*avy) // changed
 		  {
 		    legalpoints[i] = 0;
 		  }

libsrc/meshing/meshing2.hpp

@@ -135,7 +135,7 @@ protected:
 
 /** Applies 2D rules.
  Tests all 2D rules */
-  int ApplyRules (NgArray<Point2d> & lpoints, 
+  int ApplyRules (NgArray<Point<2>> & lpoints, 
 		  NgArray<int> & legalpoints,
 		  int maxlegalpoint,
 		  NgArray<INDEX_2> & llines,

libsrc/meshing/meshtype.cpp

@@ -439,7 +439,7 @@ namespace netgen
   }
 
   void Element2d :: 
-  GetTransformation (int ip, const NgArray<Point2d> & points,
+  GetTransformation (int ip, const NgArray<Point<2>> & points,
                      DenseMatrix & trans) const
   {
     int np = GetNP();
@@ -658,7 +658,7 @@ namespace netgen
 
 
   void Element2d :: 
-  GetPointMatrix (const NgArray<Point2d> & points,
+  GetPointMatrix (const NgArray<Point<2>> & points,
                   DenseMatrix & pmat) const
   {
     int np = GetNP();
@@ -673,9 +673,9 @@ namespace netgen
   
     for (int i = 1; i <= np; i++)
       {
-        const Point2d & p = points.Get(PNum(i));
+        const auto& p = points.Get(PNum(i));
-        pmat.Elem(1, i) = p.X();
+        pmat.Elem(1, i) = p[0];
-        pmat.Elem(2, i) = p.Y();
+        pmat.Elem(2, i) = p[1];
       }
   }
 
@@ -683,7 +683,7 @@ namespace netgen
 
 
 
-  double Element2d :: CalcJacobianBadness (const NgArray<Point2d> & points) const
+  double Element2d :: CalcJacobianBadness (const NgArray<Point<2>> & points) const
   {
     int i, j;
     int nip = GetNIP();
@@ -727,7 +727,7 @@ namespace netgen
     };
 
   double Element2d :: 
-  CalcJacobianBadnessDirDeriv (const NgArray<Point2d> & points,
+  CalcJacobianBadnessDirDeriv (const NgArray<Point<2>> & points,
                                int pi, Vec2d & dir, double & dd) const
   {
     if (typ == QUAD)

libsrc/meshing/meshtype.hpp

@@ -580,7 +580,7 @@ namespace netgen
     int GetNIP () const;
     void GetIntegrationPoint (int ip, Point2d & p, double & weight) const;
 
-    void GetTransformation (int ip, const NgArray<Point2d> & points,
+    void GetTransformation (int ip, const NgArray<Point<2>> & points,
 			    class DenseMatrix & trans) const;
     void GetTransformation (int ip, class DenseMatrix & pmat,
 			    class DenseMatrix & trans) const;
@@ -595,16 +595,16 @@ namespace netgen
     void GetDShapeNew (const Point<2,T> & p, class MatrixFixWidth<2,T> & dshape) const;
     
     /// matrix 2 * np
-    void GetPointMatrix (const NgArray<Point2d> & points,
+    void GetPointMatrix (const NgArray<Point<2>> & points,
 			 class DenseMatrix & pmat) const; 
 
     void ComputeIntegrationPointData () const;
   
 
-    double CalcJacobianBadness (const NgArray<Point2d> & points) const;
+    double CalcJacobianBadness (const NgArray<Point<2>> & points) const;
     double CalcJacobianBadness (const T_POINTS & points, 
 				const Vec<3> & n) const;
-    double CalcJacobianBadnessDirDeriv (const NgArray<Point2d> & points,
+    double CalcJacobianBadnessDirDeriv (const NgArray<Point<2>> & points,
 					int pi, Vec2d & dir, double & dd) const;
 
 

libsrc/meshing/netrule2.cpp

@@ -39,11 +39,11 @@ void netrule :: SetFreeZoneTransformation (const Vector & devp, int tolclass)
     {
       oldutofreearea_i[tolclass-1] -> Mult (devp, devfree);
 
-      NgArray<Point2d> & fzi = *freezone_i[tolclass-1];
+      auto& fzi = *freezone_i[tolclass-1];
       for (int i = 0; i < fzs; i++)
 	{
-	  transfreezone[i].X() = fzi[i].X() + devfree[2*i];
+	  transfreezone[i][0] = fzi[i][0] + devfree[2*i];
-	  transfreezone[i].Y() = fzi[i].Y() + devfree[2*i+1];
+	  transfreezone[i][1] = fzi[i][1] + devfree[2*i+1];
 	}
     }
   else
@@ -57,24 +57,24 @@ void netrule :: SetFreeZoneTransformation (const Vector & devp, int tolclass)
 
       for (int i = 0; i < fzs; i++)
 	{
-	  transfreezone[i].X() = lam1 * freezone[i].X() + lam2 * freezonelimit[i].X() + devfree[2*i];
+	  transfreezone[i][0] = lam1 * freezone[i][0] + lam2 * freezonelimit[i][0] + devfree[2*i];
-	  transfreezone[i].Y() = lam1 * freezone[i].Y() + lam2 * freezonelimit[i].Y() + devfree[2*i+1];
+	  transfreezone[i][1] = lam1 * freezone[i][1] + lam2 * freezonelimit[i][1] + devfree[2*i+1];
 	}
     }
 
 
   if (fzs > 0)
     {
-      fzmaxx = fzminx = transfreezone[0].X();
+      fzmaxx = fzminx = transfreezone[0][0];
-      fzmaxy = fzminy = transfreezone[0].Y();
+      fzmaxy = fzminy = transfreezone[0][1];
     }
 
   for (int i = 1; i < fzs; i++)
     {
-      if (transfreezone[i].X() > fzmaxx) fzmaxx = transfreezone[i].X();
+      if (transfreezone[i][0] > fzmaxx) fzmaxx = transfreezone[i][0];
-      if (transfreezone[i].X() < fzminx) fzminx = transfreezone[i].X();
+      if (transfreezone[i][0] < fzminx) fzminx = transfreezone[i][0];
-      if (transfreezone[i].Y() > fzmaxy) fzmaxy = transfreezone[i].Y();
+      if (transfreezone[i][1] > fzmaxy) fzmaxy = transfreezone[i][1];
-      if (transfreezone[i].Y() < fzminy) fzminy = transfreezone[i].Y();
+      if (transfreezone[i][1] < fzminy) fzminy = transfreezone[i][1];
     }
 
   for (int i = 0; i < fzs; i++)
@@ -147,8 +147,8 @@ int netrule :: IsLineInFreeZone2 (const Point2d & p1, const Point2d & p2) const
 
       for (int i = 1; i <= transfreezone.Size(); i++)
 	{
-	  bool left  = transfreezone.Get(i).X() * nx + transfreezone.Get(i).Y() * ny + c <  1e-7;
+	  bool left  = transfreezone.Get(i)[0] * nx + transfreezone.Get(i)[1] * ny + c <  1e-7;
-          bool right = transfreezone.Get(i).X() * nx + transfreezone.Get(i).Y() * ny + c > -1e-7;
+          bool right = transfreezone.Get(i)[0] * nx + transfreezone.Get(i)[1] * ny + c > -1e-7;
 	  if (!left) allleft = false;
 	  if (!right) allright = false;
 	}
@@ -187,33 +187,12 @@ float netrule :: CalcPointDist (int pi, const Point2d & p) const
 }
 */
 
-float netrule :: CalcLineError (int li, const Vec2d & v) const
+float netrule :: CalcLineError (int li, const Vec<2> & v) const
 {
-  float dx = v.X() - linevecs.Get(li).X();
+  float dx = v[0] - linevecs.Get(li)[0];
-  float dy = v.Y() - linevecs.Get(li).Y();
+  float dy = v[1] - linevecs.Get(li)[1];
 
   const threefloat * ltf = &linetolerances.Get(li);
   return ltf->f1 * dx * dx + ltf->f2 * dx * dy + ltf->f3 * dy * dy;
 }
-
+} // namespace netgen
-
-
-
-/*
-int GetNRules ()
-  {
-  return rules.Size();
-  }
-*/
-
-
-
-
-
-
-
-
-
-
-
-}

libsrc/meshing/parser2.cpp

@@ -480,8 +480,8 @@ void netrule :: LoadRule (istream & ist)
 	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<Point2d> (freezone.Size());
+      freezone_i[i] = new NgArray<Point<2>> (freezone.Size());
-      NgArray<Point2d> & 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]);
     }

libsrc/meshing/ruler2.cpp

@@ -5,13 +5,13 @@ namespace netgen
 {
 
 
-  static double CalcElementBadness (const NgArray<Point2d> & points,
+  static double CalcElementBadness (const NgArray<Point<2>> & points,
 				    const Element2d & elem)
   {
     // badness = sqrt(3) /36 * circumference^2 / area - 1 +
     //           h / li + li / h - 2
 
-    Vec2d v12, v13, v23;
+    Vec<2> v12, v13, v23;
     double l12, l13, l23, cir, area;
     static const double c = sqrt(3.0) / 36;
 
@@ -24,7 +24,7 @@ namespace netgen
     l23 = v23.Length();
 
     cir = l12 + l13 + l23;
-    area = 0.5 * (v12.X() * v13.Y() - v12.Y() * v13.X());
+    area = 0.5 * (v12[0] * v13[1] - v12[1] * v13[0]);
     if (area < 1e-6)
       {
 	return 1e8;
@@ -45,7 +45,7 @@ namespace netgen
 
 
 
-  int Meshing2 ::ApplyRules (NgArray<Point2d> & lpoints, 
+  int Meshing2 ::ApplyRules (NgArray<Point<2>> & lpoints, 
 			     NgArray<int> & legalpoints,
 			     int maxlegalpoint,
 			     NgArray<INDEX_2> & llines1,
@@ -69,7 +69,7 @@ namespace netgen
 
     NgArrayMem<int, 20> pmap, pfixed, lmap;
   
-    NgArrayMem<Point2d,100> tempnewpoints;
+    NgArrayMem<Point<2>,100> tempnewpoints;
     NgArrayMem<INDEX_2,100> tempnewlines;
     NgArrayMem<int,100> tempdellines;
     NgArrayMem<Element2d,100> tempelements;
@@ -263,7 +263,7 @@ namespace netgen
 		    ok = 1;
 
 		    INDEX_2 loclin = llines.Get(locli);
-		    Vec2d linevec = lpoints.Get(loclin.I2()) - lpoints.Get(loclin.I1());
+		    auto linevec = lpoints.Get(loclin.I2()) - lpoints.Get(loclin.I1());
 
 		    if (rule->CalcLineError (nlok, linevec) > maxerr)
 		      {

libsrc/meshing/ruler2.hpp

@@ -23,15 +23,15 @@ private:
   ///
   char * name;
   ///
-  NgArray<Point2d> points;
+  NgArray<Point<2>> points;
   ///
   NgArray<INDEX_2> lines;
   ///
-  NgArray<Point2d> freezone, freezonelimit;
+  NgArray<Point<2>> freezone, freezonelimit;
   ///
-  NgArray<NgArray<Point2d>*> freezone_i;
+  NgArray<NgArray<Point<2>>*> freezone_i;
   ///
-  NgArray<Point2d> transfreezone;
+  NgArray<Point<2>> transfreezone;
 
   ///
   NgArray<int> dellines;
@@ -49,7 +49,7 @@ private:
   MatrixFixWidth<3> freesetinequ;
 
   ///
-  NgArray<Vec2d> linevecs;
+  NgArray<Vec<2>> linevecs;
 
   ///
   int noldp, noldl;
@@ -86,7 +86,7 @@ public:
   int GetLNearness (int li) const { return lnearness.Get(li); }
 
   ///
-  const Point2d & GetPoint (int i) const { return points.Get(i); }
+  const Point<2>& GetPoint (int i) const { return points.Get(i); }
   ///
   const INDEX_2 & GetLine (int i) const { return lines.Get(i); }
   ///
@@ -103,14 +103,14 @@ public:
 
   double CalcPointDist (int pi, const Point2d & p) const
   {
-    double dx = p.X() - points.Get(pi).X();
+    double dx = p.X() - points.Get(pi)[0];
-    double dy = p.Y() - points.Get(pi).Y();
+    double dy = p.Y() - points.Get(pi)[1];
     const threefloat * tfp = &tolerances.Get(pi);
     return tfp->f1 * dx * dx + tfp->f2 * dx * dy + tfp->f3 * dy * dy;
   }
 
   ///
-  float CalcLineError (int li, const Vec2d & v) const;
+  float CalcLineError (int li, const Vec<2>& v) const;
 
   ///
   void SetFreeZoneTransformation (const Vector & u, int tolclass);
@@ -144,7 +144,7 @@ public:
   ///
   int ConvexFreeZone () const;
   ///
-  const NgArray<Point2d> & GetTransFreeZone () { return transfreezone; }
+  const NgArray<Point<2>> & GetTransFreeZone () { return transfreezone; }
 
   ///
   int GetPointNr (int ln, int endp) const { return lines.Get(ln).I(endp); }

libsrc/meshing/smoothing2.cpp

@@ -584,7 +584,7 @@ namespace netgen
     //  meshthis -> ProjectPoint (surfi, pp1);
     //  meshthis -> GetNormalVector (surfi, pp1, n);
 
-    static NgArray<Point2d> pts2d;
+    static NgArray<Point<2>> pts2d;
     pts2d.SetSize(mesh.GetNP());
 
     grad = 0;
@@ -655,7 +655,7 @@ namespace netgen
     //    pp1.Add2 (x.Get(1), t1, x.Get(2), t2);
     pp1 = ld.sp1 + x(0) * ld.t1 + x(1) * ld.t2;
 
-    static NgArray<Point2d> pts2d;
+    static NgArray<Point<2>> pts2d;
     pts2d.SetSize(mesh.GetNP());
 
     deriv = 0;

libsrc/visualization/mvdraw.hpp

@@ -97,9 +97,9 @@ namespace netgen
   {
     double scalex = 1., scaley = 1., shiftx = 0., shifty = 0.;
   public:
-    shared_ptr<NgArray<Point3d>> locpointsptr;
+    shared_ptr<NgArray<Point<3>>> locpointsptr;
     shared_ptr<NgArray<INDEX_2>> loclinesptr;
-    shared_ptr<NgArray<Point2d>> plainpointsptr;
+    shared_ptr<NgArray<Point<2>>> plainpointsptr;
     int oldnl;
     bool clearptr;
     VisualSceneSurfaceMeshing ();