diff --git a/libsrc/general/hashtabl.hpp b/libsrc/general/hashtabl.hpp
index bd825182..fb37eb38 100644
--- a/libsrc/general/hashtabl.hpp
+++ b/libsrc/general/hashtabl.hpp
@@ -87,7 +87,7 @@ public:
   ///
   void PrintStat (ostream & ost) const;
   void BaseSetSize(int s) {hash.SetSize(s);}
-protected:
+  //protected:
   ///
   int HashValue (const INDEX_2 & ind) const
     {
@@ -96,8 +96,7 @@ protected:
   ///
   int Position (int bnr, const INDEX_2 & ind) const
   {
-    int i;
-    for (i = 1; i <= hash.EntrySize (bnr); i++)
+    for (int i = 1; i <= hash.EntrySize (bnr); i++)
       if (hash.Get(bnr, i) == ind)
 	return i;
     return 0;
@@ -150,7 +149,7 @@ public:
   ///
   bool Used (const INDEX_2 & ahash) const
   {
-    return (Position (HashValue (ahash), ahash)) ? 1 : 0;
+    return Position (HashValue (ahash), ahash) > 0;
   }
  ///
   int GetNBags () const
diff --git a/libsrc/meshing/adfront2.cpp b/libsrc/meshing/adfront2.cpp
index 81bef640..8bb4cc87 100644
--- a/libsrc/meshing/adfront2.cpp
+++ b/libsrc/meshing/adfront2.cpp
@@ -211,19 +211,6 @@ namespace netgen
   }
 
 
-  /*
-  void AdFront2 :: IncrementClass (int li)
-  {
-    lines[li].IncrementClass();
-  }
-
-
-  void AdFront2 :: ResetClass (int li)
-  {
-    lines[li].ResetClass();
-  }
-  */
-
   int AdFront2 :: SelectBaseLine (Point<3>  & p1, Point<3>  & p2, 
 				  const PointGeomInfo *& geominfo1,
 				  const PointGeomInfo *& geominfo2,
@@ -474,7 +461,7 @@ namespace netgen
     DenseMatrix a(2), ainv(2);
     Vector b(2), u(2);
     
-    // random numbers:
+    // quasi-random numbers:
     n(0) = 0.123871;
     n(1) = 0.15432;
     
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index e4a6f0ae..5d52ff15 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -5,6 +5,7 @@ namespace netgen
 {
 
   Mesh :: Mesh ()
+    : surfarea(*this)
   {
     // volelements.SetName ("vol elements");
     // surfelements.SetName ("surf elements");
@@ -301,6 +302,9 @@ namespace netgen
     surfelements.Last().next = facedecoding[el.index-1].firstelement;
     facedecoding[el.index-1].firstelement = si;
 
+    if (SurfaceArea().Valid())
+      SurfaceArea().Add (el);
+
     lock.UnLock();
     return si;
   }
@@ -2024,15 +2028,15 @@ namespace netgen
 
   void Mesh :: FindOpenSegments (int surfnr)
   {
-    int i, j, k;
+    // int i, j, k;
 
-    // new version, general elemetns
+    // new version, general elements
     // hash index: pnum1-2
     // hash data : surfnr,  surfel-nr (pos) or segment nr(neg)
     INDEX_2_HASHTABLE<INDEX_2> faceht(4 * GetNSE()+GetNSeg()+1);   
 
     PrintMessage (5, "Test Opensegments");
-    for (i = 1; i <= GetNSeg(); i++)
+    for (int i = 1; i <= GetNSeg(); i++)
       {
         const Segment & seg = LineSegment (i);
 
@@ -2052,7 +2056,7 @@ namespace netgen
       }
 
 
-    for (i = 1; i <= GetNSeg(); i++)
+    for (int i = 1; i <= GetNSeg(); i++)
       {
         const Segment & seg = LineSegment (i);
 
@@ -2067,15 +2071,17 @@ namespace netgen
           }
       }
 
+    // bool buggy = false;
+    // ofstream bout("buggy.out");
 
-    for (i = 1; i <= GetNSE(); i++)
+    for (int i = 1; i <= GetNSE(); i++)
       {
         const Element2d & el = SurfaceElement(i);
         if (el.IsDeleted()) continue;
 
         if (surfnr == 0 || el.GetIndex() == surfnr)
           {
-            for (j = 1; j <= el.GetNP(); j++)
+            for (int j = 1; j <= el.GetNP(); j++)
               {
                 INDEX_2 seg (el.PNumMod(j), el.PNumMod(j+1));
                 INDEX_2 data;
@@ -2093,9 +2099,40 @@ namespace netgen
                       }
                     else
                       {
-                        PrintSysError ("hash table si not fitting for segment: ",
+			// buggy = true;
+                        PrintWarning ("hash table si not fitting for segment: ",
                                        seg.I1(), "-", seg.I2(), " other = ",
                                        data.I2());
+			// cout << "me: index = " << el.GetIndex() << ", el = " << el << endl;
+
+			/*
+			bout << "has index = " << seg << endl;
+			bout << "hash value = " << faceht.HashValue (seg) << endl;
+
+			if (data.I2() > 0)
+			  {
+			    int io = data.I2();
+			    cout << "other trig: index = " << SurfaceElement(io).GetIndex()
+				 << ", el = " << SurfaceElement(io) << endl;
+			  }
+			else
+			  {
+			    cout << "other seg " << -data.I2() << ", si = " << data.I1() << endl;
+			  }
+
+			
+			bout << "me: index = " << el.GetIndex() << ", el = " << el << endl;
+			if (data.I2() > 0)
+			  {
+			    int io = data.I2();
+			    bout << "other trig: index = " << SurfaceElement(io).GetIndex()
+				 << ", el = " << SurfaceElement(io) << endl;
+			  } 
+			else
+			  {
+			    bout << "other seg " << -data.I2() << ", si = " << data.I1() << endl;
+			  }
+			*/
                       }
                   }
                 else
@@ -2110,10 +2147,35 @@ namespace netgen
           }
       }  
 
+    /*
+    if (buggy)
+      {
+	for (int i = 1; i <= GetNSeg(); i++)
+	  bout << "seg" << i << " " << LineSegment(i) << endl;
+
+	for (int i = 1; i <= GetNSE(); i++)
+	  bout << "sel" << i << " " << SurfaceElement(i) << " ind = " 
+	       << SurfaceElement(i).GetIndex() << endl;
+
+	bout << "hashtable: " << endl;
+	for (int j = 1; j <= faceht.GetNBags(); j++)
+	  {
+	    bout << "bag " << j << ":" << endl;
+	    for (int k = 1; k <= faceht.GetBagSize(j); k++)
+	      {
+		INDEX_2 i2, data;
+		faceht.GetData (j, k, i2, data);
+		bout << "key = " << i2 << ", data = " << data << endl;
+	      }
+	  }
+	exit(1);
+      }
+    */
+
     (*testout) << "open segments: " << endl;
     opensegments.SetSize(0);
-    for (i = 1; i <= faceht.GetNBags(); i++)
-      for (j = 1; j <= faceht.GetBagSize(i); j++)
+    for (int i = 1; i <= faceht.GetNBags(); i++)
+      for (int j = 1; j <= faceht.GetBagSize(i); j++)
         {
           INDEX_2 i2;
           INDEX_2 data;
@@ -2130,7 +2192,7 @@ namespace netgen
                 {
                   // segment due to triangle
                   const Element2d & el = SurfaceElement (data.I2());
-                  for (k = 1; k <= el.GetNP(); k++)
+                  for (int k = 1; k <= el.GetNP(); k++)
                     {
                       if (seg[0] == el.PNum(k))
                         seg.geominfo[0] = el.GeomInfoPi(k);
@@ -2184,16 +2246,16 @@ namespace netgen
       ptyps.Elem(seg[1]) = EDGEPOINT;
       }
     */
-    for (i = 1; i <= points.Size(); i++)
+    for (int i = 1; i <= points.Size(); i++)
       points.Elem(i).SetType(SURFACEPOINT);
 
-    for (i = 1; i <= GetNSeg(); i++)
+    for (int i = 1; i <= GetNSeg(); i++)
       {
         const Segment & seg = LineSegment (i);
         points[seg[0]].SetType(EDGEPOINT);
         points[seg[1]].SetType(EDGEPOINT);
       }
-    for (i = 1; i <= GetNOpenSegments(); i++)
+    for (int i = 1; i <= GetNOpenSegments(); i++)
       {
         const Segment & seg = GetOpenSegment (i);
         points[seg[0]].SetType (EDGEPOINT);
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 1c6d5746..2b373191 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -676,6 +676,45 @@ namespace netgen
     const class AnisotropicClusters & GetClusters () const
     { return *clusters; }
 
+
+    class CSurfaceArea
+    {
+      const Mesh & mesh;
+      bool valid;
+      double area;
+    public:
+      CSurfaceArea (const Mesh & amesh) 
+	: mesh(amesh), valid(false) { ; }
+
+      void Add (const Element2d & sel)
+      {
+	if (sel.GetNP() == 3)
+	  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)),
+				mesh.Point (sel.PNum(3))),
+			 Vec3d (mesh.Point (sel.PNum(1)),
+				mesh.Point (sel.PNum(4)))).Length() / 2;;
+      }
+      void ReCalc ()
+      {
+	area = 0;
+	for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
+	  Add (mesh[sei]);
+	valid = true;
+      }
+
+      operator double () const { return area; }
+      bool Valid() const { return valid; }
+    };
+
+    CSurfaceArea surfarea;
+    CSurfaceArea & SurfaceArea() { return surfarea; }
+    const CSurfaceArea & SurfaceArea() const { return surfarea; }
+
+
+
     int GetTimeStamp() const { return timestamp; }
     void SetNextTimeStamp() 
     { timestamp = NextTimeStamp(); }
diff --git a/libsrc/meshing/meshing2.cpp b/libsrc/meshing/meshing2.cpp
index 55539797..499546ad 100644
--- a/libsrc/meshing/meshing2.cpp
+++ b/libsrc/meshing/meshing2.cpp
@@ -200,8 +200,15 @@ namespace netgen
     static int timer1 = NgProfiler::CreateTimer ("surface meshing1");
     static int timer2 = NgProfiler::CreateTimer ("surface meshing2");
     static int timer3 = NgProfiler::CreateTimer ("surface meshing3");
+
+    static int ts1 = NgProfiler::CreateTimer ("surface meshing start 1");
+    static int ts2 = NgProfiler::CreateTimer ("surface meshing start 2");
+    static int ts3 = NgProfiler::CreateTimer ("surface meshing start 3");
+
+
     NgProfiler::RegionTimer reg (timer);
 
+    NgProfiler::StartTimer (ts1);
 
     Array<int> pindex, lindex;
     Array<int> delpoints, dellines;
@@ -267,7 +274,6 @@ namespace netgen
     trigsonnode = 0;
     illegalpoint = 0;
   
-
     double totalarea = Area ();
     double meshedarea = 0;
 
@@ -302,14 +308,15 @@ namespace netgen
 	box.Set ( mesh[sel[0]] );
 	box.Add ( mesh[sel[1]] );
 	box.Add ( mesh[sel[2]] );
-	// surfeltree.Insert (box, i);
 	surfeltree.Insert (box, seia[i]);
       }
 
-
-
+    NgProfiler::StopTimer (ts1);
+    NgProfiler::StartTimer (ts2);
 
     if (totalarea > 0 || maxarea > 0)
+      meshedarea = mesh.SurfaceArea();
+      /*
       for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
 	{
 	  const Element2d & sel = mesh[sei];
@@ -326,9 +333,12 @@ namespace netgen
 				      mesh.Point (sel.PNum(4)))).Length() / 2;;
 	  meshedarea += trigarea;
 	}
+      */
+      // cout << "meshedarea = " << meshedarea << " =?= " 
+      // << mesh.SurfaceArea() << endl;
 
-
-
+    NgProfiler::StopTimer (ts2);
+    NgProfiler::StartTimer (ts3);
 
     const char * savetask = multithread.task;
     multithread.task = "Surface meshing";
@@ -340,6 +350,7 @@ namespace netgen
 
     double meshedarea_before = meshedarea;
 
+    NgProfiler::StopTimer (ts3);
 
     while (!adfront ->Empty() && !multithread.terminate)
       {
diff --git a/libsrc/meshing/smoothing2.cpp b/libsrc/meshing/smoothing2.cpp
index fe47f8b8..729e4525 100644
--- a/libsrc/meshing/smoothing2.cpp
+++ b/libsrc/meshing/smoothing2.cpp
@@ -568,7 +568,6 @@ namespace netgen
 	faceindex = 0;
 	return;
       }
- 
 
     static int timer = NgProfiler::CreateTimer ("MeshSmoothing 2D");
     static int timer1 = NgProfiler::CreateTimer ("MeshSmoothing 2D start");
@@ -580,7 +579,6 @@ namespace netgen
 
     Opti2dLocalData ld;
 
-    // SurfaceElementIndex sei;
 
     Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
@@ -593,23 +591,48 @@ namespace netgen
 	  break;
 	}
 
-    int loci;
-    double fact;
-    int moveisok;
-
-    PointGeomInfo ngi;
-    Point3d origp;
-
-
-
-    Vec3d n1, n2;
-    Vector x(2), xedge(1);
+    Vector x(2);
 
     Array<MeshPoint, PointIndex::BASE> savepoints(mesh.GetNP());
 
     ld.uselocalh = mp.uselocalh;
 
+    Array<int, PointIndex::BASE> compress(mesh.GetNP());
+    Array<PointIndex> icompress;
+    for (int i = 0; i < seia.Size(); i++)
+      {
+	const Element2d & el = mesh[seia[i]];
+	for (int j = 0; j < el.GetNP(); j++)
+	  compress[el[j]] = -1;
+      }
+    for (int i = 0; i < seia.Size(); i++)
+      {
+	const Element2d & el = mesh[seia[i]];
+	for (int j = 0; j < el.GetNP(); j++)
+	  if (compress[el[j]] == -1)
+	    {
+	      compress[el[j]] = icompress.Size();
+	      icompress.Append(el[j]);
+	    }
+      }
+    Array<int> cnta(icompress.Size());
+    cnta = 0;
+    for (int i = 0; i < seia.Size(); i++)
+      {
+	const Element2d & el = mesh[seia[i]];
+	for (int j = 0; j < el.GetNP(); j++)
+	  cnta[compress[el[j]]]++;
+      }
+    TABLE<SurfaceElementIndex> elementsonpoint(cnta);
+    for (int i = 0; i < seia.Size(); i++)
+      {
+	const Element2d & el = mesh[seia[i]];
+	for (int j = 0; j < el.GetNP(); j++)
+	  elementsonpoint.Add (compress[el[j]], seia[i]);
+      }
 
+    
+    /*
     Array<int, PointIndex::BASE> nelementsonpoint(mesh.GetNP());
     nelementsonpoint = 0;
     for (int i = 0; i < seia.Size(); i++)
@@ -627,6 +650,8 @@ namespace netgen
 	for (int j = 0; j < el.GetNP(); j++)
 	  elementsonpoint.Add (el[j], seia[i]);
       }
+    */
+    
 
 
     ld.loch = mp.maxh;
@@ -645,7 +670,7 @@ namespace netgen
 
     /*
     int i, j, k;
-
+    Vector xedge(1);
       if (improveedges)
       for (i = 1; i <= mesh.GetNP(); i++)
       if (mesh.PointType(i) == EDGEPOINT)
@@ -695,6 +720,7 @@ namespace netgen
 	  
       if (surfi2 && !surfi3)
       {
+      Vec3d n1, n2;
       GetNormalVector (surfi, sp1, n1);
       GetNormalVector (surfi2, sp1, n2);
       t1 = Cross (n1, n2);
@@ -717,77 +743,87 @@ namespace netgen
     if (mesh.GetNP() > 1000)
       {
 	plotchar = '+';
-	modplot = 10;
+	modplot = 100;
       }
     if (mesh.GetNP() > 10000)
       {
 	plotchar = 'o';
-	modplot = 100;
+	modplot = 1000;
+      }
+    if (mesh.GetNP() > 100000)
+      {
+	plotchar = 'O';
+	modplot = 10000;
       }
     int cnt = 0;
 
 
     NgProfiler::StopTimer (timer1);
 
-
+    /*
     for (PointIndex pi = PointIndex::BASE; pi < mesh.GetNP()+PointIndex::BASE; pi++)
       if (mesh[pi].Type() == SURFACEPOINT)
-	{
-	  if (multithread.terminate)
-	    throw NgException ("Meshing stopped");
-	
-	  cnt++;
-	  if (cnt % modplot == 0 && writestatus)
-	    {
-	      printeddot = 1;
-	      PrintDot (plotchar);
-	    }
-
-	  if (elementsonpoint[pi].Size() == 0)
-	    continue;
-
-	
-	  ld.sp1 = mesh[pi];
-
-	  Element2d & hel = mesh[elementsonpoint[pi][0]];
-
-	  int hpi = 0;
-	  for (int j = 1; j <= hel.GetNP(); j++)
-	    if (hel.PNum(j) == pi)
+    */
+    for (int hi = 0; hi < icompress.Size(); hi++)
+      {
+	PointIndex pi = icompress[hi];
+	if (mesh[pi].Type() == SURFACEPOINT)
+	  {
+	    if (multithread.terminate)
+	      throw NgException ("Meshing stopped");
+	    
+	    cnt++;
+	    if (cnt % modplot == 0 && writestatus)
 	      {
-		hpi = j;
-		break;
+		printeddot = 1;
+		PrintDot (plotchar);
 	      }
-
-	  ld.gi1 = hel.GeomInfoPi(hpi);
-	  SelectSurfaceOfPoint (ld.sp1, ld.gi1);
-	  
-	  ld.locelements.SetSize(0);
-	  ld.locrots.SetSize (0);
-	  ld.lochs.SetSize (0);
+	    
+	    // if (elementsonpoint[pi].Size() == 0) continue;
+	    if (elementsonpoint[hi].Size() == 0) continue;
 	
-	  for (int j = 0; j < elementsonpoint[pi].Size(); j++)
-	    {
-	      SurfaceElementIndex sei = elementsonpoint[pi][j];
-	      const Element2d & bel = mesh[sei];
-	      ld.surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();
+	    ld.sp1 = mesh[pi];
 	    
-	      ld.locelements.Append (sei);
+	    // Element2d & hel = mesh[elementsonpoint[pi][0]];
+	    Element2d & hel = mesh[elementsonpoint[hi][0]];
 	    
-	      for (int k = 1; k <= bel.GetNP(); k++)
-		if (bel.PNum(k) == pi)
-		  {
-		    ld.locrots.Append (k);
-		    break;
-		  }
-	      
-	      if (ld.uselocalh)
+	    int hpi = 0;
+	    for (int j = 1; j <= hel.GetNP(); j++)
+	      if (hel.PNum(j) == pi)
 		{
-		  Point3d pmid = Center (mesh[bel[0]], mesh[bel[1]], mesh[bel[2]]);
-		  ld.lochs.Append (mesh.GetH(pmid));
+		  hpi = j;
+		  break;
 		}
-	    }
+
+	    ld.gi1 = hel.GeomInfoPi(hpi);
+	    SelectSurfaceOfPoint (ld.sp1, ld.gi1);
 	  
+	    ld.locelements.SetSize(0);
+	    ld.locrots.SetSize (0);
+	    ld.lochs.SetSize (0);
+	
+	    for (int j = 0; j < elementsonpoint[hi].Size(); j++)
+	      {
+		SurfaceElementIndex sei = elementsonpoint[hi][j];
+		const Element2d & bel = mesh[sei];
+		ld.surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();
+		
+		ld.locelements.Append (sei);
+		
+		for (int k = 1; k <= bel.GetNP(); k++)
+		  if (bel.PNum(k) == pi)
+		    {
+		      ld.locrots.Append (k);
+		      break;
+		    }
+		
+		if (ld.uselocalh)
+		  {
+		    Point3d pmid = Center (mesh[bel[0]], mesh[bel[1]], mesh[bel[2]]);
+		    ld.lochs.Append (mesh.GetH(pmid));
+		  }
+	      }
+	    
 	  GetNormalVector (ld.surfi, ld.sp1, ld.gi1, ld.normal);
 	  ld.t1 = ld.normal.GetNormal ();
 	  ld.t2 = Cross (ld.normal, ld.t1);
@@ -830,10 +866,10 @@ namespace netgen
 	    }
 
 	
-	  origp = mesh[pi];
-	  loci = 1;
-	  fact = 1;
-	  moveisok = 0;
+	  Point3d origp = mesh[pi];
+	  int loci = 1;
+	  double fact = 1;
+	  int moveisok = 0;
 
 	  // restore other points
 	  for (int j = 0; j < ld.locelements.Size(); j++)
@@ -867,6 +903,7 @@ namespace netgen
 	      // ProjectPoint (surfi, mesh[pi]);
 	      // moveisok = CalcPointGeomInfo(surfi, ngi, mesh[pi]); 
 
+	      PointGeomInfo ngi;
 	      ngi = ld.gi1;
 	      moveisok = ProjectPointGI (ld.surfi, mesh[pi], ngi);
 	      // point lies on same chart in stlsurface
@@ -883,10 +920,10 @@ namespace netgen
 	    
 	    }
 	}
-
+      }
     if (printeddot)
       PrintDot ('\n');
-  
+
     CheckMeshApproximation (mesh);
     mesh.SetNextTimeStamp();
   }
diff --git a/libsrc/stlgeom/meshstlsurface.cpp b/libsrc/stlgeom/meshstlsurface.cpp
index fcbb1e02..e6bc344c 100644
--- a/libsrc/stlgeom/meshstlsurface.cpp
+++ b/libsrc/stlgeom/meshstlsurface.cpp
@@ -16,10 +16,7 @@ namespace netgen
 static void STLFindEdges (STLGeometry & geom,
 			  class Mesh & mesh)
 {
-  int i, j;
-  double h;
-
-  h = mparam.maxh;
+  double h = mparam.maxh;
 
   // mark edge points:
   //int ngp = geom.GetNP();
@@ -33,23 +30,29 @@ static void STLFindEdges (STLGeometry & geom,
 
   PrintMessage(3,"Mesh Lines");
 
-  for (i = 1; i <= geom.GetNLines(); i++)
+  cout << geom.GetNLines() << " lines" << endl;
+  double totnp = 0;
+  for (int i = 1; i <= geom.GetNLines(); i++)
+    totnp += geom.GetLine(i)->NP();
+  cout << "avg np per line " << totnp/geom.GetNLines() << endl;
+
+  for (int i = 1; i <= geom.GetNLines(); i++)
     {
       meshlines.Append(geom.GetLine(i)->Mesh(geom.GetPoints(), meshpoints, h, mesh)); 
       SetThreadPercent(100.0 * (double)i/(double)geom.GetNLines());
     }
 
+  NgProfiler::Print(stdout);
+
   geom.meshpoints.SetSize(0); //testing
   geom.meshlines.SetSize(0);  //testing
-  int pim;
-  for (i = 1; i <= meshpoints.Size(); i++)
+  for (int i = 1; i <= meshpoints.Size(); i++)
     {
       geom.meshpoints.Append(meshpoints.Get(i)); //testing
-
-      pim = mesh.AddPoint(meshpoints.Get(i));
+      mesh.AddPoint(meshpoints.Get(i));
     }
   //(++++++++++++++testing
-  for (i = 1; i <= geom.GetNLines(); i++)
+  for (int i = 1; i <= geom.GetNLines(); i++)
     {
       geom.meshlines.Append(meshlines.Get(i));
     }
@@ -57,11 +60,11 @@ static void STLFindEdges (STLGeometry & geom,
 
   PrintMessage(7,"feed with edges");
 
-  for (i = 1; i <= meshlines.Size(); i++)
+  for (int i = 1; i <= meshlines.Size(); i++)
     {
       STLLine* line = meshlines.Get(i);
       (*testout) << "store line " << i << endl;
-      for (j = 1; j <= line->GetNS(); j++)
+      for (int j = 1; j <= line->GetNS(); j++)
 	{
 	  int p1, p2;
 	  
@@ -243,7 +246,6 @@ void STLSurfaceMeshing1 (STLGeometry & geom,
 int STLSurfaceMeshing (STLGeometry & geom,
 		       class Mesh & mesh)
 {
-  int i, j;
   PrintFnStart("Do Surface Meshing");
 
   geom.PrepareSurfaceMeshing();
@@ -256,7 +258,7 @@ int STLSurfaceMeshing (STLGeometry & geom,
 
   //  mesh.Save ("mesh.edges");
   
-  for (i = 1; i <= mesh.GetNSeg(); i++)
+  for (int i = 1; i <= mesh.GetNSeg(); i++)
     {
       const Segment & seg = mesh.LineSegment (i);
       if (seg.geominfo[0].trignum <= 0 || seg.geominfo[1].trignum <= 0)
@@ -298,14 +300,14 @@ int STLSurfaceMeshing (STLGeometry & geom,
 	      if (nopen)
 		{
 		  geom.ClearMarkedSegs();
-		  for (i = 1; i <= nopen; i++)
+		  for (int i = 1; i <= nopen; i++)
 		    {
 		      const Segment & seg = mesh.GetOpenSegment (i);
 		      geom.AddMarkedSeg(mesh.Point(seg[0]),mesh.Point(seg[1]));
 		    }
 
 		  geom.InitMarkedTrigs();
-		  for (i = 1; i <= nopen; i++)
+		  for (int i = 1; i <= nopen; i++)
 		    {
 		      const Segment & seg = mesh.GetOpenSegment (i);
 		      geom.SetMarkedTrig(seg.geominfo[0].trignum,1);
@@ -361,7 +363,7 @@ int STLSurfaceMeshing (STLGeometry & geom,
 
 	      // Open edge-segments will be refined !
 	      INDEX_2_HASHTABLE<int> openseght (nopen+1);
-	      for (i = 1; i <= mesh.GetNOpenSegments(); i++)
+	      for (int i = 1; i <= mesh.GetNOpenSegments(); i++)
 		{
 		  const Segment & seg = mesh.GetOpenSegment (i);
 		  INDEX_2 i2(seg[0], seg[1]);
@@ -379,7 +381,7 @@ int STLSurfaceMeshing (STLGeometry & geom,
 	      INDEX_2_HASHTABLE<int> newpht(100);
 
 	      int nsegold = mesh.GetNSeg();
-	      for (i = 1; i <= nsegold; i++)
+	      for (int i = 1; i <= nsegold; i++)
 		{
 		  Segment seg = mesh.LineSegment(i);
 		  INDEX_2 i2(seg[0], seg[1]);
@@ -462,7 +464,7 @@ int STLSurfaceMeshing (STLGeometry & geom,
 
 	  geom.InitMarkedTrigs();
 
-	  for (i = 1; i <= mesh.GetNSE(); i++)
+	  for (int i = 1; i <= mesh.GetNSE(); i++)
 	    if (mesh.SurfaceElement(i).BadElement())
 	      {
 		int trig = mesh.SurfaceElement(i).PNum(1);
@@ -477,10 +479,10 @@ int STLSurfaceMeshing (STLGeometry & geom,
 
 	  // was commented:
 
-	  for (i = 1; i <= mesh.GetNSE(); i++)
+	  for (int i = 1; i <= mesh.GetNSE(); i++)
 	    if (mesh.SurfaceElement(i).BadElement())
 	      {
-		for (j = 1; j <= 3; j++)
+		for (int j = 1; j <= 3; j++)
 		  {
 		    refpts.Append (mesh.Point (mesh.SurfaceElement(i).PNum(j)));
 		    refh.Append (mesh.GetH (refpts.Last()) / 2);
@@ -489,7 +491,7 @@ int STLSurfaceMeshing (STLGeometry & geom,
 	      }
 	  	  
 	  // delete wrong oriented element
-	  for (i = 1; i <= mesh.GetNSE(); i++)
+	  for (int i = 1; i <= mesh.GetNSE(); i++)
 	    {
 	      const Element2d & el = mesh.SurfaceElement(i);
 	      if (!el.PNum(1))
@@ -509,7 +511,7 @@ int STLSurfaceMeshing (STLGeometry & geom,
 	    }
 	  // end comments
 
-	  for (i = 1; i <= refpts.Size(); i++)
+	  for (int i = 1; i <= refpts.Size(); i++)
 	    mesh.RestrictLocalH (refpts.Get(i), refh.Get(i));
 
 	  mesh.RemoveOneLayerSurfaceElements();
@@ -550,89 +552,121 @@ void STLSurfaceMeshing1 (STLGeometry & geom,
 			 class Mesh & mesh,
 			 int retrynr)
 {
-  int i, j;
-  double h;
-  
-  
-  h = mparam.maxh;
+  double h = mparam.maxh;
 
   mesh.FindOpenSegments();
   
   Array<int> spiralps(0);
   spiralps.SetSize(0);
-  for (i = 1; i <= geom.GetNP(); i++)
+  for (int i = 1; i <= geom.GetNP(); i++)
     {
       if (geom.GetSpiralPoint(i)) {spiralps.Append(i);}
     }
   
   PrintMessage(7,"NO spiralpoints = ", spiralps.Size());
   //int spfound;
-  int sppointnum;
   int spcnt = 0;
 
   Array<int> meshsp(mesh.GetNP());
-  for (i = 1; i <= mesh.GetNP(); i++)
-    {
-      meshsp.Elem(i) = 0;
-      for (j = 1; j <= spiralps.Size(); j++)
-	if (Dist2(geom.GetPoint(spiralps.Get(j)), mesh.Point(i)) < 1e-20) 
-	  meshsp.Elem(i) = spiralps.Get(j);
-    }
+  meshsp = 0;
+  for (int i = 1; i <= mesh.GetNP(); i++)
+    for (int j = 1; j <= spiralps.Size(); j++)
+      if (Dist2(geom.GetPoint(spiralps.Get(j)), mesh.Point(i)) < 1e-20) 
+	meshsp.Elem(i) = spiralps.Get(j);
 
+  Array<int, 1> opensegsperface(mesh.GetNFD());
+  opensegsperface = 0;
 
-  Array<int> opensegsperface(mesh.GetNFD());
-  for (i = 1; i <= mesh.GetNFD(); i++)
-    opensegsperface.Elem(i) = 0;
-  for (i = 1; i <= mesh.GetNOpenSegments(); i++)
-    {
-      int si = mesh.GetOpenSegment (i).si;
-      if (si >= 1 && si <= mesh.GetNFD())
-	{
-	  opensegsperface.Elem(si)++;
-	}
-      else
-	{
-	  cerr << "illegal face index" << endl;
-	}
-    }
-
+  for (int i = 1; i <= mesh.GetNOpenSegments(); i++)
+    opensegsperface[mesh.GetOpenSegment(i).si]++;
 
   double starttime = GetTime ();
+  mesh.SurfaceArea().ReCalc();
+
+  int oldnp = mesh.GetNP();
+  Array<int,PointIndex::BASE> compress(mesh.GetNP());
+  Array<PointIndex> icompress;
 
   for (int fnr = 1; fnr <= mesh.GetNFD(); fnr++)
-    if (opensegsperface.Get(fnr))
-      {
-	if (multithread.terminate)
-	  return;
-	
-	PrintMessage(5,"Meshing surface ", fnr, "/", mesh.GetNFD());
-	MeshingSTLSurface meshing (geom, mparam);
-	
-	meshing.SetStartTime (starttime);
-	
-	for (i = 1; i <= mesh.GetNP(); i++)
+    {
+      if (!opensegsperface[fnr]) continue;
+      if (multithread.terminate) return;
+
+      PrintMessage(5,"Meshing surface ", fnr, "/", mesh.GetNFD());
+      MeshingSTLSurface meshing (geom, mparam);
+      meshing.SetStartTime (starttime);
+      
+      compress = 0;
+      icompress.SetSize(0);
+      int cntused = 0;
+      for (int i = 1; i <= meshsp.Size(); i++)
+	if (meshsp.Elem(i))
 	  {
-	    /*
-	      spfound = 0;
-	      for (j = 1; j <= spiralps.Size(); j++)
-	      {
-	      if (Dist2(geom.GetPoint(spiralps.Get(j)),mesh.Point(i)) < 1e-20) 
-		{spfound =  1; sppointnum = spiralps.Get(j);}
+	    compress[i] = ++cntused;
+	    icompress.Append(i);
+	  }
+
+      for (int i = 1; i <= mesh.GetNOpenSegments(); i++)
+	{
+	  const Segment & seg = mesh.GetOpenSegment (i);
+	  if (seg.si == fnr)
+	    for (int j = 0; j < 2; j++)
+	      if (compress[seg[j]] == 0)
+		{
+		  compress[seg[j]] = ++cntused;
+		  icompress.Append(seg[j]);
 		}
-	    */
-	    sppointnum = 0;
-	    if (i <= meshsp.Size())
-	      sppointnum = meshsp.Get(i);
-	    
-	  //spfound = 0;
+	}
+
+      // for (int i = 1; i <= oldnp; i++) compress[i] = i;
+
+      /*
+      // for (int i = 1; i <= mesh.GetNP(); i++)
+      for (int i = 1; i <= oldnp; i++)
+	{
+	  int sppointnum = meshsp.Get(i);
 	  if (sppointnum)
 	    {
 	      MultiPointGeomInfo mgi;
-  
+	      
 	      int ntrigs = geom.NOTrigsPerPoint(sppointnum);
 	      spcnt++;
 	      
-	      for (j = 0; j < ntrigs; j++)
+	      for (int j = 0; j < ntrigs; j++)
+		{
+		  PointGeomInfo gi;
+		  gi.trignum = geom.TrigPerPoint(sppointnum, j+1);
+		  mgi.AddPointGeomInfo (gi);
+		}
+	      
+	      // Einfuegen von ConePoint: Point bekommt alle
+	      // Dreiecke (werden dann intern kopiert)
+	      // Ein Segment zum ConePoint muss vorhanden sein !!!
+	      
+	      if (compress[i])
+		meshing.AddPoint (mesh.Point(i), i, &mgi);
+	    }
+	    else
+	      {
+		if (compress[i])
+		  meshing.AddPoint (mesh.Point(i), i);
+	      }
+	  }
+      */
+
+      for (int hi = 0; hi < icompress.Size(); hi++)
+	{
+	  PointIndex i = icompress[hi];
+
+	  int sppointnum = meshsp.Get(i);
+	  if (sppointnum)
+	    {
+	      MultiPointGeomInfo mgi;
+	      
+	      int ntrigs = geom.NOTrigsPerPoint(sppointnum);
+	      spcnt++;
+	      
+	      for (int j = 0; j < ntrigs; j++)
 		{
 		  PointGeomInfo gi;
 		  gi.trignum = geom.TrigPerPoint(sppointnum, j+1);
@@ -644,35 +678,36 @@ void STLSurfaceMeshing1 (STLGeometry & geom,
 	      // Ein Segment zum ConePoint muss vorhanden sein !!!
 	      
 	      meshing.AddPoint (mesh.Point(i), i, &mgi);
-	      
 	    }
 	  else
 	    {
 	      meshing.AddPoint (mesh.Point(i), i);
 	    }
 	}
-      
-      
-      for (i = 1; i <= mesh.GetNOpenSegments(); i++)
-	{
-	  const Segment & seg = mesh.GetOpenSegment (i);
-	  if (seg.si == fnr)
-	    meshing.AddBoundaryElement (seg[0], seg[1], seg.geominfo[0], seg.geominfo[1]);
-	}
-      
-      
-      PrintMessage(3,"start meshing, trialcnt = ", retrynr);
 
-      /*
-      (*testout) << "start meshing with h = " << h << endl;
-      */
-      meshing.GenerateMesh (mesh, mparam, h, fnr);  // face index
 
-      extern void Render();
-      Render();
-    }    
+        for (int i = 1; i <= mesh.GetNOpenSegments(); i++)
+	  {
+	    const Segment & seg = mesh.GetOpenSegment (i);
+	    if (seg.si == fnr)
+	      meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], 
+					  seg.geominfo[0], seg.geominfo[1]);
+	  }
+
+	
+	PrintMessage(3,"start meshing, trialcnt = ", retrynr);
+	
+	/*
+	  (*testout) << "start meshing with h = " << h << endl;
+	*/
+	meshing.GenerateMesh (mesh, mparam, h, fnr);  // face index
+	
+	extern void Render();
+	Render();
+      }    
       
-  
+  NgProfiler::Print(stdout);
+
   mesh.CalcSurfacesOfNode();
 }
 
diff --git a/libsrc/stlgeom/stlgeom.cpp b/libsrc/stlgeom/stlgeom.cpp
index 9c537df6..cc789ec2 100644
--- a/libsrc/stlgeom/stlgeom.cpp
+++ b/libsrc/stlgeom/stlgeom.cpp
@@ -41,7 +41,7 @@ void STLMeshing (STLGeometry & geom,
     lineendpoints(), spiralpoints(), selectedmultiedge()
   */
 {
-	edgedata = new STLEdgeDataList(*this);
+  edgedata = new STLEdgeDataList(*this);
   externaledges.SetSize(0);
   Clear();
   meshchart = 0; // initialize all ?? JS
@@ -55,6 +55,7 @@ void STLMeshing (STLGeometry & geom,
   status = STL_GOOD;
   statustext = "Good Geometry";
   smoothedges = NULL;
+  area = -1;
 }
 
 STLGeometry :: ~STLGeometry()
@@ -2004,13 +2005,11 @@ void STLGeometry :: Clear()
 
 double STLGeometry :: Area()
 {
-  double ar = 0;
-  int i;
-  for (i = 1; i <= GetNT(); i++)
-    {
-      ar += GetTriangle(i).Area(points);
-    }
-  return ar;
+  if (area >= 0) return area;
+  area = 0;
+  for (int i = 1; i <= GetNT(); i++)
+    area += GetTriangle(i).Area(points);
+  return area;
 }
 
 double STLGeometry :: GetAngle(int t1, int t2)
@@ -3106,7 +3105,7 @@ int STLGeometry :: IsSmoothEdge (int pi1, int pi2) const
 
 
 
-
+/*
 //function is not used now
 int IsInArray(int n, const Array<int>& ia)
 {
@@ -3117,6 +3116,7 @@ int IsInArray(int n, const Array<int>& ia)
     }
   return 0;
 }
+*/
 
 void STLGeometry :: AddConeAndSpiralEdges()
 {
diff --git a/libsrc/stlgeom/stlgeom.hpp b/libsrc/stlgeom/stlgeom.hpp
index 05c7b014..85ae157f 100644
--- a/libsrc/stlgeom/stlgeom.hpp
+++ b/libsrc/stlgeom/stlgeom.hpp
@@ -26,8 +26,17 @@
 
 namespace netgen
 {
-  extern int IsInArray(int n, const Array<int>& ia);
-  extern int AddIfNotExists(Array<int>& list, int x);
+  inline int IsInArray(int n, const Array<int>& ia)
+  {
+    return ia.Contains(n); 
+  }
+
+  inline bool AddIfNotExists(Array<int>& list, int x)
+  {
+    if (list.Contains(x)) return false;
+    list.Append(x);
+    return true;
+  }
   
   extern DLL_HEADER MeshingParameters mparam;
   
@@ -166,6 +175,7 @@ namespace netgen
     Array<STLLine*> meshlines;
     Array<Point3d> meshpoints;
 
+    double area;
   public:
     STLGeometry();
     virtual ~STLGeometry();
diff --git a/libsrc/stlgeom/stlgeomchart.cpp b/libsrc/stlgeom/stlgeomchart.cpp
index 867c0cc2..47d2398b 100644
--- a/libsrc/stlgeom/stlgeomchart.cpp
+++ b/libsrc/stlgeom/stlgeomchart.cpp
@@ -22,7 +22,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
   int timer1 = NgProfiler::CreateTimer ("makeatlas");
   int timer2 = NgProfiler::CreateTimer ("makeatlas - part 2");
   int timer3 = NgProfiler::CreateTimer ("makeatlas - part 3");
-  int timer4 = NgProfiler::CreateTimer ("makeatlas - part 4");
 
   PushStatusF("Make Atlas");
 
@@ -33,7 +32,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 
   //speedup for make atlas
   if (GetNT() > 50000)
-    mesh.SetGlobalH(0.05*Dist (boundingbox.PMin(), boundingbox.PMax()));
+    mesh.SetGlobalH(min2 (0.05*Dist (boundingbox.PMin(), boundingbox.PMax()),
+			  mparam.maxh));
 
   atlas.SetSize(0);
   ClearSpiralPoints();
@@ -53,15 +53,15 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
   double sinchartangle = sin(chartangle);
   double sinouterchartangle = sin(outerchartangle);
 
-  Array<int> outermark(GetNT()); //marks all trigs form actual outer region
-  Array<int> outertested(GetNT()); //marks tested trigs for outer region
+  Array<int> outermark(GetNT());     //marks all trigs form actual outer region
+  Array<int> outertested(GetNT());   //marks tested trigs for outer region
   Array<int> pointstochart(GetNP()); //point in chart becomes chartnum
   Array<int> innerpointstochart(GetNP()); //point in chart becomes chartnum
-  Array<int> chartpoints; //point in chart becomes chartnum
+  Array<int> chartpoints;             //point in chart becomes chartnum
   Array<int> innerchartpoints;
   Array<int> dirtycharttrigs;
 
-  Array<int> chartdistacttrigs (GetNT()); //outercharttrigs
+  Array<int> chartdistacttrigs (GetNT());   //outercharttrigs
   chartdistacttrigs = 0;
 
 
@@ -71,8 +71,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 
   Array<int> chartpointchecked(GetNP()); //for dirty-chart-trigs
 
-  outermark.SetSize(GetNT());
-  outertested.SetSize(GetNT());
   pointstochart.SetSize(GetNP());
   innerpointstochart.SetSize(GetNP());
   chartmark.SetSize(GetNT());
@@ -117,21 +115,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 
       //find unmarked trig
       int prelastunmarked = lastunmarked;
-      /*
-      int j = lastunmarked;
-      bool found = 0;
-      while (j <= GetNT())
-	{
-	  if (!GetMarker(j)) 
-	    {
-	      found = 1; 
-	      lastunmarked = j;
-	      break;
-	    }
-	  else 
-	    j++;
-	}
-      */
+
       bool found = false;
       for (int j = lastunmarked; j <= GetNT(); j++)
 	if (!GetMarker(j))
@@ -255,15 +239,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 				  if (!accepted) {break;}
 				}
 			      
-			      /*
-				mindist = 1E50;
-				for (int ii = 1; ii <= 3; ii++)
-				{
-				tdist = Dist(GetPoint(GetTriangle(nt).PNum(ii)),startp);
-				if (tdist < mindist) {mindist = tdist;}
-				}
-				if (mindist > maxdist1) {accepted = 0;}
-			      */
 
 			      if (accepted)
 				{
@@ -337,16 +312,7 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 		  outertested.Elem(nt) = chartnum;
 		  
 		  Vec<3> n2 = GetTriangle(nt).Normal();
-		  /*
-		    double ang;
-		    ang = Angle(n2,sn);
-		    if (ang < -M_PI*0.5) {ang += 2*M_PI;}
-		    
-		    (*testout) << "ang < ocharang = " << (fabs(ang) <= outerchartangle);
-		    (*testout) << " = " << ( (n2 * sn) >= cosouterchartangle) << endl;
-		    
-		    //			      if (fabs(ang) <= outerchartangle) 
-		    */
+
 		  //abfragen, ob noch im tolerierten Winkel
 		  if ( (n2 * sn) >= cosouterchartangle )
 		    {
@@ -361,8 +327,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 		      
 		      //zurueckweisen, falls eine Spiralartige outerchart entsteht
 		      
-		      // NgProfiler::StartTimer (timer4);
-
 		      //find overlapping charts exacter: 
 		      //do not check dirty trigs!
 		      if (spiralcheckon && !isdirtytrig)
@@ -388,7 +352,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 			    if (!accepted) break;
 			  }
 		      
-		      //  NgProfiler::StopTimer (timer4);
 		      
 		      // outer chart is only small environment of
 		      //    inner chart:
@@ -409,7 +372,6 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 			      {
 				Point<3> pt = GetPoint(ntrig.PNum(k));					  
 				double h2 = sqr(mesh.GetH(pt));
-				
 				for (int l = 1; l <= innerchartpoints.Size(); l++)
 				  {
 				    double tdist = 
@@ -477,8 +439,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 	  lastunmarked = prelastunmarked;
 	}
 
-      //calculate an estimate meshsize, not to produce to large outercharts, with factor 2 larger!
+      //calculate an estimate meshsize, not to produce too large outercharts, with factor 2 larger!
       RestrictHChartDistOneChart(chartnum, chartdistacttrigs, mesh, h, 0.5, atlasminh);
+      // NgProfiler::Print(stdout);
     }
   
   NgProfiler::StopTimer (timer1);
@@ -610,6 +573,7 @@ int STLGeometry :: AtlasMade() const
 }
 
 
+/*
 //return 1 if not exists
 int AddIfNotExists(Array<int>& list, int x)
 {
@@ -621,6 +585,7 @@ int AddIfNotExists(Array<int>& list, int x)
   list.Append(x);
   return 1;
 }
+*/
 
 void STLGeometry :: GetInnerChartLimes(Array<twoint>& limes, int chartnum)
 {
diff --git a/libsrc/stlgeom/stlgeommesh.cpp b/libsrc/stlgeom/stlgeommesh.cpp
index c20c6acf..68c226cd 100644
--- a/libsrc/stlgeom/stlgeommesh.cpp
+++ b/libsrc/stlgeom/stlgeommesh.cpp
@@ -297,16 +297,13 @@ void STLGeometry :: PrepareSurfaceMeshing()
 {
   meshchart = -1; //clear no old chart
   meshcharttrigs.SetSize(GetNT());
-  int i;
-  for (i = 1; i <= GetNT(); i++) 
-    {meshcharttrigs.Elem(i) = 0;}
+  meshcharttrigs = 0;
 }
 
 void STLGeometry::GetMeshChartBoundary (Array<Point2d > & apoints,
 					Array<Point3d > & points3d,
 					Array<INDEX_2> & alines, double h)
 {
-  int i, j;
   twoint seg, newseg;
   int zone;
   Point<2> p2;
@@ -314,11 +311,11 @@ void STLGeometry::GetMeshChartBoundary (Array<Point2d > & apoints,
   const STLChart& chart = GetChart(meshchart);
 
 
-  for (i = 1; i <= chart.GetNOLimit(); i++)
+  for (int i = 1; i <= chart.GetNOLimit(); i++)
     {
       seg = chart.GetOLimit(i);
       INDEX_2 i2;
-      for (j = 1; j <= 2; j++)
+      for (int j = 1; j <= 2; j++)
 	{
 	  int pi = (j == 1) ? seg.i1 : seg.i2;
 	  int lpi;
@@ -358,7 +355,7 @@ void STLGeometry::GetMeshChartBoundary (Array<Point2d > & apoints,
       */
     }
 
-  for (i = 1; i <= chart.GetNOLimit(); i++)
+  for (int i = 1; i <= chart.GetNOLimit(); i++)
     {
       seg = chart.GetOLimit(i);
       ha_points.Elem(seg.i1) = 0;
@@ -1071,12 +1068,9 @@ void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh)
 
       //berechne minimale distanz von chart zu einem nicht-outerchart-punkt in jedem randpunkt einer chart
       
-      Array<int> acttrigs; //outercharttrigs
-      acttrigs.SetSize(GetNT());
-      for (i = 1; i <= GetNT(); i++)
-	{
-	  acttrigs.Elem(i) = 0;
-	}
+      Array<int> acttrigs(GetNT()); //outercharttrigs
+      acttrigs = 0;
+
       for (i = 1; i <= GetNOCharts(); i++)
 	{
 	  SetThreadPercent((double)i/(double)GetNOCharts()*100.);
@@ -1087,6 +1081,7 @@ void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh)
 	}
       
       PopStatus();
+      NgProfiler::Print(stdout);
     }
 
   if (stlparam.resthlinelengthenable)
@@ -1124,13 +1119,14 @@ void STLGeometry :: RestrictLocalH(class Mesh & mesh, double gh)
 void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrigs, 
 					       class Mesh & mesh, double gh, double fact, double minh)
 {
-  int i = chartnum;
-  int j;
+  static int timer1 = NgProfiler::CreateTimer ("restrictH OneChart 1");
+  static int timer2 = NgProfiler::CreateTimer ("restrictH OneChart 2");
+  static int timer3 = NgProfiler::CreateTimer ("restrictH OneChart 3");
 
+  NgProfiler::StartTimer (timer1);
   double limessafety = stlparam.resthchartdistfac*fact;  // original: 2
   double localh;
 
-  double f1,f2;
   //  mincalch = 1E10;
   //maxcalch = -1E10;  
   Array<int> limes1;
@@ -1146,8 +1142,8 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
 
   int divisions = 10;
 	  
-  int k, t, nt, np1, np2;
-  Point3d p3p1, p3p2;
+  int np1, np2;
+  // Point3d p3p1, p3p2;
   STLTriangle tt;
       
   limes1.SetSize(0);
@@ -1158,24 +1154,24 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
   plimes2trigs.SetSize(0);
   plimes1origin.SetSize(0);
 
-  STLChart& chart = GetChart(i);
+  STLChart& chart = GetChart(chartnum);
   chart.ClearOLimit();
   chart.ClearILimit();
 
-  for (j = 1; j <= chart.GetNChartT(); j++)
+  for (int j = 1; j <= chart.GetNChartT(); j++)
     {
-      t = chart.GetChartTrig(j); 
+      int t = chart.GetChartTrig(j); 
       tt = GetTriangle(t);
-      for (k = 1; k <= 3; k++)
+      for (int k = 1; k <= 3; k++)
 	{
-	  nt = NeighbourTrig(t,k); 
-	  if (GetChartNr(nt) != i)
+	  int nt = NeighbourTrig(t,k); 
+	  if (GetChartNr(nt) != chartnum)
 	    {	      
 	      tt.GetNeighbourPoints(GetTriangle(nt),np1,np2);
 	      if (!IsEdge(np1,np2) && !GetSpiralPoint(np1) && !GetSpiralPoint(np2))
 		{
-		  p3p1 = GetPoint(np1);
-		  p3p2 = GetPoint(np2);
+		  Point3d p3p1 = GetPoint(np1);
+		  Point3d p3p2 = GetPoint(np2);
 		  if (AddIfNotExists(limes1,np1)) 
 		    {
 		      plimes1.Append(p3p1); 
@@ -1192,8 +1188,8 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
 
 		  for (int di = 1; di <= divisions; di++)
 		    {
-		      f1 = (double)di/(double)(divisions+1.);
-		      f2 = (divisions+1.-(double)di)/(double)(divisions+1.);
+		      double f1 = (double)di/(double)(divisions+1.);
+		      double f2 = (divisions+1.-(double)di)/(double)(divisions+1.);
 			      
 		      plimes1.Append(Point3d(p3p1.X()*f1+p3p2.X()*f2,
 					     p3p1.Y()*f1+p3p2.Y()*f2,
@@ -1206,28 +1202,27 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
 	}
     }
 	  
-	 
-  for (j = 1; j <= chart.GetNT(); j++)
-    {
-      acttrigs.Elem(chart.GetTrig(j)) = i;
-    }
-	  
-  for (j = 1; j <= chart.GetNOuterT(); j++)
-    {
-      t = chart.GetOuterTrig(j); 
-      tt = GetTriangle(t);
-      for (k = 1; k <= 3; k++)
-	{
-	  nt = NeighbourTrig(t,k);
+  NgProfiler::StopTimer (timer1);	 
 
-	  if (acttrigs.Get(nt) != i)
+  NgProfiler::StartTimer (timer2);	 
+  for (int j = 1; j <= chart.GetNT(); j++)
+    acttrigs.Elem(chart.GetTrig(j)) = chartnum;
+
+  for (int j = 1; j <= chart.GetNOuterT(); j++)
+    {
+      int t = chart.GetOuterTrig(j); 
+      tt = GetTriangle(t);
+      for (int k = 1; k <= 3; k++)
+	{
+	  int nt = NeighbourTrig(t,k);
+	  if (acttrigs.Get(nt) != chartnum)
 	    {
 	      tt.GetNeighbourPoints(GetTriangle(nt),np1,np2);
 		      
 	      if (!IsEdge(np1,np2))
 		{
-		  p3p1 = GetPoint(np1);
-		  p3p2 = GetPoint(np2);
+		  Point3d p3p1 = GetPoint(np1);
+		  Point3d p3p2 = GetPoint(np2);
 			  
 		  if (AddIfNotExists(limes2,np1)) {plimes2.Append(p3p1); plimes2trigs.Append(t);}
 		  if (AddIfNotExists(limes2,np2)) {plimes2.Append(p3p2); plimes2trigs.Append(t);}
@@ -1235,8 +1230,8 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
 
 		  for (int di = 1; di <= divisions; di++)
 		    {
-		      f1 = (double)di/(double)(divisions+1.);
-		      f2 = (divisions+1.-(double)di)/(double)(divisions+1.);
+		      double f1 = (double)di/(double)(divisions+1.);
+		      double f2 = (divisions+1.-(double)di)/(double)(divisions+1.);
 			      
 		      plimes2.Append(Point3d(p3p1.X()*f1+p3p2.X()*f2,
 					     p3p1.Y()*f1+p3p2.Y()*f2,
@@ -1248,24 +1243,25 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
 	}
     }
 	  
-	  
+  NgProfiler::StopTimer (timer2);	 
+  NgProfiler::StartTimer (timer3);	 	  
+
   double chartmindist = 1E50;
 
   if (plimes2.Size())
     {
       Box3d bbox;
       bbox.SetPoint (plimes2.Get(1));
-      for (j = 2; j <= plimes2.Size(); j++)
+      for (int j = 2; j <= plimes2.Size(); j++)
 	bbox.AddPoint (plimes2.Get(j));
       Point3dTree stree(bbox.PMin(), bbox.PMax());
-      for (j = 1; j <= plimes2.Size(); j++)
+      for (int j = 1; j <= plimes2.Size(); j++)
 	stree.Insert (plimes2.Get(j), j);
       Array<int> foundpts;
 	  
-      for (j = 1; j <= plimes1.Size(); j++)
+      for (int j = 1; j <= plimes1.Size(); j++)
 	{
 	  double mindist = 1E50;
-	  double dist;
 
 	  const Point3d & ap1 = plimes1.Get(j);
 	  double boxs = mesh.GetH (plimes1.Get(j)) * limessafety;
@@ -1278,12 +1274,9 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
 
 	  for (int kk = 1; kk <= foundpts.Size(); kk++)
 	    {
-	      k = foundpts.Get(kk);
-	      dist = Dist2(plimes1.Get(j),plimes2.Get(k));
-	      if (dist < mindist) 
-		{
-		  mindist = dist;
-		}
+	      int k = foundpts.Get(kk);
+	      double dist = Dist2(plimes1.Get(j),plimes2.Get(k));
+	      if (dist < mindist) mindist = dist;
 	    }
 
 	  /*
@@ -1325,7 +1318,7 @@ void STLGeometry :: RestrictHChartDistOneChart(int chartnum, Array<int>& acttrig
 	    }
 	}
     }
-
+  NgProfiler::StopTimer (timer3);	 
 }
 
 
diff --git a/libsrc/stlgeom/stlline.cpp b/libsrc/stlgeom/stlline.cpp
index 2c836f4e..9a49b748 100644
--- a/libsrc/stlgeom/stlline.cpp
+++ b/libsrc/stlgeom/stlline.cpp
@@ -648,6 +648,13 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
 			 Array<Point3d>& mp, double ghi,
 			 class Mesh& mesh) const
 {
+  static int timer1a = NgProfiler::CreateTimer ("mesh stl-line 1a");
+  static int timer1b = NgProfiler::CreateTimer ("mesh stl-line 1b");
+  static int timer2 = NgProfiler::CreateTimer ("mesh stl-line 2");
+  static int timer3 = NgProfiler::CreateTimer ("mesh stl-line 3");
+
+  NgProfiler::StartTimer (timer1a);
+
   STLLine* line = new STLLine(geometry);
 
   //stlgh = ghi; //uebergangsloesung!!!!
@@ -659,8 +666,6 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
   int ind;
   Point3d p;
 
-  int i, j;
-
   Box<3> bbox;
   GetBoundingBox (ap, bbox);
   double diam = bbox.Diam();
@@ -668,7 +673,7 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
   double minh = mesh.LocalHFunction().GetMinH (bbox.PMin(), bbox.PMax());
 
   double maxseglen = 0;
-  for (i = 1; i <= GetNS(); i++)
+  for (int i = 1; i <= GetNS(); i++)
     maxseglen = max2 (maxseglen, GetSegLen (ap, i));
   
   int nph = 10+int(maxseglen / minh); //anzahl der integralauswertungen pro segment
@@ -676,11 +681,14 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
   Array<double> inthi(GetNS()*nph);
   Array<double> curvelen(GetNS()*nph);
 
+  NgProfiler::StopTimer (timer1a);
+  NgProfiler::StartTimer (timer1b);
 
-  for (i = 1; i <= GetNS(); i++)
+
+  for (int i = 1; i <= GetNS(); i++)
     {
       //double seglen = GetSegLen(ap,i);
-      for (j = 1; j <= nph; j++)
+      for (int j = 1; j <= nph; j++)
 	{
 	  p = GetPointInDist(ap,dist,ind);
 	  //h = GetH(p,dist/len);
@@ -709,7 +717,7 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
      
   double fact = inthl/(double)inthlint;
   dist = 0;
-  j = 1;
+  int j = 1;
 
 
   p = ap.Get(StartP());
@@ -721,18 +729,20 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
   line->AddRightTrig(GetRightTrig(segn));
   line->AddDist(dist);
 
+  NgProfiler::StopTimer (timer1b);
+  NgProfiler::StartTimer (timer2);
+
   inthl = 0; //restart each meshseg
-  for (i = 1; i <= inthlint; i++)
+  for (int i = 1; i <= inthlint; i++)
     {
       while (inthl < 1.000000001 && j <= inthi.Size())
-      //      while (inthl-1. < 1e-9) && j <= inthi.Size())
 	{
 	  inthl += inthi.Get(j)/fact;
 	  dist += curvelen.Get(j);
 	  j++;
 	}
 
-      //went to far:
+      //went too far:
       j--;
       double tofar = (inthl - 1)/inthi.Get(j);
       inthl -= tofar*inthi.Get(j);
@@ -759,6 +769,10 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
       j++;
     }
 
+  NgProfiler::StopTimer (timer2);
+  NgProfiler::StartTimer (timer3);
+
+
   p = ap.Get(EndP());
   pn = AddPointIfNotExists(mp, p, 1e-10*diam);
   segn = GetNS();
@@ -776,6 +790,9 @@ STLLine* STLLine :: Mesh(const Array<Point<3> >& ap,
   (*testout) << "line, " << ap.Get(StartP()) << "-" << ap.Get(EndP())
 	     << " len = " << Dist (ap.Get(StartP()), ap.Get(EndP())) << endl;
   */
+
+  NgProfiler::StopTimer (timer3);
+
   return line;
 }
 }
diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp
index 80ed764b..20cb2c61 100644
--- a/libsrc/stlgeom/stltool.cpp
+++ b/libsrc/stlgeom/stltool.cpp
@@ -15,8 +15,9 @@ namespace netgen
 //add a point into a pointlist, return pointnumber
 int AddPointIfNotExists(Array<Point3d>& ap, const Point3d& p, double eps)
 {
+  double eps2 = sqr(eps);
   for (int i = 1; i <= ap.Size(); i++)
-    if (Dist(ap.Get(i),p) <= eps ) 
+    if (Dist2(ap.Get(i),p) <= eps2 ) 
       return i;
   return ap.Append(p);
 }
diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp
index 9fa6ae7f..fd9ab2e6 100644
--- a/libsrc/visualization/vsmesh.cpp
+++ b/libsrc/visualization/vsmesh.cpp
@@ -219,8 +219,6 @@ namespace netgen
 
 #ifdef PARALLELGL
 	if (ntasks > 1 && vispar.drawtetsdomain > 0 && vispar.drawtetsdomain < ntasks)
-	  // for (int dest = 1; dest < ntasks; dest++)
-	  // if (vispar.drawtetsdomain == dest)
 	  glCallList (par_linelists[vispar.drawtetsdomain]);
 	else
 #endif
@@ -230,7 +228,6 @@ namespace netgen
 	glDisable (GL_POLYGON_OFFSET_LINE);
       }
 
-
     if (vispar.drawidentified)
       {
 	glPolygonOffset (1, -1);