From da42e35805d6a607c22a805368a425a3e839128f Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Fri, 14 Oct 2011 16:22:14 +0000
Subject: [PATCH] findElInfoAtPoint for non-convex geometries

---
 AMDiS/src/Mesh.cc | 23 +++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 63553c60..34cb2a6a 100644
--- a/AMDiS/src/Mesh.cc
+++ b/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 */
-- 
GitLab