findElInfoAtPoint for non-convex geometries

AMDiS/src/Mesh.cc

@@ -705,20 +705,31 @@ namespace AMDiS {
     // second time, what can happen with periodic boundary conditions, the point is
     // not within the mesh!
     std::set<int> macrosVisited;
-    macrosVisited.insert(mel->getIndex());
+    std::stack<MacroElement*> active;
+    
+//     macrosVisited.insert(mel->getIndex());
 
     int k;
     while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
-      if (mel->getNeighbour(k)) {
+      macrosVisited.insert(mel->getIndex());
+      if (mel->getNeighbour(k) && !macrosVisited.count(mel->getNeighbour(k)->getIndex()) {
 	mel = mel->getNeighbour(k);
-	if (macrosVisited.count(mel->getIndex()))
-	  return false;
+// 	if (macrosVisited.count(mel->getIndex())) // Nur für periodische RB interessant, muss noch implementiert werden
+// 	  return false;
 
-	macrosVisited.insert(mel->getIndex());
 	mel_info->fillMacroInfo(mel);
 	continue;
+      } else {
+	for (int i = 0; i < dim + 1; ++i) {
+	  if (i !=  k && mel->getNeighbour(i) && !macrosVisited.count(mel->getNeighbour(i)->getIndex())
+	    active.push(mel->getNeighbour(i));    
+	}
+	if (active.empty())
+	  return false;
+	mel = active.pop();
+
+	mel_info->fillMacroInfo(mel);
       }
-      break;
     }
 
     /* now, descend in tree to find leaf element at point */