mesh optimization improvements

libsrc/csg/algprim.cpp

@@ -351,6 +351,7 @@ namespace netgen
   {
     c = ac;
     r = ar;
+    invr = 1.0/r;
 
     cxx = cyy = czz = 0.5 / r;
     cxy = cxz = cyz = 0;
@@ -377,6 +378,7 @@ namespace netgen
     c(2) = coeffs.Elem(3);
     r = coeffs.Elem(4);
 
+    invr = 1.0/r;
     cxx = cyy = czz = 0.5 / r;
     cxy = cxz = cyz = 0;
     cx = - c(0) / r;
@@ -412,6 +414,10 @@ namespace netgen
   }
 
 
+  double Sphere :: CalcFunctionValue (const Point<3> & point) const
+  {
+    return 0.5* (invr * Abs2 (point-c) - r);
+  }
 
 
   int Sphere :: IsIdentic (const Surface & s2, int & inv, double eps) const

libsrc/csg/algprim.hpp

@@ -122,7 +122,7 @@ namespace netgen
     ///
     Point<3> c;
     ///
-    double r;
+    double r, invr;
   public:
     ///
     Sphere (const Point<3> & ac, double ar);
@@ -135,6 +135,8 @@ namespace netgen
     virtual Primitive * Copy () const;
     virtual void Transform (Transformation<3> & trans);
 
+    virtual double CalcFunctionValue (const Point<3> & point) const;
+
 
     virtual int IsIdentic (const Surface & s2, int & inv, double eps) const;
 

libsrc/csg/genmesh.cpp

@@ -221,16 +221,6 @@ namespace netgen
 
   static void MeshSurface (CSGeometry & geom, Mesh & mesh)
   {
-    /*
-    Point3d pmin, pmax;
-    mesh.GetBox(pmin, pmax);
-    cout << "box = " << pmin << " -  " << pmax << endl;
-    cout << "localhbox = " 
-         << mesh.LocalHFunction().GetBoundingBox().PMin() << " - " 
-         << mesh.LocalHFunction().GetBoundingBox().PMax() 
-         << endl;
-    */
-
     const char * savetask = multithread.task;
     multithread.task = "Surface meshing";
   
@@ -429,7 +419,7 @@ namespace netgen
         double eps = 1e-8 * geom.MaxSize();
 	for (PointIndex pi = PointIndex::BASE; pi < noldp+PointIndex::BASE; pi++)
 	  { 
-	    //if(surf->PointOnSurface(mesh[pi]))
+	    // if(surf->PointOnSurface(mesh[pi]))
 	    meshing.AddPoint (mesh[pi], pi, NULL,
 			      (surf->PointOnSurface(mesh[pi], eps) != 0));
 	  }
@@ -511,6 +501,11 @@ namespace netgen
 	if (segments.Size())   
 	  { 
 	    // surface was meshed, not copied
+
+	    static int timer = NgProfiler::CreateTimer ("total surface mesh optimization");
+	    NgProfiler::RegionTimer reg (timer);
+
+
 	    PrintMessage (2, "Optimize Surface");
 	    for (int i = 1; i <= mparam.optsteps2d; i++)
 	      {

libsrc/csg/solid.cpp

@@ -689,7 +689,7 @@ namespace netgen
       case UNION:
 	{
 	  int in1, in2, strin1, strin2;
-	  Solid * tansol1, * tansol2;
+	  Solid * tansol1 = 0, * tansol2 = 0;
 
 	  s1 -> RecTangentialSolid (p, tansol1, surfids, in1, strin1, eps);
 	  s2 -> RecTangentialSolid (p, tansol2, surfids, in2, strin2, eps);
@@ -703,6 +703,11 @@ namespace netgen
 	      else if (tansol2)
 		tansol = tansol2;
 	    }
+	  else
+	    {
+	      delete tansol1;
+	      delete tansol2;
+	    }
 	  in = (in1 || in2);
 	  strin = (strin1 || strin2);
 	  break;

libsrc/csg/specpoin.cpp

@@ -76,6 +76,10 @@ namespace netgen
   CalcSpecialPoints (const CSGeometry & ageometry, 
 		     Array<MeshPoint> & apoints)
   {
+    static int timer = NgProfiler::CreateTimer ("CSG: find special points");
+    NgProfiler::RegionTimer reg (timer);
+
+
     geometry = &ageometry;
     points = &apoints;
 
@@ -231,7 +235,7 @@ namespace netgen
     // explicit solution for planes only and at most one quadratic
     if (numprim <= 5)
       {
-	int nplane = 0, nquad = 0, quadi = -1;
+	int nplane = 0, nquad = 0, quadi = -1, nsphere = 0;
 	const QuadraticSurface *qsurf = 0, *qsurfi;
 
 	for (int i = 0; i < numprim; i++)
@@ -247,6 +251,9 @@ namespace netgen
 		quadi = i;
 		qsurf = qsurfi;
 	      }
+
+	    if (dynamic_cast<const Sphere*> (qsurfi))
+	      nsphere++;
 	  }
 
 	/*
@@ -374,6 +381,82 @@ namespace netgen
 	    
 	    return;
 	  }
+
+
+
+	if (nsphere == numprim) //  && calccp == false)
+	  {
+	    Array<Point<3> > pts;
+	    Array<int> surfids;
+
+
+	    
+	    for (int k1 = 0; k1 < numprim; k1++)
+	      for (int k2 = 0; k2 < k1; k2++)
+		for (int k3 = 0; k3 < k2; k3++)
+		  {
+		    ComputeCrossPoints (dynamic_cast<const Sphere*> (geometry->GetSurface(locsurf[k1])),
+					dynamic_cast<const Sphere*> (geometry->GetSurface(locsurf[k2])),
+					dynamic_cast<const Sphere*> (geometry->GetSurface(locsurf[k3])),
+					pts);
+		    
+		    for (int j = 0; j < pts.Size(); j++)
+		      if (Dist (pts[j], box.Center()) < box.Diam()/2)
+			{
+			  Solid * tansol;
+			  sol -> TangentialSolid (pts[j], tansol, surfids, 1e-9*size);
+			  
+			  if(!tansol)
+			    continue;
+			  
+			  bool ok1 = false, ok2 = false, ok3 = false;
+			  int rep1 = geometry->GetSurfaceClassRepresentant(locsurf[k1]);
+			  int rep2 = geometry->GetSurfaceClassRepresentant(locsurf[k2]);
+			  int rep3 = geometry->GetSurfaceClassRepresentant(locsurf[k3]);
+			  for(int jj=0; jj<surfids.Size(); jj++)
+			    {
+			      int actrep = geometry->GetSurfaceClassRepresentant(surfids[jj]);
+			      if(actrep == rep1) ok1 = true;
+			      if(actrep == rep2) ok2 = true;
+			      if(actrep == rep3) ok3 = true;				  
+			    }
+			  
+			  
+			  if (tansol && ok1 && ok2 && ok3)
+			    // if (sol -> IsIn (pts[j], 1e-6*size) && !sol->IsStrictIn (pts[j], 1e-6*size))
+			    {
+			      if (AddPoint (pts[j], layer))
+				(*testout) << "cross point found, 1: " << pts[j] << endl;
+			    }  
+			  delete tansol;
+			}
+		  }
+	    
+
+	    for (int k1 = 0; k1 < numprim; k1++)
+	      for (int k2 = 0; k2 < k1; k2++)
+		{
+		  ComputeExtremalPoints (dynamic_cast<const Sphere*> (geometry->GetSurface(locsurf[k1])),
+					 dynamic_cast<const Sphere*> (geometry->GetSurface(locsurf[k2])),
+					 pts);
+		  
+		  for (int j = 0; j < pts.Size(); j++)
+		    if (Dist (pts[j], box.Center()) < box.Diam()/2)
+		      {
+			Solid * tansol;
+			sol -> TangentialSolid (pts[j], tansol, surfids, 1e-9*size);
+			if (tansol)
+			  // sol -> IsIn (pts[j], 1e-6*size) && !sol->IsStrictIn (pts[j], 1e-6*size) )
+			  {
+			    if (AddPoint (pts[j], layer))
+			      (*testout) << "extremal point found, spheres: " << pts[j] << endl;
+			  }  
+			delete tansol;
+		      }
+		}
+	    
+	    return;
+	  }
       }
 
 
@@ -1143,6 +1226,83 @@ namespace netgen
 
 
 
+  void SpecialPointCalculation :: 
+  ComputeCrossPoints (const Sphere * sphere1, 
+		      const Sphere * sphere2, 
+		      const Sphere * sphere3, 
+		      Array<Point<3> > & pts)
+  {
+    Mat<2,3> mat;
+    Mat<3,2> inv;
+    Vec<2> rhs;
+    Vec<3> sol, t;
+    Point<3> p0(0,0,0);
+
+    pts.SetSize (0);
+
+
+    Point<3> c1 = sphere1 -> Center();
+    Point<3> c2 = sphere2 -> Center();
+    Point<3> c3 = sphere3 -> Center();
+    double r1 = sphere1 -> Radius();
+    double r2 = sphere2 -> Radius();
+    double r3 = sphere3 -> Radius();
+
+
+    Vec<3> a1 = c2-c1;
+    double b1 = 0.5 * (sqr(r1) - sqr(r2) - Abs2(Vec<3> (c1)) + Abs2(Vec<3> (c2)) );
+
+    Vec<3> a2 = c3-c1;
+    double b2 = 0.5 * (sqr(r1) - sqr(r3) - Abs2(Vec<3> (c1)) + Abs2(Vec<3> (c3)) );
+
+
+    for (int j = 0; j < 3; j++)
+      {
+	mat(0,j) = a1(j);
+	mat(1,j) = a2(j);
+      }
+    
+    rhs(0) = b1;
+    rhs(1) = b2;
+
+
+    CalcInverse (mat, inv);
+    sol = inv * rhs;
+    t = Cross (mat.Row(0), mat.Row(1));
+
+    if (t.Length() > 1e-8)
+      {
+	Point<3> p (sol);
+	// quadratic on  p + s t = 0
+	double quad_a;
+	Vec<3> quad_b;
+	Mat<3> quad_c;
+	
+	quad_a = sphere1 -> CalcFunctionValue(p);
+	sphere1 -> CalcGradient (p, quad_b);
+	sphere1 -> CalcHesse (p, quad_c);
+	
+	double a, b, c;
+	a = quad_a;
+	b = quad_b * t;
+	c = 0.5 * t * (quad_c * t);
+
+	// a  + s b + s^2 c = 0;
+	double disc = b*b-4*a*c;
+	if (disc > 1e-10 * fabs (b))
+	  {
+	    disc = sqrt (disc);
+	    double s1 = (-b-disc) / (2*c);
+	    double s2 = (-b+disc) / (2*c);
+
+	    pts.Append (p + s1 * t);
+	    pts.Append (p + s2 * t);
+	  }
+      }
+  }
+
+
+
 
 
 
@@ -1238,6 +1398,138 @@ namespace netgen
 
 
 
+
+  void SpecialPointCalculation :: 
+  ComputeExtremalPoints (const Sphere * sphere1,
+			 const Sphere * sphere2,
+			 Array<Point<3> > & pts)
+  {
+    // 3 equations:
+    // surf1 = 0  <===> |x-c1|^2 - r1^2 = 0;
+    // surf2 = 0  <===> |x-c2|^2 - r2^2 = 0;
+    // (grad 1 x grad 2)(i) = 0  <====> (x-p1) x (p1-p2) . e_i = 0;
+
+    pts.SetSize (0);
+
+    Point<3> c1 = sphere1 -> Center();
+    Point<3> c2 = sphere2 -> Center();
+    double r1 = sphere1 -> Radius();
+    double r2 = sphere2 -> Radius();
+
+    /*
+    *testout << "\n\ncompute extremalpoint, sphere-sphere" << endl;
+    *testout << "c1 = " << c1 << ", r1 = " << r1 << endl;
+    *testout << "c2 = " << c2 << ", r2 = " << r2 << endl;
+    *testout << "dist = " << Abs (c2-c1) << ", r1+r2 = " << r1+r2 << endl;
+    */
+
+    Vec<3> v12 = c2 - c1;
+
+    Vec<3> a1, a2;
+    double b1, b2;
+
+    // eqn: ai . x = bi
+
+    a1 = v12;
+    b1 = 0.5 * (sqr(r1) - sqr(r2) - Abs2(Vec<3> (c1)) + Abs2(Vec<3> (c2)) );
+
+    int dir = 0;
+    for (int j = 1; j < 3; j++)
+      if (fabs (v12(j)) > v12(dir))
+	dir = j;
+    
+    // *testout << "dir = " << dir << endl;
+
+    Vec<3> ei = 0.0;
+    ei(dir) = 1;
+    a2 = Cross (v12, ei);
+    b2 = Vec<3>(c1) * a2;
+    
+    
+    Point<3> p0 (0,0,0);
+    double quad_a;
+    Vec<3> quad_b;
+    Mat<3> quad_c;
+
+    quad_a = sphere1 -> CalcFunctionValue(p0);
+    sphere1 -> CalcGradient (p0, quad_b);
+    sphere1 -> CalcHesse (p0, quad_c);
+    for (int i = 0; i < 3; i++)
+      for (int j = 0; j < 3; j++)
+	quad_c(i,j) *= 0.5;
+
+    
+    // find line of two linear equations:
+    
+    Vec<2> rhs;
+    Vec<3> sol;
+    Mat<2,3> mat;
+    
+    for (int j = 0; j < 3; j++)
+      {
+	mat(0,j) = a1(j);
+	mat(1,j) = a2(j);
+      }
+    rhs(0) = b1;
+    rhs(1) = b2;
+
+
+    // *testout << "mat = " << endl << mat << endl;
+    // *testout << "rhs = " << endl << rhs << endl;
+
+    Vec<3> t = Cross (a1, a2);
+    if (Abs2(t) > 0)
+      {
+	mat.Solve (rhs, sol);
+	
+	/*
+	*testout << "sol = " << endl << sol << endl;
+
+	*testout << "a * sol = " << mat * sol << endl;
+
+	*testout << "c1-sol = " << Abs (Vec<3>(c1)-sol) << endl;
+	*testout << "c2-sol = " << Abs (Vec<3>(c2)-sol) << endl;
+	*/
+
+	// solve quadratic equation along line  sol + alpha t ....
+	double a = quad_a + quad_b * sol + sol * (quad_c * sol);
+	double b = quad_b * t + 2 * (sol * (quad_c * t));
+	double c = t * (quad_c * t);
+
+	// solve a + b alpha + c alpha^2:
+	
+	if (fabs (c) > 1e-32)
+	  {
+	    double disc = sqr (0.5*b/c) - a/c;
+	    if (disc > 0)
+	      {
+		disc = sqrt (disc);
+		double alpha1 = -0.5*b/c + disc;
+		double alpha2 = -0.5*b/c - disc;
+		
+		pts.Append (Point<3> (sol+alpha1*t));
+		pts.Append (Point<3> (sol+alpha2*t));
+
+		// *testout << "pts = " << endl << pts << endl;
+
+		/*
+		  cout << "sol1 = " << sol + alpha1 * t
+		  << ", sol2 = " << sol + alpha2 * t << endl;
+		*/
+	      }
+	  }
+      }
+  }
+
+
+
+
+
+
+
+
+
+
   /*
     bool SpecialPointCalculation :: ExtremalPointPossible (const Surface * f1, 
     const Surface * f2, 
@@ -1445,11 +1737,11 @@ namespace netgen
 	    if (tlo->GetLayer() != apoints[i].GetLayer())
 	      continue;
 	    
+
 	    Solid * locsol;
 	    sol -> TangentialSolid (p, locsol, surfind, ideps*geomsize);
 
 
-
 	    rep_surfind.SetSize (surfind.Size());
 	    int num_indep_surfs = 0;
 	    
@@ -1471,6 +1763,7 @@ namespace netgen
 
 	    if (!locsol) continue;
 
+	  
 	    // get all surface indices, 
 	    if (surf)
 	      {

libsrc/csg/specpoin.hpp

@@ -161,6 +161,11 @@ namespace netgen
 				const QuadraticSurface * quadric, 
 				Array<Point<3> > & pts);
 
+    void ComputeExtremalPoints (const Sphere * sphere1, 
+				const Sphere * sphere2, 
+				Array<Point<3> > & pts);
+
+
     void ComputeCrossPoints (const Plane * plane1, 
 			     const Plane * plane2, 
 			     const Plane * plane3, 
@@ -170,6 +175,11 @@ namespace netgen
 			     const Plane * plane2, 
 			     const QuadraticSurface * quadratic, 
 			     Array<Point<3> > & pts);
+
+    void ComputeCrossPoints (const Sphere * sphere1, 
+			     const Sphere * sphere2, 
+			     const Sphere * sphere3, 
+			     Array<Point<3> > & pts);
   };
 
 }

libsrc/meshing/adfront2.cpp

@@ -94,6 +94,7 @@ namespace netgen
     if (mgi)
       cpointsearchtree.Insert (p, pi);
 
+    if (pointonsurface)
       pointsearchtree.Insert (p, pi);
     
     return pi;
@@ -287,7 +288,9 @@ namespace netgen
 			     Array<INDEX> & lindex,
 			     double xh)
   {
-    // baselineindex += 1-lines.Begin();
+    static int timer = NgProfiler::CreateTimer ("adfront2::GetLocals");
+    NgProfiler::RegionTimer reg (timer);
+
 
     int pstind;
     Point<3>  midp, p0;
@@ -296,13 +299,12 @@ namespace netgen
     p0 = points[pstind].P();
 
     loclines.Append(lines[baselineindex].L());
-    lindex.Append(baselineindex);  //  +1-lines.Begin());
+    lindex.Append(baselineindex);  
 
-    static Array<int> nearlines;
+    ArrayMem<int, 1000> nearlines(0);
-    nearlines.SetSize(0);
+    ArrayMem<int, 1000> nearpoints(0);
-    static Array<int> nearpoints;
-    nearpoints.SetSize(0);
 
+    // dominating costs !!
     linesearchtree.GetIntersecting (p0 - Vec3d(xh, xh, xh),
 				    p0 + Vec3d(xh, xh, xh),
 				    nearlines);
@@ -311,11 +313,10 @@ namespace netgen
                                      p0 + Vec3d(xh, xh, xh),
                                      nearpoints);
     
-
+    for (int ii = 0; ii < nearlines.Size(); ii++)
-    for (int ii = 1; ii <= nearlines.Size(); ii++)
       {
-	int i = nearlines.Get(ii);
+	int i = nearlines[ii];
-	if (lines[i].Valid() && i != baselineindex) //  + 1-lines.Begin())
+	if (lines[i].Valid() && i != baselineindex) 
 	  {
             loclines.Append(lines[i].L());
             lindex.Append(i);

libsrc/meshing/improve2.cpp

@@ -12,42 +12,50 @@
 namespace netgen
 {
 
-class Neighbour
+  class Neighbour
-{
+  {
     int nr[3];
     int orient[3];
 
-public:
+  public:
-  Neighbour () { nr[0] = nr[1] = nr[2] = -1; orient[0] = orient[1] = orient[2] = 0; }
+    Neighbour () { ; } 
 
-  void SetNr (int side, int anr) { nr[side-1] = anr; }
+    void SetNr (int side, int anr) { nr[side] = anr; }
-  int GetNr (int side) { return nr[side-1]; }
+    int GetNr (int side) { return nr[side]; }
 
-  void SetOrientation (int side, int aorient) { orient[side-1] = aorient; }
+    void SetOrientation (int side, int aorient) { orient[side] = aorient; }
-  int GetOrientation (int side) { return orient[side-1]; }
+    int GetOrientation (int side) { return orient[side]; }
-};
+
+
+
+    /*
+      void SetNr1 (int side, int anr) { nr[side-1] = anr; }
+      int GetNr1 (int side) { return nr[side-1]; }
+
+      void SetOrientation1 (int side, int aorient) { orient[side-1] = aorient; }
+      int GetOrientation1 (int side) { return orient[side-1]; }
+    */
+  };
 
 
 
 
-class trionedge
+  class trionedge
-{
+  {
-public:
+  public:
     int tnr;
     int sidenr;
 
     trionedge () { tnr = 0; sidenr = 0; }
     trionedge (int atnr, int asidenr)
     { tnr = atnr; sidenr = asidenr; }
-};
+  };
 
 
 
  
-void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
+  void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
-{
+  {
-  //  return; 
-
     if (!faceindex)
       {
 	if (usemetric)
@@ -72,15 +80,14 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
     static int timer = NgProfiler::CreateTimer ("EdgeSwapping 2D");
     NgProfiler::RegionTimer reg1 (timer);
 
+    static int timerstart = NgProfiler::CreateTimer ("EdgeSwapping 2D start");
+    NgProfiler::StartTimer (timerstart);
 
-  int i, i2, j, j2;
-  bool should;
-  PointIndex pi;
 
     Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
 
-  for (i = 0; i < seia.Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
       if (mesh[seia[i]].GetNP() != 3)
 	{
 	  GenericImprove (mesh);
@@ -97,27 +104,27 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
     Array<int,PointIndex::BASE> pdef(mesh.GetNP());
     Array<double,PointIndex::BASE> pangle(mesh.GetNP());
 
-  SurfaceElementIndex t1, t2;
-  int o1, o2;
 
-  PointIndex pi1, pi2, pi3, pi4;
+    // int e;
-  PointGeomInfo gi1, gi2, gi3, gi4;
+    // double d;
+    // Vec3d nv1, nv2;
 
-
+    // double loch(-1);
-  int nswaps = 0;
-  int e, done;
-  double d;
-  Vec3d nv1, nv2;
-  double horder;
-  double loch(-1);
     static const double minangle[] = { 0, 1.481, 2.565, 3.627, 4.683, 5.736, 7, 9 };
 
-  pangle = 0;
 
-  for (i = 0; i < seia.Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
       {
 	const Element2d & sel = mesh[seia[i]];
-      for (j = 0; j < 3; j++)
+	for (int j = 0; j < 3; j++)
+	  pangle[sel[j]] = 0.0;
+      }
+    // pangle = 0;
+
+    for (int i = 0; i < seia.Size(); i++)
+      {
+	const Element2d & sel = mesh[seia[i]];
+	for (int j = 0; j < 3; j++)
 	  {
 	    POINTTYPE typ = mesh[sel[j]].Type();
 	    if (typ == FIXEDPOINT || typ == EDGEPOINT)
@@ -129,31 +136,38 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 	  }
       }
 
-  for (pi = PointIndex::BASE; 
+    // for (PointIndex pi = PointIndex::BASE; 
-       pi < mesh.GetNP()+PointIndex::BASE; pi++)
+    // pi < mesh.GetNP()+PointIndex::BASE; pi++)
+    
+    // pdef = 0;
+    for (int i = 0; i < seia.Size(); i++)
       {
+	const Element2d & sel = mesh[seia[i]];
+	for (int j = 0; j < 3; j++)
+	  {
+	    PointIndex pi = sel[j];
 	    if (mesh[pi].Type() == INNERPOINT || mesh[pi].Type() == SURFACEPOINT)
 	      pdef[pi] = -6;
 	    else
-	for (j = 0; j < 8; j++)
+	      for (int j = 0; j < 8; j++)
 		if (pangle[pi] >= minangle[j])
 		  pdef[pi] = -1-j;
 	  }
+      }
 
-  for (i = 0; i < seia.Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
       {
 	const Element2d & sel = mesh[seia[i]];
-      for (j = 0; j < 3; j++)
+	for (int j = 0; j < 3; j++)
 	  pdef[sel[j]]++;
       }
 
-  for (i = 0; i < seia.Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
       {
-      //const Element2d & sel = mesh[seia[i]];
+	for (int j = 0; j < 3; j++)
-      for (j = 0; j < 3; j++)
 	  {
-	  neighbors[seia[i]].SetNr (j+1, -1);
+	    neighbors[seia[i]].SetNr (j, -1);
-	  neighbors[seia[i]].SetOrientation (j+1, 0);
+	    neighbors[seia[i]].SetOrientation (j, 0);
 	  }
       }
 
@@ -174,17 +188,16 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
       }
     */	    
 
-  
+    for (int i = 0; i < seia.Size(); i++)
-  for (i = 0; i < seia.Size(); i++)
       {
 	const Element2d & sel = mesh[seia[i]];
 
-      for (j = 1; j <= 3; j++)
+	for (int j = 0; j < 3; j++)
 	  {
-	  pi1 = sel.PNumMod(j+1);
+	    PointIndex pi1 = sel.PNumMod(j+2);
-	  pi2 = sel.PNumMod(j+2);
+	    PointIndex pi2 = sel.PNumMod(j+3);
 	  
-	  loch = mesh.GetH(mesh[pi1]);
+	    //	    double loch = mesh.GetH(mesh[pi1]);
 	    
 	    INDEX_2 edge(pi1, pi2);
 	    edge.Sort();
@@ -202,8 +215,8 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 		// INDEX_2 i2s(ii2);
 		// i2s.Sort();
 	      
-	      i2 = other.Get(ii2).tnr;
+		int i2 = other.Get(ii2).tnr;
-	      j2 = other.Get(ii2).sidenr;
+		int j2 = other.Get(ii2).sidenr;
 		
 		neighbors[seia[i]].SetNr (j, i2);
 		neighbors[seia[i]].SetOrientation (j, j2);
@@ -217,17 +230,20 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 	  }
       }
 
-  for (i = 0; i < seia.Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
       swapped[seia[i]] = 0;
 
+    NgProfiler::StopTimer (timerstart);
+  
+
 
     int t = 4;
-  done = 0;
+    int done = 0;
     while (!done && t >= 2)
       {
-      for (i = 0; i < seia.Size(); i++)
+	for (int i = 0; i < seia.Size(); i++)
 	  {
-	  t1 = seia[i];
+	    SurfaceElementIndex t1 = seia[i];
 
 	    if (mesh[t1].IsDeleted())
 	      continue;
@@ -238,39 +254,40 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 	    if (multithread.terminate)
 	      throw NgException ("Meshing stopped");
 
-	  for (o1 = 1; o1 <= 3; o1++)
+	    for (int o1 = 0; o1 < 3; o1++)
 	      {
-	      t2 = neighbors[t1].GetNr (o1);
+		bool should;
-	      o2 = neighbors[t1].GetOrientation (o1);
+
+
+		SurfaceElementIndex t2 = neighbors[t1].GetNr (o1);
+		int o2 = neighbors[t1].GetOrientation (o1);
 
 		if (t2 == -1) continue;
 		if (swapped[t1] || swapped[t2]) continue;
 	      
 
-	      pi1 = mesh[t1].PNumMod(o1+1);
+		PointIndex pi1 = mesh[t1].PNumMod(o1+1+1);
-	      pi2 = mesh[t1].PNumMod(o1+2);
+		PointIndex pi2 = mesh[t1].PNumMod(o1+1+2);
-	      pi3 = mesh[t1].PNumMod(o1);
+		PointIndex pi3 = mesh[t1].PNumMod(o1+1);
-	      pi4 = mesh[t2].PNumMod(o2);
+		PointIndex pi4 = mesh[t2].PNumMod(o2+1);
 	      
-	      gi1 = mesh[t1].GeomInfoPiMod(o1+1);
+		PointGeomInfo gi1 = mesh[t1].GeomInfoPiMod(o1+1+1);
-	      gi2 = mesh[t1].GeomInfoPiMod(o1+2);
+		PointGeomInfo gi2 = mesh[t1].GeomInfoPiMod(o1+1+2);
-	      gi3 = mesh[t1].GeomInfoPiMod(o1);
+		PointGeomInfo gi3 = mesh[t1].GeomInfoPiMod(o1+1);
-	      gi4 = mesh[t2].GeomInfoPiMod(o2);
+		PointGeomInfo gi4 = mesh[t2].GeomInfoPiMod(o2+1);
 	    
 		bool allowswap = true;
 
-	      Vec3d auxvec1,auxvec2;
+		Vec<3> auxvec1 = mesh[pi3]-mesh[pi4];
+		Vec<3> auxvec2 = mesh[pi1]-mesh[pi4];
 
-	      auxvec1 = mesh.Point(pi3)-mesh.Point(pi4);
-	      auxvec2 = mesh.Point(pi1)-mesh.Point(pi4);
 		allowswap = allowswap && fabs(1.-(auxvec1*auxvec2)/(auxvec1.Length()*auxvec2.Length())) > 1e-4;
 
 		if(!allowswap)
 		  continue;
 
 		// normal of new
-	      nv1 = Cross (auxvec1, 
+		Vec<3> nv1 = Cross (auxvec1, auxvec2);
-			   auxvec2);
 
 		auxvec1 = mesh.Point(pi4)-mesh.Point(pi3);
 		auxvec2 = mesh.Point(pi2)-mesh.Point(pi3);
@@ -280,16 +297,12 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 		if(!allowswap)
 		  continue;
 
-	      nv2 = Cross (auxvec1, 
+		Vec<3> nv2 = Cross (auxvec1, auxvec2);
-			   auxvec2);
 
 	      
 		// normals of original
-	      Vec3d nv3, nv4;
+		Vec<3> nv3 = Cross (mesh[pi1]-mesh[pi4], mesh[pi2]-mesh[pi4]);
-	      nv3 = Cross (mesh.Point(pi1)-mesh.Point(pi4), 
+		Vec<3> nv4 = Cross (mesh[pi2]-mesh[pi3], mesh[pi1]-mesh[pi3]);
-			   mesh.Point(pi2)-mesh.Point(pi4));
-	      nv4 = Cross (mesh.Point(pi2)-mesh.Point(pi3), 
-			   mesh.Point(pi1)-mesh.Point(pi3));
 	      
 		nv3 *= -1;
 		nv4 *= -1;
@@ -322,7 +335,7 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 		  (nvp4 * nv4 > critval);
 	      
 
-	      horder = Dist (mesh.Point(pi1), mesh.Point(pi2));
+		double horder = Dist (mesh.Point(pi1), mesh.Point(pi2));
 
 		if ( // nv1 * nv2 >= 0 &&
 		    nv1.Length() > 1e-3 * horder * horder &&
@@ -331,8 +344,8 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 		  {
 		    if (!usemetric)
 		      {
-		      e = pdef[pi1] + pdef[pi2] - pdef[pi3] - pdef[pi4];
+			int e = pdef[pi1] + pdef[pi2] - pdef[pi3] - pdef[pi4];
-		      d = 
+			double d = 
 			  Dist2 (mesh.Point(pi1), mesh.Point(pi2)) - 
 			  Dist2 (mesh.Point(pi3), mesh.Point(pi4));
 		      
@@ -340,6 +353,7 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 		      }
 		    else
 		      {
+			double loch = mesh.GetH(mesh[pi1]);
 			should = 
 			  CalcTriangleBadness (mesh.Point(pi4), mesh.Point(pi3), mesh.Point(pi1), 
 					       metricweight, loch) +
@@ -371,13 +385,6 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 		      {
 			// do swapping !
 		      
-		      // cout << "swap " << endl;
-
-		      nswaps ++;
-		      
-		      // testout << "nv1 = " << nv1 << "   nv2 = " << nv2 << endl;
-		     
-
 			done = 1;
 		      
 			mesh[t1].PNum(1) = pi1;
@@ -411,7 +418,7 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
       }
 
     mesh.SetNextTimeStamp();
-}
+  }
 
 
 
@@ -420,8 +427,8 @@ void MeshOptimize2d :: EdgeSwapping (Mesh & mesh, int usemetric)
 
 
  
-void MeshOptimize2d :: CombineImprove (Mesh & mesh)
+  void MeshOptimize2d :: CombineImprove (Mesh & mesh)
-{
+  {
     if (!faceindex)
       {
 	PrintMessage (3, "Combine improve");
@@ -441,18 +448,25 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
     static int timer = NgProfiler::CreateTimer ("Combineimprove 2D");
     NgProfiler::RegionTimer reg (timer);
 
+    static int timerstart = NgProfiler::CreateTimer ("Combineimprove 2D start");
+    NgProfiler::StartTimer  (timerstart);
 
 
-  int i, j, k, l;
+    static int timerstart1 = NgProfiler::CreateTimer ("Combineimprove 2D start1");
-  PointIndex pi;
+    NgProfiler::StartTimer  (timerstart1);
-  SurfaceElementIndex sei;
+
+
+
+    // int i, j, k, l;
+    // PointIndex pi;
+    // SurfaceElementIndex sei;
 
 
     Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
 
 
-  for (i = 0; i < seia.Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
       if (mesh[seia[i]].GetNP() != 3)
 	return;
 
@@ -463,8 +477,8 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
       surfnr = mesh.GetFaceDescriptor (faceindex).SurfNr();
 
 
-  PointIndex pi1, pi2;
+    // PointIndex pi1, pi2;
-  MeshPoint p1, p2, pnew;
+    // MeshPoint p1, p2, pnew;
     double bad1, bad2;
     Vec<3> nv;
 
@@ -474,65 +488,79 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
     TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonnode(np); 
     Array<SurfaceElementIndex> hasonepi, hasbothpi;
 
-  for (i = 0; i < seia.Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
       {
 	Element2d & el = mesh[seia[i]];
-      for (j = 0; j < el.GetNP(); j++)
+	for (int j = 0; j < el.GetNP(); j++)
-	{
 	  elementsonnode.Add (el[j], seia[i]);
       }
-    }
-
 
     Array<bool,PointIndex::BASE> fixed(np);
     fixed = false;
 
-  SegmentIndex si;
+    NgProfiler::StopTimer  (timerstart1);
-  for (si = 0; si < mesh.GetNSeg(); si++)
+
+    /*
+    for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
       {
 	INDEX_2 i2(mesh[si][0], mesh[si][1]);
 	fixed[i2.I1()] = true;
 	fixed[i2.I2()] = true;
       }
+    */
 
-  for(i = 0; i<mesh.LockedPoints().Size(); i++)
+    for (int i = 0; i < seia.Size(); i++)
+      {
+	Element2d & sel = mesh[seia[i]];
+	for (int j = 0; j < sel.GetNP(); j++)
+	  {
+	    PointIndex pi1 = sel.PNumMod(j+2);
+	    PointIndex pi2 = sel.PNumMod(j+3);
+	    if (mesh.IsSegment (pi1, pi2))
+	      {	
+		fixed[pi1] = true;
+		fixed[pi2] = true;
+	      }
+	  }
+      }
+
+
+
+    for(int i = 0; i < mesh.LockedPoints().Size(); i++)
       fixed[mesh.LockedPoints()[i]] = true;
 
 
+
     Array<Vec<3>,PointIndex::BASE> normals(np);
 
-  for (pi = PointIndex::BASE; 
+    for (PointIndex pi = PointIndex::BASE; 
 	 pi < np + PointIndex::BASE; pi++)
       {
 	if (elementsonnode[pi].Size())
 	  {
 	    Element2d & hel = mesh[elementsonnode[pi][0]];
-	  for (k = 0; k < 3; k++)
+	    for (int k = 0; k < 3; k++)
 	      if (hel[k] == pi)
 		{
 		  SelectSurfaceOfPoint (mesh[pi], hel.GeomInfoPi(k+1));
 		  GetNormalVector (surfnr, mesh[pi], hel.GeomInfoPi(k+1), normals[pi]);
 		  break;
 		}
-	  if (k == 3)
+	  }
+      }
+
+    NgProfiler::StopTimer  (timerstart);
+
+    for (int i = 0; i < seia.Size(); i++)
       {
-	      cerr << "Neuer Fehler von Joachim, code 17121" << endl;
+	SurfaceElementIndex sei = seia[i];
-	    }
-	}
-    }
-
-
-  for (i = 0; i < seia.Size(); i++)
-    {
-
-      sei = seia[i];
 	Element2d & elem = mesh[sei];
 	if (elem.IsDeleted()) continue;
 
-      for (j = 0; j < 3; j++)
+	for (int j = 0; j < 3; j++)
 	  {
-	  pi1 = elem[j];
+	    PointIndex pi1 = elem[j];
-	  pi2 = elem[(j+1) % 3];
+	    PointIndex pi2 = elem[(j+1) % 3];
 
 	    if (pi1 < PointIndex::BASE || 
 		pi2 < PointIndex::BASE)
@@ -581,7 +609,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 	    hasonepi.SetSize(0);
 	    hasbothpi.SetSize(0);
 
-	  for (k = 0; k < elementsonnode[pi1].Size(); k++)
+	    for (int k = 0; k < elementsonnode[pi1].Size(); k++)
 	      {
 		const Element2d & el2 = mesh[elementsonnode[pi1][k]];
 
@@ -601,7 +629,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 
 
 	    Element2d & hel = mesh[hasbothpi[0]];
-	  for (k = 0; k < 3; k++)
+	    for (int k = 0; k < 3; k++)
 	      if (hel[k] == pi1)
 		{
 		  SelectSurfaceOfPoint (mesh[pi1],
@@ -609,17 +637,12 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 		  GetNormalVector (surfnr, mesh[pi1], hel.GeomInfoPi(k+1), nv);
 		  break;
 		}
-	  if (k == 3)
-	    {
-	      cerr << "Neuer Fehler von Joachim, code 32434" << endl;
-	    }
-
 
 	    //	  nv = normals.Get(pi1);
 
 
 
-	  for (k = 0; k < elementsonnode[pi2].Size(); k++)
+	    for (int k = 0; k < elementsonnode[pi2].Size(); k++)
 	      {
 		const Element2d & el2 = mesh[elementsonnode[pi2][k]];
 		if (el2.IsDeleted()) continue;
@@ -632,7 +655,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 
 	    bad1 = 0;
 	    int illegal1 = 0, illegal2 = 0;
-	  for (k = 0; k < hasonepi.Size(); k++)
+	    for (int k = 0; k < hasonepi.Size(); k++)
 	      {
 		const Element2d & el = mesh[hasonepi[k]];
 		bad1 += CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]],
@@ -640,7 +663,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 		illegal1 += 1-mesh.LegalTrig(el);
 	      }
 	  
-	  for (k = 0; k < hasbothpi.Size(); k++)
+	    for (int k = 0; k < hasbothpi.Size(); k++)
 	      {
 		const Element2d & el = mesh[hasbothpi[k]];
 		bad1 += CalcTriangleBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]],
@@ -649,15 +672,15 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 	      }
 	    bad1 /= (hasonepi.Size()+hasbothpi.Size());
 
-	  p1 = mesh[pi1];
+	    MeshPoint p1 = mesh[pi1];
-	  p2 = mesh[pi2];
+	    MeshPoint p2 = mesh[pi2];
 
-	  pnew = p1;
+	    MeshPoint pnew = p1;
 	    mesh[pi1] = pnew;
 	    mesh[pi2] = pnew;
 
 	    bad2 = 0;
-	  for (k = 0; k < hasonepi.Size(); k++)
+	    for (int k = 0; k < hasonepi.Size(); k++)
 	      {
 		Element2d & el = mesh[hasonepi[k]];
 		double err = 
@@ -672,7 +695,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 		if (hnv * nv < 0)
 		  bad2 += 1e10;
               
-	      for (l = 0; l < 3; l++)
+		for (int l = 0; l < 3; l++)
 		  if ( (normals[el[l]] * nv) < 0.5)
 		    bad2 += 1e10;
 
@@ -710,7 +733,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 	      
 	      
 		Element2d *el1p(NULL);
-	      l=0;
+		int l = 0;
 		while(mesh[elementsonnode[pi1][l]].IsDeleted() && l<elementsonnode.EntrySize(pi1)) l++;
 		if(l<elementsonnode.EntrySize(pi1))
 		  el1p = &mesh[elementsonnode[pi1][l]];
@@ -725,7 +748,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 		    }
 
 		// (*testout) << "Connect point " << pi2 << " to " << pi1 << "\n";
-	      for (k = 0; k < elementsonnode[pi2].Size(); k++)
+		for (int k = 0; k < elementsonnode[pi2].Size(); k++)
 		  {
 		    Element2d & el = mesh[elementsonnode[pi2][k]];
 		    if(el.IsDeleted()) continue;
@@ -737,7 +760,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 			haspi1 = 1;
 		    if (haspi1) continue;
 
-		  for (l = 0; l < el.GetNP(); l++)
+		    for (int l = 0; l < el.GetNP(); l++)
 		      {
 			if (el[l] == pi2)
 			  {
@@ -759,7 +782,7 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 		  }
 		*/
 
-	      for (k = 0; k < hasbothpi.Size(); k++)
+		for (int k = 0; k < hasbothpi.Size(); k++)
 		  {
 		    mesh[hasbothpi[k]].Delete();
 		    /*
@@ -774,11 +797,11 @@ void MeshOptimize2d :: CombineImprove (Mesh & mesh)
 
     //  mesh.Compress();
     mesh.SetNextTimeStamp();
-}
+  }
 
 
-void MeshOptimize2d :: CheckMeshApproximation (Mesh & mesh)
+  void MeshOptimize2d :: CheckMeshApproximation (Mesh & mesh)
-{
+  {
     // Check angles between elements and normals at corners
     /*
   
@@ -827,5 +850,5 @@ void MeshOptimize2d :: CheckMeshApproximation (Mesh & mesh)
     }
     }
     */
-}
+  }
 }

libsrc/meshing/meshclass.cpp

@@ -896,7 +896,6 @@ namespace netgen
               }
           }
 
-
         if (strcmp (str, "volumeelements") == 0)
           {
             infile >> n;

libsrc/meshing/meshing2.cpp

@@ -196,6 +196,14 @@ namespace netgen
 
   MESHING2_RESULT Meshing2 :: GenerateMesh (Mesh & mesh, double gh, int facenr)
   {
+    static int timer = NgProfiler::CreateTimer ("surface meshing");
+
+    static int timer1 = NgProfiler::CreateTimer ("surface meshing1");
+    static int timer2 = NgProfiler::CreateTimer ("surface meshing2");
+    static int timer3 = NgProfiler::CreateTimer ("surface meshing3");
+    NgProfiler::RegionTimer reg (timer);
+
+
     Array<int> pindex, lindex;
     Array<int> delpoints, dellines;
 
@@ -205,7 +213,6 @@ namespace netgen
     Array<Element2d> locelements;
 
     int z1, z2, oldnp(-1);
-    SurfaceElementIndex sei;
     bool found;
     int rulenr(-1);
     int globind;
@@ -254,7 +261,9 @@ namespace netgen
     double totalarea = Area ();
     double meshedarea = 0;
 
+
     // search tree for surface elements:
+    /*
     for (sei = 0; sei < mesh.GetNSE(); sei++)
       {
 	const Element2d & sel = mesh[sei];
@@ -269,6 +278,31 @@ namespace netgen
 	    box.Add ( mesh[sel[2]] );
 	    surfeltree.Insert (box, sei);
 	  }
+      }
+    */
+    Array<SurfaceElementIndex> seia;
+    mesh.GetSurfaceElementsOfFace (facenr, seia);
+    for (int i = 0; i < seia.Size(); i++)
+      {
+	const Element2d & sel = mesh[seia[i]];
+
+	if (sel.IsDeleted()) continue;
+
+	Box<3> box;
+	box.Set ( mesh[sel[0]] );
+	box.Add ( mesh[sel[1]] );
+	box.Add ( mesh[sel[2]] );
+	surfeltree.Insert (box, i);
+      }
+
+
+
+
+    if (totalarea > 0 || maxarea > 0)
+      for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
+	{
+	  const Element2d & sel = mesh[sei];
+	  if (sel.IsDeleted()) continue;
 	
 	  double trigarea = Cross ( mesh[sel[1]]-mesh[sel[0]],
 				    mesh[sel[2]]-mesh[sel[0]] ).Length() / 2;
@@ -283,6 +317,8 @@ namespace netgen
 	}
 
 
+
+
     const char * savetask = multithread.task;
     multithread.task = "Surface meshing";
 
@@ -293,8 +329,12 @@ namespace netgen
 
     double meshedarea_before = meshedarea;
 
+
     while (!adfront ->Empty() && !multithread.terminate)
       {
+	NgProfiler::RegionTimer reg1 (timer1);
+
+
 	if (multithread.terminate)
 	  throw NgException ("Meshing stopped");
 
@@ -370,6 +410,10 @@ namespace netgen
 
 	adfront ->GetLocals (baselineindex, locpoints, mpgeominfo, loclines, 
 			     pindex, lindex, 2*hinner);
+
+
+	NgProfiler::RegionTimer reg2 (timer2);
+
 	//(*testout) << "h for locals: " << 2*hinner << endl;
 	
 
@@ -707,6 +751,9 @@ namespace netgen
 	      }
 	  }
       
+	NgProfiler::RegionTimer reg3 (timer3);
+
+
 	for (int i = 1; i <= locelements.Size() && found; i++)
 	  {
 	    const Element2d & el = locelements.Get(i);
@@ -1451,14 +1498,15 @@ namespace netgen
 
     PrintMessage (3, "Surface meshing done");
 
+
     adfront->PrintOpenSegments (*testout);
 
     multithread.task = savetask;
 
 
-    //  cout << "surfeltree.depth = " << surfeltree.Tree().Depth() << endl;
     EndMesh ();
 
+
     if (!adfront->Empty())
       return MESHING2_GIVEUP;
     

libsrc/meshing/meshtype.hpp

@@ -232,7 +232,8 @@ class MeshPoint : public Point<3>
 #endif
 
 public:
-  MeshPoint () : layer(1), singular(0.), type(INNERPOINT) 
+  MeshPoint () 
+  // : layer(1), singular(0.), type(INNERPOINT)   // would unnecessarily initialize large arrays !
   { 
 #ifdef PARALLEL
   isghost = 0;

libsrc/meshing/smoothing2.cpp

@@ -667,14 +667,15 @@ namespace netgen
       }
  
 
-
     static int timer = NgProfiler::CreateTimer ("MeshSmoothing 2D");
-    NgProfiler::RegionTimer reg (timer);
+    static int timer1 = NgProfiler::CreateTimer ("MeshSmoothing 2D start");
 
+    NgProfiler::RegionTimer reg (timer);
+    NgProfiler::StartTimer (timer1);
 
     CheckMeshApproximation (mesh);
 
-    SurfaceElementIndex sei;
+    // SurfaceElementIndex sei;
 
     Array<SurfaceElementIndex> seia;
     mesh.GetSurfaceElementsOfFace (faceindex, seia);
@@ -694,14 +695,17 @@ namespace netgen
     PointGeomInfo ngi;
     Point3d origp;
 
+
+
     Vec3d n1, n2;
     Vector x(2), xedge(1);
 
     Array<MeshPoint, PointIndex::BASE> savepoints(mesh.GetNP());
+
     uselocalh = mparam.uselocalh;
 
-    Array<int, PointIndex::BASE> nelementsonpoint(mesh.GetNP());
 
+    Array<int, PointIndex::BASE> nelementsonpoint(mesh.GetNP());
     nelementsonpoint = 0;
     for (int i = 0; i < seia.Size(); i++)
       {
@@ -711,6 +715,7 @@ namespace netgen
       }
 
     TABLE<SurfaceElementIndex,PointIndex::BASE> elementsonpoint(nelementsonpoint);
+
     for (int i = 0; i < seia.Size(); i++)
       {
 	const Element2d & el = mesh[seia[i]];
@@ -718,6 +723,7 @@ namespace netgen
 	  elementsonpoint.Add (el[j], seia[i]);
       }
 
+
     loch = mparam.maxh;
     locmetricweight = metricweight;
     meshthis = this;
@@ -813,6 +819,10 @@ namespace netgen
       }
     int cnt = 0;
 
+
+    NgProfiler::StopTimer (timer1);
+
+
     for (PointIndex pi = PointIndex::BASE; pi < mesh.GetNP()+PointIndex::BASE; pi++)
       if (mesh[pi].Type() == SURFACEPOINT)
 	{
@@ -829,6 +839,7 @@ namespace netgen
 	  if (elementsonpoint[pi].Size() == 0)
 	    continue;
 
+	
 	  sp1 = mesh[pi];
 
 	  Element2d & hel = mesh[elementsonpoint[pi][0]];
@@ -850,7 +861,7 @@ namespace netgen
 	
 	  for (int j = 0; j < elementsonpoint[pi].Size(); j++)
 	    {
-	      sei = elementsonpoint[pi][j];
+	      SurfaceElementIndex sei = elementsonpoint[pi][j];
 	      const Element2d & bel = mesh[sei];
 	      surfi = mesh.GetFaceDescriptor(bel.GetIndex()).SurfNr();