2d meshing improvements

libsrc/general/array.hpp

@@ -393,6 +393,16 @@ namespace netgen
       Array<T>::operator= (val);
       return *this;
     }
+
+    /// array copy
+    ArrayMem & operator= (const FlatArray<T> & a2)
+    {
+      SetSize (a2.Size());
+      for (int i = 0; i < this->size; i++)
+	(*this)[i] = a2[i];
+      return *this;
+    }
+
   };
 
 

libsrc/meshing/netrule2.cpp

@@ -95,18 +95,12 @@ void netrule :: SetFreeZoneTransformation (const Vector & devp, int tolclass)
 	}
       else
 	{
-	  vn /= sqrt (len2);    // should not be necessary
+	  vn /= sqrt (len2);    // scaling necessary ?
 
 	  freesetinequ(i,0) = vn.X(); 
 	  freesetinequ(i,1) = vn.Y(); 
 	  freesetinequ(i,2) = -(p1.X() * vn.X() + p1.Y() * vn.Y());
 	}
-
-      /*
-      freesetinequ(i,0) = vn.X(); 
-      freesetinequ(i,1) = vn.Y(); 
-      freesetinequ(i,2) = -(p1.X() * vn.X() + p1.Y() * vn.Y());
-      */
     }
 }
 

libsrc/meshing/ruler2.cpp

@@ -48,160 +48,193 @@ static double CalcElementBadness (const Array<Point2d> & points,
   int Meshing2 ::ApplyRules (Array<Point2d> & lpoints, 
 			     Array<int> & legalpoints,
 			     int maxlegalpoint,
-			   Array<INDEX_2> & llines,
+			     Array<INDEX_2> & llines1,
 			     int maxlegalline,
 			     Array<Element2d> & elements,
 			     Array<INDEX> & dellines, int tolerance)
   {
-  int i, j, ri, nlok, npok, incnpok, refpi, locli = 0;
+    static int timer = NgProfiler::CreateTimer ("meshing2::ApplyRules");
+    NgProfiler::RegionTimer reg (timer);
+
+
 
     double maxerr = 0.5 + 0.3 * tolerance;
     double minelerr = 2 + 0.5 * tolerance * tolerance;
 
+    int noldlp = lpoints.Size();
+    int noldll = llines1.Size();
 
-  bool ok;
-  int found;   // rule number
-  Vector oldu, newu;
-  Point2d np;
-  Vec2d linevec;
-  int oldnp;
-  INDEX_2 loclin;
-  double hf, elerr;
-  int noldlp, noldll;
-  int loctestmode;
 
-  static Array<int> pused, pmap, pfixed;
+    ArrayMem<int,100> pused(maxlegalpoint), lused(maxlegalline);
-  static Array<int, 1> lmap, lused;
+    ArrayMem<int,100> pnearness(noldlp), lnearness(llines1.Size());
-  static Array<int> pnearness, lnearness;
 
-  static Array<Point2d> tempnewpoints;
+    ArrayMem<int, 20> pmap, pfixed, lmap;
-  static Array<INDEX_2> tempnewlines;
-  static Array<int> tempdellines;
-  static Array<Element2d> tempelements;
   
+    ArrayMem<Point2d,100> tempnewpoints;
+    ArrayMem<INDEX_2,100> tempnewlines;
+    ArrayMem<int,100> tempdellines;
+    ArrayMem<Element2d,100> tempelements;
 
 
     elements.SetSize (0);
     dellines.SetSize (0);
 
-  noldlp = lpoints.Size();
-  noldll = llines.Size();
-
-  pused.SetSize (maxlegalpoint);
-  lused.SetSize (maxlegalline);
-  pnearness.SetSize (noldlp);
-  lnearness.SetSize (llines.Size());
-
-
     testmode = debugparam.debugoutput;
-  loctestmode = testmode;
+
+#ifdef LOCDEBUG
+    int loctestmode = testmode;
 
     if (loctestmode)
       {
 	(*testout) << endl << endl << "Check new environment" << endl;
 	(*testout) << "tolerance = " << tolerance << endl;
-      for (i = 1; i <= lpoints.Size(); i++)
+	for (int i = 1; i <= lpoints.Size(); i++)
 	  (*testout) << "P" << i << " = " << lpoints.Get(i) << endl;
 	(*testout) << endl;
-      for (i = 1; i <= llines.Size(); i++)
+	for (int i = 1; i <= llines1.Size(); i++)
-	(*testout) << "(" << llines.Get(i).I1() << "-" << llines.Get(i).I2() << ")" << endl;
+	  (*testout) << "(" << llines1.Get(i).I1() << "-" << llines1.Get(i).I2() << ")" << endl;
       }
+#endif
 
     // check every rule
 
-  found = 0;
+    int found = 0;   // rule number
+
+    pnearness = 1000;
+  
+    for (int j = 0; j < 2; j++)
+      pnearness.Set(llines1[0][j], 0);
 
 
-  for (i = 1; i <= noldlp; i++)
-    pnearness.Set(i, 1000);
 
-  for (j = 1; j <= 2; j++)
+    enum { MAX_NEARNESS = 3 };
-    pnearness.Set(llines.Get(1).I(j), 0);
 
-
+    for (int cnt = 0; cnt < MAX_NEARNESS; cnt++)
-  do
       {
-      ok = 1;
+	bool ok = true;
-      for (i = 1; i <= maxlegalline; i++)
+	for (int i = 0; i < maxlegalline; i++)
 	  {
-	  const INDEX_2 & hline = llines.Get(i);
+	    const INDEX_2 & hline = llines1[i];
 
-	  /*
+	    int minn = min2 (pnearness.Get(hline[0]),  pnearness.Get(hline[1]));
-	  int minn = INT_MAX-1;
-	  for (j = 1; j <= 2; j++)
-	    {
-	      int hi = pnearness.Get(hline.I(j));
-	      if (hi < minn) minn = hi;
-	    }
-	  */
-	  int minn = pnearness.Get(hline.I1());
-	  int minn2 = pnearness.Get(hline.I2());
-	  if (minn2 < minn)
-	    minn = minn2;
 
-	  /*
+	    for (int j = 0; j < 2; j++)
-	  for (j = 1; j <= 2; j++)
+	      if (pnearness.Get(hline[j]) > minn+1)
 		{
-	      int hpi = hline.I(j);
+		  ok = false;
-	      if (pnearness.Get(hpi) > minn+1)
+		  pnearness.Set(hline[j], minn+1);
-		{
-		  ok = 0;
-		  pnearness.Set(hpi, minn+1);
 		}
 	  }
-	  */
+	if (!ok) break;
-	  int hpi = hline.I1();
-	  if (pnearness.Get(hpi) > minn+1)
-	    {
-	      ok = 0;
-	      pnearness.Set(hpi, minn+1);
-	    }
-	  hpi = hline.I2();
-	  if (pnearness.Get(hpi) > minn+1)
-	    {
-	      ok = 0;
-	      pnearness.Set(hpi, minn+1);
-	    }
-	}
-    }
-  while (!ok);
-  
-  for (i = 1; i <= maxlegalline /* lnearness.Size() */; i++)
-    {
-      lnearness.Set(i, 0);
-      for (j = 1; j <= 2; j++)
-	lnearness.Elem(i) += pnearness.Get(llines.Get(i).I(j));
       }
 
 
-  for (ri = 1; ri <= rules.Size(); ri++)
+    for (int i = 0; i < maxlegalline; i++)
+      lnearness[i] = pnearness.Get(llines1[i][0]) + pnearness.Get(llines1[i][1]);
+
+
+    // resort lines after lnearness
+    Array<INDEX_2> llines(llines1.Size());
+    Array<int> sortlines(llines1.Size());
+    int lnearness_class[MAX_NEARNESS];
+
+    for (int j = 0; j < MAX_NEARNESS; j++)
+      lnearness_class[j] = 0;
+    for (int i = 0; i < maxlegalline; i++)
+      if (lnearness[i] < MAX_NEARNESS)
+	lnearness_class[lnearness[i]]++;
+    
+    int cumm = 0;
+    for (int j = 0; j < MAX_NEARNESS; j++)
       {
+	int hcnt = lnearness_class[j];
+	lnearness_class[j] = cumm;
+	cumm += hcnt;
+      }
+
+    for (int i = 0; i < maxlegalline; i++)
+      if (lnearness[i] < MAX_NEARNESS)
+	{
+	  llines[lnearness_class[lnearness[i]]] = llines1[i];
+	  sortlines[lnearness_class[lnearness[i]]] = i+1;
+	  lnearness_class[lnearness[i]]++;
+	}
+      else
+	{
+	  llines[cumm] = llines1[i];
+	  sortlines[cumm] = i+1;
+	  cumm++;
+	}
+
+    for (int i = maxlegalline; i < llines1.Size(); i++)
+      {
+	llines[cumm] = llines1[i];
+	sortlines[cumm] = i+1;
+	cumm++;
+      }
+
+    for (int i = 0; i < maxlegalline; i++)
+      lnearness[i] = pnearness.Get(llines[i][0]) + pnearness.Get(llines[i][1]);
+
+
+
+
+    static bool firsttime = true;
+    static int timers[100];
+    static int timers2[100];
+    static int timers3[100];
+    if (firsttime)
+      {
+	for (int ri = 0; ri < rules.Size(); ri++)
+	  timers[ri] = NgProfiler::CreateTimer (string("netrule ")+rules[ri]->Name());
+	for (int ri = 0; ri < rules.Size(); ri++)
+	  timers2[ri] = NgProfiler::CreateTimer (string("netrule,mapped ")+rules[ri]->Name());
+	for (int ri = 0; ri < rules.Size(); ri++)
+	  timers3[ri] = NgProfiler::CreateTimer (string("netrule,lines mapped ")+rules[ri]->Name());
+	firsttime = false;
+      }
+
+    lused = 0;
+    pused = 0;
+
+
+    static int timer1 = NgProfiler::CreateTimer ("meshing2::ApplyRules 1");
+    NgProfiler::RegionTimer reg1 (timer1);
+
+
+    for (int ri = 1; ri <= rules.Size(); ri++)
+      {
+	NgProfiler::RegionTimer reg(timers[ri-1]);
 	netrule * rule = rules.Get(ri);
 
+#ifdef LOCDEBUG
 	if (loctestmode)
 	  (*testout) << "Rule " << rule->Name() << endl;
+#endif
 
 	if (rule->GetQuality() > tolerance) continue;
 
 	pmap.SetSize (rule->GetNP());
 	lmap.SetSize (rule->GetNL());
       
-      lused = 0;
-      pused = 0;
 	pmap = 0;
 	lmap = 0;
 
-      lused[1] = 1;   // .Set (1, 1);
+	lused[0] = 1; 
-      lmap[1] = 1;    // .Set (1, 1);
+	lmap[0] = 1;  
 
-      for (j = 0; j < 2; j++)
+	for (int j = 0; j < 2; j++)
 	  {
 	    pmap.Elem(rule->GetLine(1)[j]) = llines[0][j];
 	    pused.Elem(llines[0][j])++;
 	  }
 
 
-      nlok = 2;
+
+	int nlok = 2;
+
+
+	bool ok = false;
 
 	while (nlok >= 2)
 	  {
@@ -210,38 +243,48 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 
 	      {
 		ok = 0;
-	      while (!ok && lmap.Get(nlok) < maxlegalline  /* llines.Size() */)
+		
+		int maxline = (rule->GetLNearness(nlok) < MAX_NEARNESS) ? lnearness_class[rule->GetLNearness(nlok)] : maxlegalline;
+		// int maxline = maxlegalline;
+
+		while (!ok && lmap.Get(nlok) < maxline)
 		  {
 		    lmap.Elem(nlok)++;
-		  locli = lmap.Get(nlok);
+		    int locli = lmap.Get(nlok);
+
+		    if (lnearness.Get(locli) > rule->GetLNearness (nlok) ) continue;
+		    if (lused.Get(locli)) continue;
+
 
-		  if (!lused.Get(locli) && 
-		      lnearness.Get(locli) <= rule->GetLNearness (nlok) )
-		    {
 		    ok = 1;
 
-		      loclin = llines.Get(locli);
+		    INDEX_2 loclin = llines.Get(locli);
-                      linevec = lpoints.Get(loclin.I2()) - lpoints.Get(loclin.I1());
+		    Vec2d linevec = lpoints.Get(loclin.I2()) - lpoints.Get(loclin.I1());
 
 		    if (rule->CalcLineError (nlok, linevec) > maxerr)
 		      {
 			ok = 0;
+#ifdef LOCDEBUG
 			if(loctestmode)
 			  (*testout) << "not ok pos1" << endl;
+#endif
+			continue;
 		      }
 
-		      for (j = 0; j < 2 && ok; j++)
+		    for (int j = 0; j < 2; j++)
 		      {
-			  // refpi = rule->GetPointNr (nlok, j);
+			int refpi = rule->GetLine(nlok)[j];
-                          refpi = rule->GetLine(nlok)[j];
 
 			if (pmap.Get(refpi) != 0)
 			  {
 			    if (pmap.Get(refpi) != loclin[j])
 			      {
 				ok = 0;
+#ifdef LOCDEBUG
 				if(loctestmode)
 				  (*testout) << "not ok pos2" << endl;
+#endif
+				break;
 			      }
 			  }
 			else
@@ -251,6 +294,7 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 				|| pused.Get(loclin[j]))
 			      {
 				ok = 0;
+#ifdef LOCDEBUG
 				if(loctestmode)
 				  {
 				    (*testout) << "nok pos3" << endl;
@@ -261,7 +305,8 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 				    //if(pused.Get(loclin[j]))
 				    //(*testout) << "r3" << endl;
 				  }
-				}
+#endif
+				break;
 			      }
 			  }
 		      }
@@ -269,8 +314,11 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 
 		if (ok)
 		  {
+		    int locli = lmap.Get(nlok);
+		    INDEX_2 loclin = llines.Get(locli);
+
 		    lused.Elem (locli) = 1;
-		  for (j = 0; j < 2; j++)
+		    for (int j = 0; j < 2; j++)
 		      {
 			pmap.Set(rule->GetLine (nlok)[j], loclin[j]);
 			pused.Elem(loclin[j])++;
@@ -284,7 +332,7 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 		    nlok--;
 
 		    lused.Elem (lmap.Get(nlok)) = 0;
-		  for (j = 0; j < 2; j++)
+		    for (int j = 0; j < 2; j++)
 		      {
 			pused.Elem(llines.Get(lmap.Get(nlok))[j]) --;
 			if (! pused.Get (llines.Get (lmap.Get (nlok))[j]))
@@ -296,16 +344,17 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 	    else
 
 	      {
+		NgProfiler::RegionTimer reg(timers3[ri-1]);
 
 		// all lines are mapped !!
 
 		// map also all points:
 
-	      npok = 1;
+		int npok = 1;
-	      incnpok = 1;
+		int incnpok = 1;
 
 		pfixed.SetSize (pmap.Size());
-	      for (i = 0; i < pmap.Size(); i++)
+		for (int i = 0; i < pmap.Size(); i++)
 		  pfixed[i] = (pmap[i] >= 1);
  
 		while (npok >= 1)
@@ -370,32 +419,42 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 		    else
 
 		      {
+			NgProfiler::RegionTimer reg(timers2[ri-1]);
+
+			npok = rule->GetNOldP();
+			incnpok = 0;
+
 			if (ok)
 			  foundmap.Elem(ri)++; 
 
+#ifdef LOCDEBUG
 			if (loctestmode)
 			  (*testout) << "lines and points mapped" << endl;
-
+#endif
 
 			ok = 1;
 
 			// check orientations
 
-		      for (i = 1; i <= rule->GetNOrientations() && ok; i++)
+			for (int i = 1; i <= rule->GetNOrientations(); i++)
 			  {
 			    if (CW (lpoints.Get(pmap.Get(rule->GetOrientation(i).i1)),
 				    lpoints.Get(pmap.Get(rule->GetOrientation(i).i2)),
 				    lpoints.Get(pmap.Get(rule->GetOrientation(i).i3))) )
 			      {
 				ok = 0;
+#ifdef LOCDEBUG
 				if (loctestmode)
 				  (*testout) << "Orientation " << i << " not ok" << endl;
+#endif
+				break;
 			      }
 			  }
 
-		      if (ok)
+
-			{
+			if (!ok) continue;
-			  oldu.SetSize (2 * rule->GetNOldP());
+
+			Vector oldu (2 * rule->GetNOldP());
 		      
 			for (int i = 1; i <= rule->GetNOldP(); i++)
 			  {
@@ -405,15 +464,16 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 			  }
 		      
 			rule -> SetFreeZoneTransformation (oldu, tolerance);
-			}
 
 		      
-		      if (ok && !rule->ConvexFreeZone())
+			if (!ok) continue;
+			if (!rule->ConvexFreeZone())
 			  {
 			    ok = 0;
+#ifdef LOCDEBUG
 			    if (loctestmode) 
 			      (*testout) << "freezone not convex" << endl;
-
+#endif
 			    /*
 			      static int cnt = 0;
 			      cnt++;
@@ -429,50 +489,67 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 			  }
 
 			// check freezone:
-
+			if (!ok) continue;
-		      for (i = 1; i <= maxlegalpoint && ok; i++)
+			for (int i = 1; i <= maxlegalpoint && ok; i++)
 			  {
 			    if ( !pused.Get(i) &&
 				 rule->IsInFreeZone (lpoints.Get(i)) )
 			      {
 				ok = 0;
+#ifdef LOCDEBUG
 				if (loctestmode)
 				  (*testout) << "Point " << i << " in freezone" << endl;
+#endif
+				break;
 			      }
 			  }
 
-
+			if (!ok) continue;
-		      for (i = maxlegalpoint+1; i <= lpoints.Size() && ok; i++)
+			for (int i = maxlegalpoint+1; i <= lpoints.Size(); i++)
 			  {
 			    if ( rule->IsInFreeZone (lpoints.Get(i)) )
 			      {
 				ok = 0;
+#ifdef LOCDEBUG
 				if (loctestmode)
 				  (*testout) << "Point " << i << " in freezone" << endl;
+#endif
+				break;
 			      }
 			  }
 
-		      for (i = 1; i <= maxlegalline && ok; i++)
+
+			if (!ok) continue;
+			for (int i = 1; i <= maxlegalline; i++)
 			  {
 			    if (!lused.Get(i) && 
 				rule->IsLineInFreeZone (lpoints.Get(llines.Get(i).I1()),
 							lpoints.Get(llines.Get(i).I2())))
 			      {
 				ok = 0;
+#ifdef LOCDEBUG
 				if (loctestmode)
 				  (*testout) << "line " << llines.Get(i).I1() << "-"
 					     << llines.Get(i).I2() << " in freezone" << endl;
+#endif
+				break;
 			      }
 			  }
-		      for (i = maxlegalline+1; i <= llines.Size() && ok; i++)
+
+			if (!ok) continue;
+
+			for (int i = maxlegalline+1; i <= llines.Size(); i++)
 			  {
 			    if (rule->IsLineInFreeZone (lpoints.Get(llines.Get(i).I1()),
 							lpoints.Get(llines.Get(i).I2())))
 			      {
 				ok = 0;
+#ifdef LOCDEBUG
 				if (loctestmode)
 				  (*testout) << "line " << llines.Get(i).I1() << "-"
 					     << llines.Get(i).I2() << " in freezone" << endl;
+#endif
+				break;
 			      }
 			  }
 
@@ -494,34 +571,33 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 			*/
 
 
-		      if (ok)
+			if (!ok) continue;
-			{
+
+#ifdef LOCDEBUG
 			if (loctestmode)
 			  (*testout) << "rule ok" << endl;
-
+#endif
-			  //			  newu = rule->GetOldUToNewU() * oldu;
-			  if (rule->GetNOldP() < rule->GetNP())
-			    {
-			      newu.SetSize (rule->GetOldUToNewU().Height());
-			      rule->GetOldUToNewU().Mult (oldu, newu);
-			    }
 
 			// Setze neue Punkte:
-
+			if (rule->GetNOldP() < rule->GetNP())
-			  oldnp = rule->GetNOldP();
-
-			  for (i = oldnp + 1; i <= rule->GetNP(); i++)
 			  {
-			      np = rule->GetPoint(i);
+			    Vector newu(rule->GetOldUToNewU().Height());
+			    rule->GetOldUToNewU().Mult (oldu, newu);
+			    
+			    int oldnp = rule->GetNOldP();
+			    for (int i = oldnp + 1; i <= rule->GetNP(); i++)
+			      {
+				Point2d np = rule->GetPoint(i);
 				np.X() += newu (2 * (i-oldnp) - 2);
 				np.Y() += newu (2 * (i-oldnp) - 1);
 				
 				pmap.Elem(i) = lpoints.Append (np);
 			      }
+			  }
 
 			// Setze neue Linien:
 
-			  for (i = rule->GetNOldL() + 1; i <= rule->GetNL(); i++)
+			for (int i = rule->GetNOldL() + 1; i <= rule->GetNL(); i++)
 			  {
 			    llines.Append (INDEX_2 (pmap.Get(rule->GetLine (i)[0]),
 						    pmap.Get(rule->GetLine (i)[1])));
@@ -529,60 +605,65 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 
 
 			// delete old lines:
-                          for (i = 1; i <= rule->GetNDelL(); i++)
+			for (int i = 1; i <= rule->GetNDelL(); i++)
-                            dellines.Append (lmap.Get(rule->GetDelLine(i)));
+			  dellines.Append (sortlines.Elem (lmap.Get(rule->GetDelLine(i))));
+			// dellines.Append (lmap.Get(rule->GetDelLine(i))));
 
-                          // dellines.Append (lmap[rule->GetDelLines()]);
+			// dellines.Append (lmap.Elem(rule->GetDelLines()));
 			// lmap[rule->GetDelLines()];
 
 
 			// insert new elements:
 
-			  for (i = 1; i <= rule->GetNE(); i++)
+			for (int i = 1; i <= rule->GetNE(); i++)
 			  {
 			    elements.Append (rule->GetElement(i));
-			      for (j = 1; j <= elements.Get(i).GetNP(); j++)
+			    for (int j = 1; j <= elements.Get(i).GetNP(); j++)
 			      elements.Elem(i).PNum(j) = pmap.Get(elements.Get(i).PNum(j));
 			  }
 
 
-			  elerr = 0;
+			double elerr = 0;
-			  for (i = 1; i <= elements.Size(); i++)
+			for (int i = 1; i <= elements.Size(); i++)
 			  {
+			    double hf;
 			    if (!mparam.quad)
 			      hf = CalcElementBadness (lpoints, elements.Get(i));
 			    else
 			      hf = elements.Get(i).CalcJacobianBadness (lpoints) * 5;
+#ifdef LOCDEBUG
 			    if (loctestmode)
 			      (*testout) << "r " << rule->Name() << "bad = " << hf << endl;
+#endif
 			    if (hf > elerr) elerr = hf;
 			  }
 
+#ifdef LOCDEBUG
 			if (loctestmode)
 			  (*testout) << "error = " << elerr;
-
+#endif
 
 			canuse.Elem(ri) ++;
 
 			if (elerr < 0.99*minelerr)
 			  {
-
+#ifdef LOCDEBUG
 			    if (loctestmode)
 			      {
 				(*testout) << "rule = " << rule->Name() << endl;
 				(*testout) << "class = " << tolerance << endl;
 				(*testout) << "lpoints: " << endl;
-				  for (i = 1; i <= lpoints.Size(); i++)
+				for (int i = 1; i <= lpoints.Size(); i++)
 				  (*testout) << lpoints.Get(i) << endl;
 				(*testout) << "llines: " << endl;
-				  for (i = 1; i <= llines.Size(); i++)
+				for (int i = 1; i <= llines.Size(); i++)
 				  (*testout) << llines.Get(i).I1() << " " << llines.Get(i).I2() << endl;
 
 				(*testout) << "Freezone: ";
-				  for (i = 1; i <= rule -> GetTransFreeZone().Size(); i++)
+				for (int i = 1; i <= rule -> GetTransFreeZone().Size(); i++)
 				  (*testout) << rule->GetTransFreeZone().Get(i) << endl;
 			      }
-
+#endif
 
 			    minelerr = elerr;
 			    found = ri;
@@ -599,19 +680,15 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
 			elements.SetSize (0);
 			ok = 0;
 		      }
-
-		      npok = rule->GetNOldP();
-		      incnpok = 0;
-		    }
 		  }
 
 		nlok = rule->GetNOldL();
 
 		lused.Set (lmap.Get(nlok), 0);
 
-	      for (j = 1; j <= 2; j++)
+		for (int j = 1; j <= 2; j++)
 		  {
-		  refpi = rule->GetPointNr (nlok, j);
+		    int refpi = rule->GetPointNr (nlok, j);
 		    pused.Elem(pmap.Get(refpi))--;
 
 		    if (pused.Get(pmap.Get(refpi)) == 0)
@@ -625,7 +702,7 @@ int Meshing2 ::ApplyRules (Array<Point2d> & lpoints,
     if (found)
       {
 	lpoints.Append (tempnewpoints);
-      llines.Append (tempnewlines);
+	llines1.Append (tempnewlines);
 	dellines.Append (tempdellines);
 	elements.Append (tempelements);
       }

libsrc/meshing/topology.cpp

@@ -211,7 +211,7 @@ namespace netgen
 		edgeflag[i2] = i;
 		edgenr[i2] = ednr;
 	      }
-	    for (int j = 0; j << vert2vertcoarse[i].Size(); j++)
+	    for (int j = 0; j < vert2vertcoarse[i].Size(); j++)      // fix by Markus
 	      {
 		int v2 = vert2vertcoarse[i][j];
 		if (edgeflag[v2] < i)

libsrc/occ/occgenmesh.cpp

@@ -20,6 +20,8 @@ namespace netgen
 
    bool merge_solids = 1;
 
+
+  // can you please explain what you intend to compute here (JS) !!!
    double Line :: Dist (Line l)
    {
       Vec<3> n = p1-p0;

libsrc/visualization/soldata.hpp

@@ -26,45 +26,71 @@ namespace netgen
     virtual ~SolutionData ()
     { ; }
 
-    int GetComponents() { return components; }
+    int GetComponents() 
-    bool IsComplex() { return iscomplex; }
+    { 
+      return components; 
+    }
+
+    bool IsComplex() 
+    {
+      return iscomplex; 
+    }
 
     virtual bool GetValue (int /* elnr */, 
                            double /* lam1 */, double /* lam2 */, double /* lam3 */,
                            double * /* values */) 
-    { return false; }
+    { 
+      return false; 
+    }
 
     virtual bool GetValue (int selnr,
                            const double xref[], const double x[], const double dxdxref[],
                            double * values)
-    { return GetValue (selnr, xref[0], xref[1], xref[2], values); }
+    {
+      return GetValue (selnr, xref[0], xref[1], xref[2], values); 
+    }
 
     virtual bool GetMultiValue (int elnr, int npts,
 				const double * xref, int sxref,
 				const double * x, int sx,
 				const double * dxdxref, int sdxdxref,
-				double * values, int svalues);
+				double * values, int svalues)
+    {
+      bool res = false;
+      for (int i = 0; i < npts; i++)
+	res = GetValue (elnr, &xref[i*sxref], &x[i*sx], &dxdxref[i*sdxdxref], &values[i*svalues]);
+      return res;
+    }
 
 
 
     virtual bool GetSurfValue (int /* selnr */,
                                double /* lam1 */, double /* lam2 */, 
                                double * /* values */)
-    { return false; }
+    { 
+      return false; 
+    }
 
 
     virtual bool GetSurfValue (int selnr,
                                const double xref[], const double x[], const double dxdxref[],
                                double * values)
-    { return GetSurfValue (selnr, xref[0], xref[1], values); }
+    { 
+      return GetSurfValue (selnr, xref[0], xref[1], values); 
+    }
 
 
     virtual bool GetMultiSurfValue (int selnr, int npts,
                                     const double * xref, int sxref,
                                     const double * x, int sx,
                                     const double * dxdxref, int sdxdxref,
-                                    double * values, int svalues);
+                                    double * values, int svalues)
-
+    {
+      bool res = false;
+      for (int i = 0; i < npts; i++)
+	res = GetSurfValue (selnr, &xref[i*sxref], &x[i*sx], &dxdxref[i*sdxdxref], &values[i*svalues]);
+      return res;
+    }
 
     void SetMultiDimComponent (int mc)
     { multidimcomponent = mc; }

libsrc/visualization/vssolution.cpp

@@ -29,6 +29,7 @@ namespace netgen
     delete solclass;
   }
 
+  /*
   bool SolutionData :: GetMultiValue (int elnr, int npts,
 				     const double * xref, int sxref,
 				     const double * x, int sx,
@@ -53,6 +54,7 @@ namespace netgen
       res = GetSurfValue (selnr, &xref[i*sxref], &x[i*sx], &dxdxref[i*sdxdxref], &values[i*svalues]);
     return res;
   }
+  */
   
   VisualSceneSolution :: VisualSceneSolution ()
     : VisualScene()