2d Delaunay, array-iterators

libsrc/general/array.hpp

@@ -16,6 +16,54 @@ namespace netgen
 
 
 
+
+  template <typename TSIZE>
+  class ArrayRangeIterator
+  {
+    TSIZE ind;
+  public:
+    ArrayRangeIterator (TSIZE ai) : ind(ai) { ; }
+    ArrayRangeIterator operator++ (int) { return ind++; }
+    ArrayRangeIterator operator++ () { return ++ind; }
+    TSIZE operator*() const { return ind; }
+    bool operator != (ArrayRangeIterator d2) { return ind != d2.ind; }
+  };
+
+  /// a range of intergers
+  template <typename T>
+  class T_Range
+  {
+    T first, next;
+  public: 
+    T_Range (T f, T n) : first(f), next(n) {;} 
+    T Size() const { return next-first; }
+    T operator[] (T i) const { return first+i; }
+    bool Contains (T i) const { return ((i >= first) && (i < next)); }
+
+    ArrayRangeIterator<T> begin() const { return first; }
+    ArrayRangeIterator<T> end() const { return next; }
+  };
+
+
+  template <typename T, int BASE = 0, typename TIND = int>
+  class FlatArray;
+
+  template <typename T, int BASE, typename TIND>
+  class ArrayIterator
+  {
+    FlatArray<T,BASE,TIND> ar;
+    TIND ind;
+  public:
+    ArrayIterator (FlatArray<T,BASE,TIND> aar, TIND ai) : ar(aar), ind(ai) { ; }
+    ArrayIterator operator++ (int)  { return ArrayIterator(ar, ind++); }
+    ArrayIterator operator++ ()   { return ArrayIterator(ar, ++ind); }
+    T operator*() const { return ar[ind]; }
+    T & operator*() { return ar[ind]; }
+    bool operator != (ArrayIterator d2) { return ind != d2.ind; }
+  };
+
+
+
   /**
      A simple array container.
      Array represented by size and data-pointer.
@@ -24,7 +72,7 @@ namespace netgen
      Optional range check by macro RANGE_CHECK
   */
 
-  template <typename T, int BASE = 0, typename TIND = int>
+  template <typename T, int BASE, typename TIND>
   class FlatArray
   {
   protected:
@@ -42,8 +90,14 @@ namespace netgen
     /// the size
     int Size() const { return size; }
 
+    ArrayIterator<T,BASE,TIND> begin() const
+    { return ArrayIterator<T,BASE,TIND> (*this, BASE); }
+    ArrayIterator<T,BASE,TIND> end() const
+    { return ArrayIterator<T,BASE,TIND> (*this, BASE+size); }
+
     TIND Begin() const { return TIND(BASE); }
     TIND End() const { return TIND(size+BASE); }
+    T_Range<TIND> Range() const { return T_Range<TIND>(BASE, size+BASE); }
 
     /// Access array. BASE-based
     T & operator[] (TIND i) const

libsrc/geom2d/genmesh2d.cpp

@@ -428,6 +428,7 @@ namespace netgen
     mesh->CalcLocalH(mp.grading);
 
     int bnp = mesh->GetNP(); // boundary points
+    auto BndPntRange = mesh->Points().Range();
 
     int hquad = mp.quad;
 
@@ -532,10 +533,10 @@ namespace netgen
 	Array<int, PointIndex::BASE> compress(bnp);
 	compress = -1;
 	int cnt = 0;
-	for (PointIndex pi = PointIndex::BASE; pi < bnp+PointIndex::BASE; pi++)
+        for (PointIndex pi : BndPntRange)
 	  if ( (*mesh)[pi].GetLayer() == geometry.GetDomainLayer(domnr))
 	    {
-	      meshing.AddPoint ( (*mesh)[pi], pi);
+	      meshing.AddPoint ((*mesh)[pi], pi);
 	      cnt++;
 	      compress[pi] = cnt;
 	    }
@@ -546,35 +547,29 @@ namespace netgen
 	  {
 	    if ( (*mesh)[si].domin == domnr)
 	      {
-		meshing.AddBoundaryElement ( compress[(*mesh)[si][0]], 
+		meshing.AddBoundaryElement (compress[(*mesh)[si][0]], 
-					     compress[(*mesh)[si][1]], gi, gi);
+                                            compress[(*mesh)[si][1]], gi, gi);
 	      }
 	    if ( (*mesh)[si].domout == domnr)
 	      {
-		meshing.AddBoundaryElement ( compress[(*mesh)[si][1]],
+		meshing.AddBoundaryElement (compress[(*mesh)[si][1]],
-					     compress[(*mesh)[si][0]], gi, gi);
+                                            compress[(*mesh)[si][0]], gi, gi);
-		
 	      }
-
 	  }
 
-
+        // not complete, use at own risk ...
-	mp.checkoverlap = 0;
+        // meshing.Delaunay(*mesh, domnr, mp);
-
+        mp.checkoverlap = 0;
-	meshing.GenerateMesh (*mesh, mp, h, domnr);
+        meshing.GenerateMesh (*mesh, mp, h, domnr);
 
 	for (SurfaceElementIndex sei = oldnf; sei < mesh->GetNSE(); sei++)
 	  (*mesh)[sei].SetIndex (domnr);
 
-
 	// astrid
 	char * material;
-	geometry.GetMaterial( domnr, material );
+	geometry.GetMaterial (domnr, material);
-	if ( material )
+	if (material)
-	  {
+          mesh->SetMaterial (domnr, material);
-	    (*mesh).SetMaterial ( domnr,  material );
-	  }
-
       }
 
     mp.quad = hquad;

libsrc/meshing/adfront2.cpp

@@ -492,4 +492,72 @@ namespace netgen
     
     return ((cnt % 2) != 0);
   }
+
+  bool AdFront2 :: SameSide (const Point<2> & lp1, const Point<2> & lp2, 
+                             const Array<int> * testfaces) const
+  {
+    int cnt = 0;
+
+    if (testfaces)
+      {
+        for (int ii = 0; ii < testfaces->Size(); ii++)
+          if (lines[(*testfaces)[ii]].Valid())
+            {
+              int i = (*testfaces)[ii];
+              const Point<3> & p13d = points[lines[i].L().I1()].P();
+              const Point<3> & p23d = points[lines[i].L().I2()].P();
+              
+              Point<2> p1(p13d(0), p13d(1));
+              Point<2> p2(p23d(0), p23d(1));
+              
+              // p1 + alpha v = lp1 + beta vl
+              Vec<2> v = p2-p1;
+              Vec<2> vl = lp2 - lp1;
+              Mat<2,2> mat, inv;
+              Vec<2> rhs, sol;
+              mat(0,0) = v(0);
+              mat(1,0) = v(1);
+              mat(0,1) = -vl(0);
+              mat(1,1) = -vl(1);
+              rhs = lp1-p1;
+              
+              if (Det(mat) == 0) continue;
+              CalcInverse (mat, inv);
+              sol = inv * rhs;
+              if (sol(0) >= 0 && sol(0) <= 1 & sol(1) >= 0 && sol(1) <= 1)
+                { cnt++; }
+            }
+
+      }
+    else
+      {
+        for (int i = 0; i < lines.Size(); i++)
+          if (lines[i].Valid())
+            {
+              const Point<3> & p13d = points[lines[i].L().I1()].P();
+              const Point<3> & p23d = points[lines[i].L().I2()].P();
+              
+              Point<2> p1(p13d(0), p13d(1));
+              Point<2> p2(p23d(0), p23d(1));
+              
+              // p1 + alpha v = lp1 + beta vl
+              Vec<2> v = p2-p1;
+              Vec<2> vl = lp2 - lp1;
+              Mat<2,2> mat, inv;
+              Vec<2> rhs, sol;
+              mat(0,0) = v(0);
+              mat(1,0) = v(1);
+              mat(0,1) = -vl(0);
+              mat(1,1) = -vl(1);
+              rhs = lp1-p1;
+              
+              if (Det(mat) == 0) continue;
+              CalcInverse (mat, inv);
+              sol = inv * rhs;
+              if (sol(0) >= 0 && sol(0) <= 1 & sol(1) >= 0 && sol(1) <= 1)
+                { cnt++; }
+            }
+      }
+    return ((cnt % 2) == 0);
+  }
 }

libsrc/meshing/adfront2.hpp

@@ -262,11 +262,12 @@ public:
   bool Inside (const Point<2> & p) const;
 
   bool SameSide (const Point<2> & lp1, const Point<2> & lp2, 
-                 const Array<int> * /* testfaces */ = NULL) const
+                 const Array<int> * /* testfaces */ = NULL) const;
+  /*
   {
     return Inside (lp1) == Inside (lp2);
   }
-
+  */
 
   ///
   void SetStartFront ();

libsrc/meshing/boundarylayer.cpp

@@ -44,12 +44,12 @@ namespace netgen
       }
 
       for (i = 1; i <= np; i++)
-      {
+        {
-         if (bndnodes.Test(i))
+          if (bndnodes.Test(i))
             mapto.Elem(i) = mesh.AddPoint (mesh.Point (i));
-         else
+          else
             mapto.Elem(i) = 0;
-      }
+        }
 
       for (i = 1; i <= mesh.GetNSE(); i++)
       {
@@ -116,8 +116,8 @@ namespace netgen
 
       Vec3d Vect_A,Vect_B;
       
-      Vect_A = mesh.Point(el.PNum(Vertex_A)) - mesh.Point(el.PNum(Vertex));
+      Vect_A = mesh[el.PNum(Vertex_A)] - mesh[el.PNum(Vertex)];
-      Vect_B = mesh.Point(el.PNum(Vertex_B)) - mesh.Point(el.PNum(Vertex));
+      Vect_B = mesh[el.PNum(Vertex_B)] - mesh[el.PNum(Vertex)];
 
       SurfaceNormal = Cross(Vect_A,Vect_B);
       SurfaceNormal.Normalize();

libsrc/meshing/delaunay2d.cpp

@@ -6,8 +6,68 @@
 namespace netgen
 {
 
+
+
+  class DelaunayTrig
+  {
+    PointIndex pnums[3];
+    Point<3> c;
+    double r;
+    double rad2;
+  public:    
+    DelaunayTrig () { ; }
+    DelaunayTrig (int p1, int p2, int p3)
+    { pnums[0] = p1; pnums[1] = p2; pnums[2] = p3; }
+
+    PointIndex & operator[] (int j) { return pnums[j]; }
+    const PointIndex & operator[] (int j) const { return pnums[j]; }
+
+    void CalcCenter (Mesh & mesh)
+    {
+      Point<3> p1 = mesh[pnums[0]];
+      Point<3> p2 = mesh[pnums[1]];
+      Point<3> p3 = mesh[pnums[2]];
+      Vec<3> v1 = p2-p1;
+      Vec<3> v2 = p3-p1;
+      Mat<2,2> mat, inv;
+      mat(0,0) = v1*v1;
+      mat(0,1) = v1*v2;
+      mat(1,0) = v2*v1;
+      mat(1,1) = v2*v2;
+      Vec<2> rhs, sol;
+      rhs(0) = 0.5 * v1*v1;
+      rhs(1) = 0.5 * v2*v2;
+      CalcInverse (mat, inv);
+      sol = inv * rhs;
+      
+      c = p1 + sol(0) * v1 + sol(1) * v2;
+      rad2 = Dist2(c, p1);
+      r = sqrt(rad2); 
+    }
+
+    Point<3> Center() const { return c; }
+    double Radius2() const { return rad2; }
+    Box<3> BoundingBox() const { return Box<3> (c-Vec<3>(r,r,0.1), c+Vec<3>(r,r,0.1)); }
+  };
+
+  ostream & operator<< (ostream & ost, DelaunayTrig trig)
+  {
+    ost << trig[0] << "-" << trig[1] << "-" << trig[2] << endl;
+    return ost;
+  }
+
+  
   void Meshing2 :: BlockFillLocalH (Mesh & mesh, const MeshingParameters & mp)
   {
+    static int timer = NgProfiler::CreateTimer ("Meshing2::BlockFill");
+    static int timer1 = NgProfiler::CreateTimer ("Meshing2::BlockFill 1");
+    static int timer2 = NgProfiler::CreateTimer ("Meshing2::BlockFill 2");
+    static int timer3 = NgProfiler::CreateTimer ("Meshing2::BlockFill 3");
+    static int timer4 = NgProfiler::CreateTimer ("Meshing2::BlockFill 4");
+    NgProfiler::RegionTimer reg (timer);
+
+    NgProfiler::StartTimer (timer1);
+
     double filldist = mp.filldist;
     
     cout << "blockfill local h" << endl;
@@ -15,7 +75,7 @@ namespace netgen
     PrintMessage (3, "blockfill local h");
 
     Array<Point<3> > npoints;
-  
+    
     // adfront -> CreateTrees();
 
     Box<3> bbox ( Box<3>::EMPTY_BOX );
@@ -28,8 +88,7 @@ namespace netgen
 	const Point<3> & p1 = adfront->GetPoint(line.L().I1());
 	const Point<3> & p2 = adfront->GetPoint(line.L().I2());
 	
-	double hi = Dist (p1, p2);
+        maxh = max (maxh, Dist (p1, p2));
-	if (hi > maxh) maxh = hi;
 	
 	bbox.Add (p1);
 	bbox.Add (p2);
@@ -42,8 +101,12 @@ namespace netgen
     // Point<3> mpc = bbox.Center();
     bbox.Increase (bbox.Diam()/2);
     Box<3> meshbox = bbox;
-    
+
-    LocalH loch2 (bbox, 1);
+    NgProfiler::StopTimer (timer1);    
+    NgProfiler::StartTimer (timer2);    
+
+
+    LocalH loch2 (bbox, 1, 2);
     
     if (mp.maxh < maxh) maxh = mp.maxh;
     
@@ -75,7 +138,7 @@ namespace netgen
 	changed = false;
 	for (int i = 0; i < npoints.Size(); i++)
 	  {
-	    if (mesh.LocalHFunction().GetH(npoints[i]) > 1.5 * maxh)
+	    if (mesh.LocalHFunction().GetH(npoints[i]) > 1.2 * maxh)
 	      {
 		mesh.LocalHFunction().SetH (npoints[i], maxh);
 		changed = true;
@@ -84,6 +147,10 @@ namespace netgen
       }
     while (changed);
 
+    NgProfiler::StopTimer (timer2);    
+    NgProfiler::StartTimer (timer3);    
+
+
     if (debugparam.slowchecks)
       (*testout) << "Blockfill with points: " << endl;
     *testout << "loch = " << mesh.LocalHFunction() << endl;
@@ -112,6 +179,8 @@ namespace netgen
 	  }
       }
     
+    NgProfiler::StopTimer (timer3);    
+    NgProfiler::StartTimer (timer4);    
   
 
   // find outer points
@@ -141,7 +210,8 @@ namespace netgen
       }
     
     loch2.FindInnerBoxes (adfront, NULL);
-    
+
+      // outer points : smooth mesh-grading
     npoints.SetSize(0);
     loch2.GetOuterPoints (npoints);
     
@@ -154,21 +224,199 @@ namespace netgen
 	  }
       }  
 
+    NgProfiler::StopTimer (timer4);    
   }
+
+
+
   
   void Meshing2 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp)
   {
+    static int timer = NgProfiler::CreateTimer ("Meshing2::Delaunay - total");
+    static int timerstart = NgProfiler::CreateTimer ("Meshing2::Delaunay - start");
+    static int timerfinish = NgProfiler::CreateTimer ("Meshing2::Delaunay - finish");
+    static int timer1 = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental");
+    static int timer1a = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental a");
+    static int timer1b = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental b");
+    static int timer1c = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental c");
+    static int timer1d = NgProfiler::CreateTimer ("Meshing2::Delaunay - incremental d");
+    NgProfiler::RegionTimer reg (timer);
+
+
+
     cout << "2D Delaunay meshing (in progress)" << endl;
 
-    // int oldnp = mesh.GetNP();
-
-    cout << "np, old = " << mesh.GetNP() << endl;
 
     BlockFillLocalH (mesh, mp);
 
+    NgProfiler::StartTimer (timerstart);
 
-    cout << "np, now = " << mesh.GetNP() << endl;
+    // do the delaunay
+
+
+    // face bounding box:
+    Box<3> bbox (Box<3>::EMPTY_BOX);
+
+    for (int i = 0; i < adfront->GetNFL(); i++)
+      {
+	const FrontLine & line = adfront->GetLine(i);
+        bbox.Add (Point<3> (adfront->GetPoint (line.L()[0])));
+        bbox.Add (Point<3> (adfront->GetPoint (line.L()[1])));
+      }
+
+    for (int i = 0; i < mesh.LockedPoints().Size(); i++)
+      bbox.Add (mesh.Point (mesh.LockedPoints()[i]));
+
+    cout << "bbox = " << bbox << endl;
+
+    // external point
+    Vec<3> vdiag = bbox.PMax()-bbox.PMin();
     
+    auto old_points = mesh.Points().Range();
+    DelaunayTrig startel;
+    startel[0] = mesh.AddPoint (bbox.PMin() + Vec<3> (-8*vdiag(0), -8*vdiag(1), 0));
+    startel[1] = mesh.AddPoint (bbox.PMin() + Vec<3> (+8*vdiag(0), -8*vdiag(1), 0));
+    startel[2] = mesh.AddPoint (bbox.PMin() + Vec<3> (0, 8*vdiag(1), 0));
+
+    Box<3> hbox;
+    hbox.Set (mesh[startel[0]]);
+    hbox.Add (mesh[startel[1]]);
+    hbox.Add (mesh[startel[2]]);
+    Point<3> hp = mesh[startel[0]];
+    hp(2) = 1;  hbox.Add (hp);
+    hp(2) = -1; hbox.Add (hp);
+    Box3dTree searchtree(hbox);
+
+    Array<DelaunayTrig> tempels;
+    startel.CalcCenter (mesh);
+
+    tempels.Append (startel);
+    searchtree.Insert(startel.BoundingBox(), 0);
+
+    Array<int> closeels;
+    Array<int> intersecting;
+    Array<INDEX_2> edges;
+
+
+
+
+    // reorder points
+    Array<PointIndex, PointIndex::BASE, PointIndex> mixed(old_points.Size());
+    int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
+    int prim;
+  
+    {
+      int i = 0;
+      while (old_points.Size() % prims[i] == 0) i++;
+      prim = prims[i];
+    }
+
+    for (PointIndex pi : old_points)
+      mixed[pi] = PointIndex ( (prim * pi) % old_points.Size() + PointIndex::BASE );
+    
+    NgProfiler::StopTimer (timerstart);
+    NgProfiler::StartTimer (timer1);
+
+
+    for (PointIndex i1 : old_points)
+      {
+        PointIndex i = mixed[i1];
+
+        NgProfiler::StartTimer (timer1a);
+        Point<3> newp = mesh[i];
+        intersecting.SetSize(0);
+        edges.SetSize(0);
+
+        searchtree.GetIntersecting (newp, newp, closeels);
+        // for (int jj = 0; jj < closeels.Size(); jj++)
+        // for (int j = 0; j < tempels.Size(); j++)
+        for (int j : closeels)
+          {
+            if (tempels[j][0] < 0) continue;
+            Point<3> c = tempels[j].Center();
+            double r2 = tempels[j].Radius2();
+
+            bool inside = Dist2(mesh[i], c) < r2;
+            if (inside) intersecting.Append (j);
+          }
+
+        NgProfiler::StopTimer (timer1a);
+        NgProfiler::StartTimer (timer1b);
+
+        // find outer edges
+        for (auto j : intersecting)
+          {
+            const DelaunayTrig & trig = tempels[j];
+            for (int k = 0; k < 3; k++)
+              {
+                int p1 = trig[k];
+                int p2 = trig[(k+1)%3];
+                INDEX_2 edge(p1,p2);
+                edge.Sort();
+                bool found = false;
+                for (int l = 0; l < edges.Size(); l++)
+                  if (edges[l] == edge)
+                    {
+                      edges.Delete(l); 
+                      found = true;
+                      break;
+                    }
+                if (!found) edges.Append (edge);
+              }
+          }
+
+        NgProfiler::StopTimer (timer1b);
+        NgProfiler::StartTimer (timer1c);
+
+        /*
+        for (int j = intersecting.Size()-1; j >= 0; j--)
+          tempels.Delete (intersecting[j]);
+        */
+        for (int j : intersecting)
+          {
+            searchtree.DeleteElement (j);
+            tempels[j][0] = -1;
+            tempels[j][1] = -1;
+            tempels[j][2] = -1;
+          }
+
+        NgProfiler::StopTimer (timer1c);
+        NgProfiler::StartTimer (timer1d);
+
+        for (auto edge : edges)
+          {
+            DelaunayTrig trig (edge[0], edge[1], i);
+            trig.CalcCenter (mesh);
+            tempels.Append (trig);
+            searchtree.Insert(trig.BoundingBox(), tempels.Size()-1);
+          }
+
+        NgProfiler::StopTimer (timer1d);
+      }
+
+    NgProfiler::StopTimer (timer1);
+    NgProfiler::StartTimer (timerfinish);
+
+    for (DelaunayTrig & trig : tempels)
+      {
+        if (trig[0] < 0) continue;
+
+        Point<3> c = Center (mesh[trig[0]], mesh[trig[1]], mesh[trig[2]]);
+        if (!adfront->Inside (Point<2> (c(0),c(1)))) continue;
+
+        Vec<3> n = Cross (mesh[trig[1]]-mesh[trig[0]], 
+                          mesh[trig[2]]-mesh[trig[0]]);
+        if (n(2) < 0) Swap (trig[1], trig[2]);
+
+        Element2d el(trig[0], trig[1], trig[2]);
+        el.SetIndex (domainnr);
+        mesh.AddSurfaceElement (el);
+      }
+
+    for (PointIndex pi : mesh.Points().Range())
+      *testout << pi << ": " << mesh[pi].Type() << endl;
+
+    NgProfiler::StopTimer (timerfinish);
   }
 
 }

libsrc/meshing/improve2.cpp

@@ -722,9 +722,12 @@ namespace netgen
 
 	    if (should)
 	      {
-		// (*testout) << "combine !" << endl;
+                /*
-		// (*testout) << "bad1 = " << bad1 << ", bad2 = " << bad2 << endl;
+                (*testout) << "combine !" << endl;
-
+                (*testout) << "bad1 = " << bad1 << ", bad2 = " << bad2 << endl;
+                (*testout) << "illegal1 = " << illegal1 << ", illegal2 = " << illegal2 << endl;
+                (*testout) << "loch = " << loch << endl;
+                */
 
 		mesh[pi1] = pnew;
 		PointGeomInfo gi;

libsrc/meshing/localh.cpp

@@ -52,69 +52,36 @@ namespace netgen
   }
   
 
-  LocalH :: LocalH (const Point3d & pmin, const Point3d & pmax, double agrading)
+  LocalH :: LocalH (Point<3> pmin, Point<3> pmax, double agrading, int adimension)
+    : dimension(adimension)
   {
     double x1[3], x2[3];
     double hmax;
 
-    boundingbox = Box3d (pmin, pmax);
+    boundingbox = Box<3> (pmin, pmax);
     grading = agrading;
 
     // a small enlargement, non-regular points 
     double val = 0.0879;
-    for (int i = 1; i <= 3; i++)
+    for (int i = 0; i < dimension; i++)
       {
-	x1[i-1] = (1 + val * i) * pmin.X(i) - val * i * pmax.X(i);
+	x1[i] = (1 + val * (i+1)) * pmin(i) - val * (i+1) * pmax(i);
-	x2[i-1] = 1.1 * pmax.X(i) - 0.1 * pmin.X(i);
+	x2[i] = 1.1 * pmax(i) - 0.1 * pmin(i);
       }
+    for (int i = dimension; i < 3; i++)
+      x1[i] = x2[i] = 0;
 
     hmax = x2[0] - x1[0];
-    for (int i = 1; i <= 2; i++)
+    for (int i = 1; i < dimension; i++)
-      if (x2[i] - x1[i] > hmax)
+      hmax = max2(x2[i]-x1[i], hmax);
-	hmax = x2[i] - x1[i];
 
-    for (int i = 0; i <= 2; i++)
+    for (int i = 0; i < dimension; i++)
       x2[i] = x1[i] + hmax;
 
     root = new GradingBox (x1, x2);
     boxes.Append (root);
   }
 
-
-  LocalH :: LocalH (const Box<3> & box, double agrading)
-  {
-    Point3d pmin = box.PMin();
-    Point3d pmax = box.PMax();
-
-    double x1[3], x2[3];
-    double hmax;
-
-    boundingbox = Box3d (pmin, pmax);
-    grading = agrading;
-
-    // a small enlargement, non-regular points 
-    double val = 0.0879;
-    for (int i = 1; i <= 3; i++)
-      {
-	x1[i-1] = (1 + val * i) * pmin.X(i) - val * i * pmax.X(i);
-	x2[i-1] = 1.1 * pmax.X(i) - 0.1 * pmin.X(i);
-      }
-
-    hmax = x2[0] - x1[0];
-    for (int i = 1; i <= 2; i++)
-      if (x2[i] - x1[i] > hmax)
-	hmax = x2[i] - x1[i];
-
-    for (int i = 0; i <= 2; i++)
-      x2[i] = x1[i] + hmax;
-
-    root = new GradingBox (x1, x2);
-    boxes.Append (root);
-  }
-
-
-
-
   LocalH :: ~LocalH ()
   {
     root->DeleteChilds();
@@ -126,147 +93,227 @@ namespace netgen
     root->DeleteChilds();
   }
 
-  void LocalH :: SetH (const Point3d & p, double h)
+  void LocalH :: SetH (Point<3> p, double h)
   {
-    /*
+    if (dimension == 2)
-      (*testout) << "Set h at " << p << " to " << h << endl;
-      if (h < 1e-8)
       {
-      cout << "do not set h to " << h << endl;
+        if (fabs (p(0) - root->xmid[0]) > root->h2 ||
-      return;
+            fabs (p(1) - root->xmid[1]) > root->h2)
-      }
+          return;
-    */
+        
+        if (GetH(p) <= 1.2 * h) return;
+        
+        GradingBox * box = root;
+        GradingBox * nbox = root;
+        GradingBox * ngb;
+        int childnr;
+        double x1[3], x2[3];
+        
+        while (nbox)
+          {
+            box = nbox;
+            childnr = 0;
+            if (p(0) > box->xmid[0]) childnr += 1;
+            if (p(1) > box->xmid[1]) childnr += 2;
+            nbox = box->childs[childnr];
+          };
+        
+        while (2 * box->h2 > h)
+          {
+            childnr = 0;
+            if (p(0) > box->xmid[0]) childnr += 1;
+            if (p(1) > box->xmid[1]) childnr += 2;
+            
+            double h2 = box->h2;
+            if (childnr & 1)
+              {
+                x1[0] = box->xmid[0];
+                x2[0] = x1[0]+h2;   // box->x2[0];
+              }
+            else
+              {
+                x2[0] = box->xmid[0];
+                x1[0] = x2[0]-h2;   // box->x1[0];
+              }
+            
+            if (childnr & 2)
+              {
+                x1[1] = box->xmid[1];
+                x2[1] = x1[1]+h2;   // box->x2[1];
+              }
+            else
+              {
+                x2[1] = box->xmid[1];
+                x1[1] = x2[1]-h2;   // box->x1[1];
+              }
+            x1[2] = x2[2] = 0;
 
-    if (fabs (p.X() - root->xmid[0]) > root->h2 ||
+            ngb = new GradingBox (x1, x2);
-	fabs (p.Y() - root->xmid[1]) > root->h2 ||
+            box->childs[childnr] = ngb;
-	fabs (p.Z() - root->xmid[2]) > root->h2)
+            ngb->father = box;
-      return;
+            
+            boxes.Append (ngb);
+            box = box->childs[childnr];
+          }
 
-    /*      
+        box->hopt = h;
-	    if (p.X() < root->x1[0] || p.X() > root->x2[0] ||
+        
-	    p.Y() < root->x1[1] || p.Y() > root->x2[1] ||
+        double hbox = 2 * box->h2;  // box->x2[0] - box->x1[0];
-	    p.Z() < root->x1[2] || p.Z() > root->x2[2])
+        double hnp = h + grading * hbox;
-	    return;
+        
-    */
+        Point<3> np;
+        for (int i = 0; i < 2; i++)
+          {
+            np = p;
+            np(i) = p(i) + hbox;
+            SetH (np, hnp);
+            
+            np(i) = p(i) - hbox;
+            SetH (np, hnp);
+          }
 
-
-    if (GetH(p) <= 1.2 * h) return;
-
-
-    GradingBox * box = root;
-    GradingBox * nbox = root;
-    GradingBox * ngb;
-    int childnr;
-    double x1[3], x2[3];
-
-    while (nbox)
-      {
-	box = nbox;
-	childnr = 0;
-	if (p.X() > box->xmid[0]) childnr += 1;
-	if (p.Y() > box->xmid[1]) childnr += 2;
-	if (p.Z() > box->xmid[2]) childnr += 4;
-	nbox = box->childs[childnr];
-      };
-
-
-    while (2 * box->h2 > h)
-      {
-	childnr = 0;
-	if (p.X() > box->xmid[0]) childnr += 1;
-	if (p.Y() > box->xmid[1]) childnr += 2;
-	if (p.Z() > box->xmid[2]) childnr += 4;
-
-	double h2 = box->h2;
-	if (childnr & 1)
-	  {
-	    x1[0] = box->xmid[0];
-	    x2[0] = x1[0]+h2;   // box->x2[0];
-	  }
-	else
-	  {
-	    x2[0] = box->xmid[0];
-	    x1[0] = x2[0]-h2;   // box->x1[0];
-	  }
-
-	if (childnr & 2)
-	  {
-	    x1[1] = box->xmid[1];
-	    x2[1] = x1[1]+h2;   // box->x2[1];
-	  }
-	else
-	  {
-	    x2[1] = box->xmid[1];
-	    x1[1] = x2[1]-h2;   // box->x1[1];
-	  }
-
-	if (childnr & 4)
-	  {
-	    x1[2] = box->xmid[2];
-	    x2[2] = x1[2]+h2;  // box->x2[2];
-	  }
-	else
-	  {
-	    x2[2] = box->xmid[2];
-	    x1[2] = x2[2]-h2;  // box->x1[2];
-	  }
-
-	ngb = new GradingBox (x1, x2);
-	box->childs[childnr] = ngb;
-	ngb->father = box;
-
-	boxes.Append (ngb);
-	box = box->childs[childnr];
       }
 
-    box->hopt = h;
+    else
-
-
-    double hbox = 2 * box->h2;  // box->x2[0] - box->x1[0];
-    double hnp = h + grading * hbox;
-
-    Point3d np;
-    for (int i = 1; i <= 3; i++)
       {
-	np = p;
+        if (fabs (p(0) - root->xmid[0]) > root->h2 ||
-	np.X(i) = p.X(i) + hbox;
+            fabs (p(1) - root->xmid[1]) > root->h2 ||
-	SetH (np, hnp);
+            fabs (p(2) - root->xmid[2]) > root->h2)
+          return;
+        
+        if (GetH(p) <= 1.2 * h) return;
+        
+        GradingBox * box = root;
+        GradingBox * nbox = root;
+        GradingBox * ngb;
+        int childnr;
+        double x1[3], x2[3];
+        
+        while (nbox)
+          {
+            box = nbox;
+            childnr = 0;
+            if (p(0) > box->xmid[0]) childnr += 1;
+            if (p(1) > box->xmid[1]) childnr += 2;
+            if (p(2) > box->xmid[2]) childnr += 4;
+            nbox = box->childs[childnr];
+          };
+        
+        
+        while (2 * box->h2 > h)
+          {
+            childnr = 0;
+            if (p(0) > box->xmid[0]) childnr += 1;
+            if (p(1) > box->xmid[1]) childnr += 2;
+            if (p(2) > box->xmid[2]) childnr += 4;
+            
+            double h2 = box->h2;
+            if (childnr & 1)
+              {
+                x1[0] = box->xmid[0];
+                x2[0] = x1[0]+h2;   // box->x2[0];
+              }
+            else
+              {
+                x2[0] = box->xmid[0];
+                x1[0] = x2[0]-h2;   // box->x1[0];
+              }
+            
+            if (childnr & 2)
+              {
+                x1[1] = box->xmid[1];
+                x2[1] = x1[1]+h2;   // box->x2[1];
+              }
+            else
+              {
+                x2[1] = box->xmid[1];
+                x1[1] = x2[1]-h2;   // box->x1[1];
+              }
+            
+            if (childnr & 4)
+              {
+                x1[2] = box->xmid[2];
+                x2[2] = x1[2]+h2;  // box->x2[2];
+              }
+            else
+              {
+                x2[2] = box->xmid[2];
+                x1[2] = x2[2]-h2;  // box->x1[2];
+              }
+            
+            ngb = new GradingBox (x1, x2);
+            box->childs[childnr] = ngb;
+            ngb->father = box;
+            
+            boxes.Append (ngb);
+            box = box->childs[childnr];
+          }
 
-	np.X(i) = p.X(i) - hbox;
+        box->hopt = h;
-	SetH (np, hnp);
+        
+        double hbox = 2 * box->h2;  // box->x2[0] - box->x1[0];
+        double hnp = h + grading * hbox;
+        
+        Point<3> np;
+        for (int i = 0; i < 3; i++)
+          {
+            np = p;
+            np(i) = p(i) + hbox;
+            SetH (np, hnp);
+            
+            np(i) = p(i) - hbox;
+            SetH (np, hnp);
+          }
       }
   }
 
 
 
-  double LocalH :: GetH (const Point3d & x) const
+  double LocalH :: GetH (Point<3> x) const
   {
     const GradingBox * box = root;
-
+    if (dimension == 2)
-    while (1)
       {
-	int childnr = 0;
+        while (1)
-	if (x.X() > box->xmid[0]) childnr += 1;
+          {
-	if (x.Y() > box->xmid[1]) childnr += 2;
+            int childnr = 0;
-	if (x.Z() > box->xmid[2]) childnr += 4;
+            if (x(0) > box->xmid[0]) childnr += 1;
-
+            if (x(1) > box->xmid[1]) childnr += 2;
-	if (box->childs[childnr])
+            
-	  box = box->childs[childnr];
+            if (box->childs[childnr])
-	else
+              box = box->childs[childnr];
-	  return box->hopt;
+            else
+              return box->hopt;
+          }
       }
+    else
+      {
+        while (1)
+          {
+            int childnr = 0;
+            if (x(0) > box->xmid[0]) childnr += 1;
+            if (x(1) > box->xmid[1]) childnr += 2;
+            if (x(2) > box->xmid[2]) childnr += 4;
+            
+            if (box->childs[childnr])
+              box = box->childs[childnr];
+            else
+              return box->hopt;
+          }
+      }
+      
   }
 
 
   /// minimal h in box (pmin, pmax)
-  double LocalH :: GetMinH (const Point3d & pmin, const Point3d & pmax) const
+  double LocalH :: GetMinH (Point<3> pmin, Point<3> pmax) const
   { 
-    Point3d pmin2, pmax2;
+    Point<3> pmin2, pmax2;
-    for (int j = 1; j <= 3; j++)
+    for (int j = 0; j < 3; j++)
-      if (pmin.X(j) < pmax.X(j))
+      if (pmin(j) < pmax(j))
-	{ pmin2.X(j) = pmin.X(j); pmax2.X(j) = pmax.X(j); }
+	{ pmin2(j) = pmin(j); pmax2(j) = pmax(j); }
       else
-	{ pmin2.X(j) = pmax.X(j); pmax2.X(j) = pmin.X(j); }
+	{ pmin2(j) = pmax(j); pmax2(j) = pmin(j); }
 
     return GetMinHRec (pmin2, pmax2, root); 
   }
@@ -275,31 +322,64 @@ namespace netgen
   double LocalH :: GetMinHRec (const Point3d & pmin, const Point3d & pmax,
 			       const GradingBox * box) const
   {
-    double h2 = box->h2;
+    if (dimension == 2)
-    if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 ||
+      {
-	pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 ||
+        double h2 = box->h2;
-	pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2)
+        if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 ||
-      return 1e8;
+            pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2)
-      
+          return 1e8;
-    double hmin = 2 * box->h2; // box->x2[0] - box->x1[0];
+        
-  
+        double hmin = 2 * box->h2; // box->x2[0] - box->x1[0];
-    for (int i = 0; i < 8; i++)
+        
-      if (box->childs[i])
+        for (int i = 0; i < 8; i++)
-	hmin = min2 (hmin, GetMinHRec (pmin, pmax, box->childs[i]));
+          if (box->childs[i])
-
+            hmin = min2 (hmin, GetMinHRec (pmin, pmax, box->childs[i]));
-    return hmin;
+        
+        return hmin;
+      }
+    else
+      {
+        double h2 = box->h2;
+        if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 ||
+            pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 ||
+            pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2)
+          return 1e8;
+        
+        double hmin = 2 * box->h2; // box->x2[0] - box->x1[0];
+        
+        for (int i = 0; i < 8; i++)
+          if (box->childs[i])
+            hmin = min2 (hmin, GetMinHRec (pmin, pmax, box->childs[i]));
+        
+        return hmin;
+      }
   }
 
 
+
+
+
+
+
+
+
   void LocalH :: CutBoundaryRec (const Point3d & pmin, const Point3d & pmax,
 				 GradingBox * box)
   {
     double h2 = box->h2;
-    if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 ||
+    if (dimension == 2)
-	pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 ||
+      {
-	pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2)
+        if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 ||
-      return;
+            pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2)
-
+          return;
+      }
+    else
+      {
+        if (pmax.X() < box->xmid[0]-h2 || pmin.X() > box->xmid[0]+h2 ||
+            pmax.Y() < box->xmid[1]-h2 || pmin.Y() > box->xmid[1]+h2 ||
+            pmax.Z() < box->xmid[2]-h2 || pmin.Z() > box->xmid[2]+h2)
+          return;
+      }
 
     box->flags.cutboundary = 1;
     for (int i = 0; i < 8; i++)
@@ -309,10 +389,13 @@ namespace netgen
 
 
 
-
   void LocalH :: FindInnerBoxes (AdFront3 * adfront,
 				 int (*testinner)(const Point3d & p1))
   {
+    static int timer = NgProfiler::CreateTimer ("LocalH::FindInnerBoxes");
+    NgProfiler::RegionTimer reg (timer);
+
+
     int nf = adfront->GetNF();
 
     for (int i = 0; i < boxes.Size(); i++)
@@ -327,11 +410,11 @@ namespace netgen
 
 
     root->flags.pinner = !adfront->SameSide (rpmid, rx2);
-  
+    
     if (testinner)
       (*testout) << "inner = " << root->flags.pinner << " =?= " 
 		 << testinner(Point3d(root->xmid[0], root->xmid[1], root->xmid[2])) << endl;
-  
+
     Array<int> faceinds(nf);
     Array<Box3d> faceboxes(nf);
 
@@ -461,13 +544,11 @@ namespace netgen
 
 
 
-
-
-
   void LocalH :: FindInnerBoxes (AdFront2 * adfront,
 				 int (*testinner)(const Point<2> & p1))
   {
-    int nf = adfront->GetNFL();
+    static int timer = NgProfiler::CreateTimer ("LocalH::FindInnerBoxes 2d");
+    NgProfiler::RegionTimer reg (timer);
 
     for (int i = 0; i < boxes.Size(); i++)
       boxes[i] -> flags.isinner = 0;
@@ -485,7 +566,9 @@ namespace netgen
     if (testinner)
       (*testout) << "inner = " << root->flags.pinner << " =?= "
 		 << testinner(rpmid) << endl;
-  
+
+
+    int nf = adfront->GetNFL();
     Array<int> faceinds(nf);
     Array<Box<3> > faceboxes(nf);
 
@@ -497,7 +580,6 @@ namespace netgen
 	const FrontLine & line = adfront->GetLine(i);
 	faceboxes[i].Set (adfront->GetPoint (line.L().I1()));
 	faceboxes[i].Add (adfront->GetPoint (line.L().I2()));
-
       }
   
     for (int i = 0; i < 8; i++)
@@ -512,43 +594,40 @@ namespace netgen
 		      Array<int> & faceinds, int nfinbox)
   {
     if (!box) return;
-  
+    
     GradingBox * father = box -> father;
   
-    Point3d c(box->xmid[0], box->xmid[1], box->xmid[2]);
+    Point3d c(box->xmid[0], box->xmid[1], 0); // box->xmid[2]);
     Vec3d v(box->h2, box->h2, box->h2);
     Box3d boxc(c-v, c+v);
 
-    Point3d fc(father->xmid[0], father->xmid[1], father->xmid[2]);
+    Point3d fc(father->xmid[0], father->xmid[1], 0); // father->xmid[2]);
     Vec3d fv(father->h2, father->h2, father->h2);
     Box3d fboxc(fc-fv, fc+fv);
-
     Box3d boxcfc(c,fc);
 
     ArrayMem<int, 100> faceused;
     ArrayMem<int, 100> faceused2;
     ArrayMem<int, 100> facenotused;
 
-
+    for (int j = 0; j < nfinbox; j++)
-    for (int j = 1; j <= nfinbox; j++)
       {
 	//      adfront->GetFaceBoundingBox (faceinds.Get(j), facebox);
-	const Box3d & facebox = faceboxes.Get(faceinds.Get(j));
+	const Box3d & facebox = faceboxes[faceinds[j]];
   
 	if (boxc.Intersect (facebox))
-	  faceused.Append(faceinds.Get(j));
+	  faceused.Append(faceinds[j]);
 	else
-	  facenotused.Append(faceinds.Get(j));
+	  facenotused.Append(faceinds[j]);
 
 	if (boxcfc.Intersect (facebox))
-	  faceused2.Append (faceinds.Get(j));
+	  faceused2.Append (faceinds[j]);
       }
   
-    for (int j = 1; j <= faceused.Size(); j++)
+    for (int j = 0; j < faceused.Size(); j++)
-      faceinds.Elem(j) = faceused.Get(j);
+      faceinds[j] = faceused[j];
-    for (int j = 1; j <= facenotused.Size(); j++)
+    for (int j = 0; j < facenotused.Size(); j++)
-      faceinds.Elem(j+faceused.Size()) = facenotused.Get(j);
+      faceinds[j+faceused.Size()] = facenotused[j];
-
   
     if (!father->flags.cutboundary)
       {
@@ -560,12 +639,15 @@ namespace netgen
 	Point3d cf(father->xmid[0], father->xmid[1], father->xmid[2]);
       
 	if (father->flags.isinner)
-	  box->flags.pinner = 1;
+          {
+            box->flags.pinner = 1;
+          }
 	else
 	  {
 	    Point<2> c2d (c.X(), c.Y());
 	    Point<2> cf2d (cf.X(), cf.Y());
-	    if (adfront->SameSide (c2d, cf2d, &faceused2))
+            bool sameside = adfront->SameSide (c2d, cf2d, &faceused2);
+            if (sameside)
 	      box->flags.pinner = father->flags.pinner;
 	    else
 	      box->flags.pinner = 1 - father->flags.pinner;
@@ -655,9 +737,19 @@ namespace netgen
 
   void LocalH :: GetInnerPoints (Array<Point<3> > & points)
   {
-    for (int i = 0; i < boxes.Size(); i++)
+    if (dimension == 2)
-      if (boxes[i] -> flags.isinner)
+      {
-	points.Append ( boxes[i] -> PMid() );
+        for (int i = 0; i < boxes.Size(); i++)
+          if (boxes[i] -> flags.isinner && boxes[i] -> HasChilds())
+            points.Append ( boxes[i] -> PMid() );
+      }
+    else
+      {
+        for (int i = 0; i < boxes.Size(); i++)
+          if (boxes[i] -> flags.isinner)
+            points.Append ( boxes[i] -> PMid() );
+      }
+          
   }
 
 

libsrc/meshing/localh.hpp

@@ -45,6 +45,13 @@ namespace netgen
     Point<3> PMid() const { return Point<3> (xmid[0], xmid[1], xmid[2]); }
     double H2() const { return h2; }
 
+    bool HasChilds() const
+    {
+      for (int i = 0; i < 8; i++)
+        if (childs[i]) return true;
+      return false;
+    }
+    
     friend class LocalH;
 
     static BlockAllocator ball;
@@ -67,12 +74,15 @@ namespace netgen
     ///
     Array<GradingBox*> boxes;
     ///
-    Box3d boundingbox;
+    Box<3> boundingbox;
+    /// octree or quadtree
+    int dimension;
   public:
     ///
-    LocalH (const Point3d & pmin, const Point3d & pmax, double grading);
+    LocalH (Point<3> pmin, Point<3> pmax, double grading, int adimension = 3);
     ///
-    LocalH (const Box<3> & box, double grading);
+    LocalH (const Box<3> & box, double grading, int adimension = 3)
+      : LocalH (box.PMin(), box.PMax(), grading, adimension) { ; }
     ///
     ~LocalH();
     ///
@@ -80,15 +90,16 @@ namespace netgen
     ///
     void SetGrading (double agrading) { grading = agrading; }
     ///
-    void SetH (const Point3d & x, double h);
+    void SetH (Point<3> x, double h);
     ///
-    double GetH (const Point3d & x) const;
+    double GetH (Point<3> x) const;
     /// minimal h in box (pmin, pmax)
-    double GetMinH (const Point3d & pmin, const Point3d & pmax) const;
+    double GetMinH (Point<3> pmin, Point<3> pmax) const;
 
     /// mark boxes intersecting with boundary-box
     // void CutBoundary (const Point3d & pmin, const Point3d & pmax)
     // { CutBoundaryRec (pmin, pmax, root); }
+
     void CutBoundary (const Box<3> & box)
     { CutBoundaryRec (box.PMin(), box.PMax(), root); }
   
@@ -117,7 +128,7 @@ namespace netgen
     void Convexify ();
     ///
     int GetNBoxes () { return boxes.Size(); } 
-    const Box3d & GetBoundingBox () const
+    const Box<3> & GetBoundingBox () const
     { return boundingbox; }
     ///
     void PrintMemInfo (ostream & ost) const;

libsrc/meshing/meshclass.cpp

@@ -2416,19 +2416,19 @@ namespace netgen
 
 
 
-  void Mesh :: SetLocalH (const Point3d & pmin, const Point3d & pmax, double grading)
+  void Mesh :: SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading)
   {
-    Point3d c = Center (pmin, pmax);
+    using netgen::Point;
-    double d = max3 (pmax.X()-pmin.X(),
+    Point<3> c = Center (pmin, pmax);
-                     pmax.Y()-pmin.Y(),
+    double d = max3 (pmax(0)-pmin(0),
-                     pmax.Z()-pmin.Z());
+                     pmax(1)-pmin(1),
+                     pmax(2)-pmin(2));
     d /= 2;
-    Point3d pmin2 = c - Vec3d (d, d, d);
+    Point<3> pmin2 = c - Vec<3> (d, d, d);
-    Point3d pmax2 = c + Vec3d (d, d, d);
+    Point<3> pmax2 = c + Vec<3> (d, d, d);
-
 
     delete lochfunc;
-    lochfunc = new LocalH (pmin2, pmax2, grading);
+    lochfunc = new LocalH (pmin2, pmax2, grading, dimension);
   }
 
   void Mesh :: RestrictLocalH (const Point3d & p, double hloc)

libsrc/meshing/meshclass.hpp

@@ -376,7 +376,7 @@ namespace netgen
     /// Calculates localh 
     DLL_HEADER void CalcLocalH (double grading);
     ///
-    DLL_HEADER void SetLocalH (const Point3d & pmin, const Point3d & pmax, double grading);
+    DLL_HEADER void SetLocalH (netgen::Point<3> pmin, netgen::Point<3> pmax, double grading);
     ///
     DLL_HEADER void RestrictLocalH (const Point3d & p, double hloc);
     ///
@@ -682,10 +682,8 @@ namespace netgen
 	  area += Cross ( mesh[sel[1]]-mesh[sel[0]],
 			  mesh[sel[2]]-mesh[sel[0]] ).Length() / 2;
 	else
-	  area += Cross (Vec3d (mesh.Point (sel.PNum(1)),
+	  area += Cross (Vec3d (mesh[sel.PNum(1)], mesh[sel.PNum(3)]),
-				mesh.Point (sel.PNum(3))),
+			 Vec3d (mesh[sel.PNum(1)], mesh[sel.PNum(4)])).Length() / 2;;
-			 Vec3d (mesh.Point (sel.PNum(1)),
-				mesh.Point (sel.PNum(4)))).Length() / 2;;
       }
       void ReCalc ()
       {

libsrc/meshing/meshfunc.cpp

@@ -172,9 +172,12 @@ namespace netgen
 
            mesh3d.FindOpenElements(k);
 
-
+           /*
            for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
               meshing.AddPoint (mesh3d[pi], pi);
+           */
+           for (PointIndex pi : mesh3d.Points().Range())
+              meshing.AddPoint (mesh3d[pi], pi);
 
            for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
               meshing.AddBoundaryElement (mesh3d.OpenElement(i));

libsrc/meshing/meshfunc2d.cpp

@@ -23,28 +23,28 @@ namespace netgen
 	    case 's': 
 	      {  // topological swap
 		MeshOptimize2d meshopt;
-		meshopt.SetMetricWeight (0);
+                meshopt.SetMetricWeight (mp.elsizeweight);
 		meshopt.EdgeSwapping (mesh, 0);
 		break;
 	      }
 	    case 'S': 
 	      {  // metric swap
 		MeshOptimize2d meshopt;
-		meshopt.SetMetricWeight (0);
+                meshopt.SetMetricWeight (mp.elsizeweight);
 		meshopt.EdgeSwapping (mesh, 1);
 		break;
 	      }
 	    case 'm': 
 	      {
 		MeshOptimize2d meshopt;
-		meshopt.SetMetricWeight (1);
+                meshopt.SetMetricWeight (mp.elsizeweight);
 		meshopt.ImproveMesh(mesh, mp);
 		break;
 	      }
 	    case 'c': 
 	      {
 		MeshOptimize2d meshopt;
-		meshopt.SetMetricWeight (0.2);
+                meshopt.SetMetricWeight (mp.elsizeweight);
 		meshopt.CombineImprove(mesh);
 		break;
 	      }

libsrc/meshing/meshtype.hpp

@@ -120,12 +120,7 @@ namespace netgen
     PointIndex () { ; }
     PointIndex (int ai) : i(ai) { ; }
     PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
-    // PointIndex & operator= (int ai) { i = ai; return *this; }
     operator int () const { return i; }
-    // int GetInt () const { return i; }
-    // PointIndex operator+ (int i2) { return PointIndex (i+i2); }
-    // PointIndex operator++ (int) { int hi = i; i++; return PointIndex(hi); }
-    // PointIndex operator-- (int) { int hi = i; i--; return PointIndex(hi); }
     PointIndex operator++ (int) { PointIndex hi(*this); i++; return hi; }
     PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; }
     PointIndex operator++ () { i++; return *this; }

libsrc/meshing/smoothing2.cpp

@@ -743,7 +743,6 @@ namespace netgen
 
     Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
-
     bool mixed = 0;
     for (int i = 0; i < seia.Size(); i++)
       if (mesh[seia[i]].GetNP() != 3)