Changes in molecuilder/src/tesselation.cpp [246a3c:591f15]

File:

• molecuilder/src/tesselation.cpp (modified) (74 diffs)

molecuilder/src/tesselation.cpp

@@ -35 +35 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
   LinesCount(0),
   node(Walker),
@@ -61 +61 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
 {
         Info FunctionInfo(__func__);
@@ -105 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
+BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
 {
         Info FunctionInfo(__func__);
@@ -115 +115 @@
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
+  // set skipped to false
+  skipped = false;
+  // clear triangles list
+  Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
+};
+
+/** Constructor of BoundaryLineSet with two endpoints.
+ * Adds line automatically to each endpoints' LineMap
+ * \param *Point1 first boundary point
+ * \param *Point2 second boundary point
+ * \param number number of the list
+ */
+BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number)
+{
+  Info FunctionInfo(__func__);
+  // set number
+  Nr = number;
+  // set endpoints in ascending order
+  SetEndpointsOrdered(endpoints, Point1, Point2);
+  // add this line to the hash maps of both endpoints
+  Point1->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
+  Point2->AddLine(this); //
   // set skipped to false
   skipped = false;
@@ -171 +193 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
 {
         Info FunctionInfo(__func__);
@@ -182 +204 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
+bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
 {
         Info FunctionInfo(__func__);
@@ -197 +219 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion()
+bool BoundaryLineSet::CheckConvexityCriterion() const
 {
         Info FunctionInfo(__func__);
@@ -221 +243 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
     //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     NormalCheck.AddVector(&runner->second->NormalVector);
@@ -264 +286 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
         Info FunctionInfo(__func__);
@@ -277 +299 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
 {
         Info FunctionInfo(__func__);
@@ -317 +339 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
   Nr(number)
 {
@@ -376 +398 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
 {
         Info FunctionInfo(__func__);
@@ -388 +410 @@
 };
 
-/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
  * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
@@ -400 +422 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
-{
-        Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
+{
+  Info FunctionInfo(__func__);
   Vector CrossPoint;
   Vector helper;
@@ -411 +433 @@
   }
 
+  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
+  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
+
+  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl;
+    return true;
+  }
   // Calculate cross point between one baseline and the line from the third endpoint to intersection
   int i=0;
@@ -417 +453 @@
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
+      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+      const double s = CrossPoint.ScalarProduct(&helper)/helper.NormSquared();
+      Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl;
+      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
+        Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl;
+        i=4;
+        break;
+      }
+      i++;
     } else 
-      i++;
-    if (i>2)
       break;
-  } while (CrossPoint.NormSquared() < MYEPSILON);
+  } while (i<3);
   if (i==3) {
-    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
-  }
-  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-
-  // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
     return true;
-  } else { // outside!
-    Intersection->Zero();
+  } else {
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
     return false;
   }
+};
+
+/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
+ * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+ * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
+ * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
+ * the first two basepoints) or not.
+ * \param *x point
+ * \param *ClosestPoint desired closest point inside triangle to \a *x, is absolute vector
+ * \return Distance squared between \a *x and closest point inside triangle
+ */
+double BoundaryTriangleSet::GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const
+{
+  Info FunctionInfo(__func__);
+  Vector Direction;
+
+  // 1. get intersection with plane
+  Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
+  GetCenter(&Direction);
+  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
+    ClosestPoint->CopyVector(x);
+  }
+
+  // 2. Calculate in plane part of line (x, intersection)
+  Vector InPlane;
+  InPlane.CopyVector(x);
+  InPlane.SubtractVector(ClosestPoint);  // points from plane intersection to straight-down point
+  InPlane.ProjectOntoPlane(&NormalVector);
+  InPlane.AddVector(ClosestPoint);
+
+  Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl;
+  Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl;
+
+  // Calculate cross point between one baseline and the desired point such that distance is shortest
+  double ShortestDistance = -1.;
+  bool InsideFlag = false;
+  Vector CrossDirection[3];
+  Vector CrossPoint[3];
+  Vector helper;
+  for (int i=0;i<3;i++) {
+    // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
+    Direction.CopyVector(endpoints[(i+1)%3]->node->node);
+    Direction.SubtractVector(endpoints[i%3]->node->node);
+    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
+    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
+    CrossDirection[i].CopyVector(&CrossPoint[i]);
+    CrossDirection[i].SubtractVector(&InPlane);
+    CrossPoint[i].SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+    const double s = CrossPoint[i].ScalarProduct(&Direction)/Direction.NormSquared();
+    Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl;
+    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
+      CrossPoint[i].AddVector(endpoints[i%3]->node->node);  // make cross point absolute again
+      Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i%3]->node->node << " and " << *endpoints[(i+1)%3]->node->node << "." << endl;
+      const double distance = CrossPoint[i].DistanceSquared(x);
+      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
+        ShortestDistance = distance;
+        ClosestPoint->CopyVector(&CrossPoint[i]);
+      }
+    } else
+      CrossPoint[i].Zero();
+  }
+  InsideFlag = true;
+  for (int i=0;i<3;i++) {
+    const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+1)%3]);
+    const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+2)%3]);;
+    if ((sign > -MYEPSILON) && (othersign > -MYEPSILON))  // have different sign
+      InsideFlag = false;
+  }
+  if (InsideFlag) {
+    ClosestPoint->CopyVector(&InPlane);
+    ShortestDistance = InPlane.DistanceSquared(x);
+  } else {  // also check endnodes
+    for (int i=0;i<3;i++) {
+      const double distance = x->DistanceSquared(endpoints[i]->node->node);
+      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
+        ShortestDistance = distance;
+        ClosestPoint->CopyVector(endpoints[i]->node->node);
+      }
+    }
+  }
+  Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl;
+  return ShortestDistance;
 };
 
@@ -440 +562 @@
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
+bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
 {
         Info FunctionInfo(__func__);
@@ -453 +575 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
         Info FunctionInfo(__func__);
@@ -466 +588 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
 {
         Info FunctionInfo(__func__);
@@ -479 +601 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
-{
-        Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+        Log() << Verbose(1) << "INFO: Checking " << Points[0] << ","  << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl;
   return (((endpoints[0] == Points[0])
             || (endpoints[0] == Points[1])
@@ -501 +624 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
 {
         Info FunctionInfo(__func__);
@@ -523 +646 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
 {
         Info FunctionInfo(__func__);
@@ -540 +663 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector *center)
+void BoundaryTriangleSet::GetCenter(Vector * const center) const
 {
         Info FunctionInfo(__func__);
@@ -547 +670 @@
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
+  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -1422 +1546 @@
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
@@ -1437 +1561 @@
     Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
@@ -1587 +1711 @@
   bool insertNewLine = true;
 
-  if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  if (FindLine != a->lines.end()) {
+    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
+
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
 
     for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
@@ -1915 +2040 @@
   double maxCoordinate[NDIM];
   BoundaryLineSet BaseLine;
-  Vector Oben;
   Vector helper;
   Vector Chord;
   Vector SearchDirection;
-
-  Oben.Zero();
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector SphereCenter;
+  Vector NormalVector;
+
+  NormalVector.Zero();
 
   for (i = 0; i < 3; i++) {
@@ -1955 +2083 @@
   BTS = NULL;
   for (int k=0;k<NDIM;k++) {
-    Oben.Zero();
-    Oben.x[k] = 1.;
+    NormalVector.Zero();
+    NormalVector.x[k] = 1.;
     BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
     Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
@@ -1963 +2091 @@
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
 
-    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, Oben, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     if (Temporary == NULL)  // have we found a second point?
       continue;
     BaseLine.endpoints[1] = new BoundaryPointSet(Temporary);
 
-    helper.CopyVector(BaseLine.endpoints[0]->node->node);
-    helper.SubtractVector(BaseLine.endpoints[1]->node->node);
-    helper.Normalize();
-    Oben.ProjectOntoPlane(&helper);
-    Oben.Normalize();
-    helper.VectorProduct(&Oben);
+    // construct center of circle
+    CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node);
+    CircleCenter.AddVector(BaseLine.endpoints[1]->node->node);
+    CircleCenter.Scale(0.5);
+
+    // construct normal vector of circle
+    CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node);
+    CirclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node);
+
+    double radius = CirclePlaneNormal.NormSquared();
+    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+
+    NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
+    NormalVector.Normalize();
     ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
 
-    Chord.CopyVector(BaseLine.endpoints[0]->node->node); // bring into calling function
-    Chord.SubtractVector(BaseLine.endpoints[1]->node->node);
-    double radius = Chord.ScalarProduct(&Chord);
-    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
-    helper.CopyVector(&Oben);
-    helper.Scale(CircleRadius);
-    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+    SphereCenter.CopyVector(&NormalVector);
+    SphereCenter.Scale(CircleRadius);
+    SphereCenter.AddVector(&CircleCenter);
+    // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
@@ -1993 +2126 @@
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
     CandidateForTesselation OptCandidates(&BaseLine);
-    FindThirdPointForTesselation(Oben, SearchDirection, helper, OptCandidates, NULL, RADIUS, LC);
+    FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
     Log() << Verbose(0) << "List of third Points is:" << endl;
     for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) {
@@ -2167 +2300 @@
   Vector CircleCenter;
   Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
+  Vector RelativeSphereCenter;
   Vector SearchDirection;
   Vector helper;
@@ -2174 +2307 @@
   double radius, CircleRadius;
 
-  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " of triangle " << T << "." << endl;
   for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node))
+    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) {
       ThirdNode = T.endpoints[i]->node;
+      break;
+    }
+  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdNode " << *ThirdNode << " of triangle " << T << "." << endl;
 
   // construct center of circle
@@ -2191 +2326 @@
   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
   if (radius/4. < RADIUS*RADIUS) {
+    // construct relative sphere center with now known CircleCenter
+    RelativeSphereCenter.CopyVector(&T.SphereCenter);
+    RelativeSphereCenter.SubtractVector(&CircleCenter);
+
     CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
     Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
-    // construct old center
-    GetCenterofCircumcircle(&OldSphereCenter, *T.endpoints[0]->node->node, *T.endpoints[1]->node->node, *T.endpoints[2]->node->node);
-    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
-    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
-    helper.Scale(sqrt(RADIUS*RADIUS - radius));
-    OldSphereCenter.AddVector(&helper);
-    OldSphereCenter.SubtractVector(&CircleCenter);
-    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-
-    // construct SearchDirection
-    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl;
+
+    // construct SearchDirection and an "outward pointer"
+    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
+    helper.CopyVector(&CircleCenter);
     helper.SubtractVector(ThirdNode->node);
     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
     Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
@@ -2219 +2349 @@
 
     // add third point
-    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
 
   } else {
@@ -2332 +2462 @@
 
   // fill the set of neighbours
-  Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  Center.SubtractVector(TurningPoint->node);
-  set<TesselPoint*> SetOfNeighbours;
+  TesselPointSet SetOfNeighbours;
   SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
-  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
+  Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl;
+  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
+    Log() << Verbose(0) << **TesselRunner << endl;
   TesselPointList::iterator Runner = connectedClosestPoints->begin();
   TesselPointList::iterator Sprinter = Runner;
@@ -2350 +2481 @@
     AddTesselationPoint((*Sprinter), 2);
 
-    Center.CopyVector(&CandidateLine.OptCenter);
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -2359 +2489 @@
     BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     AddTesselationTriangle();
-    Center.Scale(-1.);
+    BTS->GetCenter(&Center);
+    Center.SubtractVector(&CandidateLine.OptCenter);
+    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
     BTS->GetNormalVector(Center);
 
@@ -2365 +2497 @@
     Runner = Sprinter;
     Sprinter++;
+    Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl;
+    if (Sprinter != connectedClosestPoints->end())
+      Log() << Verbose(0) << " There are still more triangles to add." << endl;
   }
   delete(connectedClosestPoints);
@@ -2766 +2901 @@
   Vector NewNormalVector;   // normal vector of the Candidate's triangle
   Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
+  Vector RelativeOldSphereCenter;
+  Vector NewPlaneCenter;
   double CircleRadius; // radius of this circle
   double radius;
+  double otherradius;
   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
@@ -2783 +2921 @@
   CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
 
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+
   // calculate squared radius TesselPoint *ThirdNode,f circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
+  radius = CirclePlaneNormal.NormSquared()/4.;
+  if (radius < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius;
     CirclePlaneNormal.Normalize();
-    //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
     // test whether old center is on the band's plane
-    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      eLog() << Verbose(1) << "Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
-      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
-    }
-    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+      Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl;
 
       // check SearchDirection
-      //Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
         eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
       }
@@ -2832 +2973 @@
 
                 // check for three unique points
-                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << *(Candidate->node) << "." << endl;
+                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl;
                 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){
 
-                  // construct both new centers
-                  GetCenterofCircumcircle(&NewSphereCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
-                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
-
-                  if ((NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node))
-                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                  // find center on the plane
+                  GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
+                  Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;
+
+                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)
+                  && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
                   ) {
-                    helper.CopyVector(&NewNormalVector);
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
+                    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+                    Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+                    Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl;
                     if (radius < RADIUS*RADIUS) {
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
+                      if (fabs(radius - otherradius) > HULLEPSILON) {
+                        eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl;
+                      }
+                      // construct both new centers
+                      NewSphereCenter.CopyVector(&NewPlaneCenter);
+                      OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
+                      helper.CopyVector(&NewNormalVector);
                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
-                      Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
+                      Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl;
                       NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-
                       // OtherNewSphereCenter is created by the same vector just in the other direction
                       helper.Scale(-1.);
                       OtherNewSphereCenter.AddVector(&helper);
-                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
 
@@ -2861 +3009 @@
                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                       alpha = min(alpha, Otheralpha);
+
                       // if there is a better candidate, drop the current list and add the new candidate
                       // otherwise ignore the new candidate and keep the list
@@ -2892 +3041 @@
                         }
                       }
-
                     } else {
                       Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
@@ -2936 +3084 @@
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  map<int, class BoundaryPointSet *> OrderMap;
-  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  PointMap OrderMap;
+  PointTestPair OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
@@ -2957 +3105 @@
 };
 
+/** Finds the boundary points that are closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
+ */
+DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+  PointMap::const_iterator FindPoint;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  if (LinesOnBoundary.empty()) {
+    eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl;
+    return NULL;
+  }
+
+  // gather all points close to the desired one
+  LC->SetIndexToVector(x); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  DistanceToPointMap * points = new DistanceToPointMap;
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //Log() << Verbose(1) << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        const LinkedNodes *List = LC->GetCurrentCell();
+        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            FindPoint = PointsOnBoundary.find((*Runner)->nr);
+            if (FindPoint != PointsOnBoundary.end()) {
+              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
+            }
+          }
+        } else {
+          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  // check whether we found some points
+  if (points->empty()) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    delete(points);
+    return NULL;
+  }
+  return points;
+};
+
+/** Finds the boundary line that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return closest BoundaryLineSet or NULL in degenerate case.
+ */
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+
+  // get closest points
+  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    return NULL;
+  }
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl;
+  BoundaryLineSet *ClosestLine = NULL;
+  double MinDistance = -1.;
+  Vector helper;
+  Vector Center;
+  Vector BaseLine;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+      // calculate closest point on line to desired point
+      helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
+      helper.Scale(0.5);
+      Center.CopyVector(x);
+      Center.SubtractVector(&helper);
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      Center.ProjectOntoPlane(&BaseLine);
+      const double distance = Center.NormSquared();
+      if ((ClosestLine == NULL) || (distance < MinDistance)) {
+        // additionally calculate intersection on line (whether it's on the line section or not)
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthA = helper.ScalarProduct(&BaseLine);
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthB = helper.ScalarProduct(&BaseLine);
+        if (lengthB*lengthA < 0) {  // if have different sign
+          ClosestLine = LineRunner->second;
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
+        } else {
+          Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
+        }
+      } else {
+        Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
+      }
+    }
+  }
+  delete(points);
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLine == NULL) {
+    Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
+    return NULL;
+  }
+  return ClosestLine;
+};
+
+
 /** Finds the triangle that is closest to a given Vector \a *x.
  * \param *out output stream for debugging
  * \param *x Vector to look from
- * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
- */
-list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const
-{
-        Info FunctionInfo(__func__);
-  TesselPoint *trianglePoints[3];
-  TesselPoint *SecondPoint = NULL;
-  list<BoundaryTriangleSet*> *triangles = NULL;
-
-  if (LinesOnBoundary.empty()) {
-    eLog() << Verbose(1) << "Error: There is no tesselation structure to compare the point with, please create one first.";
+ * \return BoundaryTriangleSet of nearest triangle or NULL.
+ */
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+{
+        Info FunctionInfo(__func__);
+
+        // get closest points
+        DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
     return NULL;
   }
-  Log() << Verbose(1) << "Finding closest Tesselpoint to " << *x << " ... " << endl;
-  trianglePoints[0] = FindClosestPoint(x, SecondPoint, LC);
-  
-  // check whether closest point is "too close" :), then it's inside
-  if (trianglePoints[0] == NULL) {
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
+  LineSet ClosestLines;
+  double MinDistance = 1e+16;
+  Vector BaseLineIntersection;
+  Vector Center;
+  Vector BaseLine;
+  Vector BaseLineCenter;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthBase = BaseLine.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+      const double lengthEndA = BaseLineIntersection.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthEndB = BaseLineIntersection.NormSquared();
+
+      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
+        const double lengthEnd = Min(lengthEndA, lengthEndB);
+        if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
+          ClosestLines.clear();
+          ClosestLines.insert(LineRunner->second);
+          MinDistance = lengthEnd;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl;
+        } else if  (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
+          ClosestLines.insert(LineRunner->second);
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl;
+        } else { // line is worse
+          Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
+        }
+      } else { // intersection is closer, calculate
+        // calculate closest point on line to desired point
+        BaseLineIntersection.CopyVector(x);
+        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        Center.CopyVector(&BaseLineIntersection);
+        Center.ProjectOntoPlane(&BaseLine);
+        BaseLineIntersection.SubtractVector(&Center);
+        const double distance = BaseLineIntersection.NormSquared();
+        if (Center.NormSquared() > BaseLine.NormSquared()) {
+          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
+        }
+        if ((ClosestLines.empty()) || (distance < MinDistance)) {
+          ClosestLines.insert(LineRunner->second);
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl;
+        } else {
+          Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
+        }
+      }
+    }
+  }
+  delete(points);
+
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLines.empty()) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
-  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
-    Log() << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
-    PointMap::const_iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
-    triangles = new list<BoundaryTriangleSet*>;
-    if (PointRunner != PointsOnBoundary.end()) {
-      for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
-        for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
-          triangles->push_back(TriangleRunner->second);
-      triangles->sort();
-      triangles->unique();
-    } else {
-      PointRunner = PointsOnBoundary.find(SecondPoint->nr);
-      trianglePoints[0] = SecondPoint;
-      if (PointRunner != PointsOnBoundary.end()) {
-        for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
-          for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
-            triangles->push_back(TriangleRunner->second);
-        triangles->sort();
-        triangles->unique();
-      } else {
-        eLog() << Verbose(1) << "I cannot find a boundary point to the tessel point " << *trianglePoints[0] << "." << endl;
-        return NULL;
-      }
-    }
-  } else {
-    set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(trianglePoints[0]);
-    TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(connectedPoints, trianglePoints[0], x);
-    delete(connectedPoints);
-    if (connectedClosestPoints != NULL) {
-      trianglePoints[1] = connectedClosestPoints->front();
-      trianglePoints[2] = connectedClosestPoints->back();
-      for (int i=0;i<3;i++) {
-        if (trianglePoints[i] == NULL) {
-          eLog() << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
-        }
-        //Log() << Verbose(1) << "List of triangle points:" << endl;
-        //Log() << Verbose(2) << *trianglePoints[i] << endl;
-      }
-
-      triangles = FindTriangles(trianglePoints);
-      Log() << Verbose(1) << "List of possible triangles:" << endl;
-      for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-        Log() << Verbose(2) << **Runner << endl;
-
-      delete(connectedClosestPoints);
-    } else {
-      triangles = NULL;
-      eLog() << Verbose(2) << "There is no circle of connected points!" << endl;
-    }
-  }
-  
-  if ((triangles == NULL) || (triangles->empty())) {
-    eLog() << Verbose(1) << "There is no nearest triangle. Please check the tesselation structure.";
-    delete(triangles);
-    return NULL;
-  } else
-    return triangles;
+  TriangleList * candidates = new TriangleList;
+  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
+    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
+    candidates->push_back(Runner->second);
+  }
+  return candidates;
 };
 
@@ -3046 +3318 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
 {
         Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
-  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
+  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
+  TriangleList candidates;
   Vector Center;
-
-  if (triangles == NULL)
+  Vector helper;
+
+  if ((triangles == NULL) || (triangles->empty()))
     return NULL;
 
-  if (triangles->size() == 1) { // there is no degenerate case
-    result = triangles->front();
-    Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
-  } else {
-    result = triangles->front();
-    result->GetCenter(&Center);
-    Center.SubtractVector(x);
-    Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
-    if (Center.ScalarProduct(&result->NormalVector) < 0) {
-      result = triangles->back();
-      Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
-      if (Center.ScalarProduct(&result->NormalVector) < 0) {
-        eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
-      }
+  // go through all and pick the one with the best alignment to x
+  double MinAlignment = 2.*M_PI;
+  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
+    (*Runner)->GetCenter(&Center);
+    helper.CopyVector(x);
+    helper.SubtractVector(&Center);
+    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
+    if (Alignment < MinAlignment) {
+      result = *Runner;
+      MinAlignment = Alignment;
+      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
+    } else {
+      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
     }
   }
   delete(triangles);
+
   return result;
 };
 
-/** Checks whether the provided Vector is within the tesselation structure.
+/** Checks whether the provided Vector is within the Tesselation structure.
+ * Basically calls Tesselation::GetDistanceToSurface() and checks the sign of the return value.
+ * @param point of which to check the position
+ * @param *LC LinkedCell structure
+ *
+ * @return true if the point is inside the Tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
+{
+  return (GetDistanceSquaredToSurface(Point, LC) < MYEPSILON);
+}
+
+/** Returns the distance to the surface given by the tesselation.
+ * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
+ * towards or away from the given \a &Point. Additionally, we check whether it's normal to the normal vector, i.e. on the
+ * closest triangle's plane. Then, we have to check whether \a Point is inside the triangle or not to determine whether it's
+ * an inside or outside point. This is done by calling BoundaryTriangleSet::GetIntersectionInsideTriangle().
+ * In the end we additionally find the point on the triangle who was smallest distance to \a Point:
+ *  -# Separate distance from point to center in vector in NormalDirection and on the triangle plane.
+ *  -# Check whether vector on triangle plane points inside the triangle or crosses triangle bounds.
+ *  -# If inside, take it to calculate closest distance
+ *  -# If not, take intersection with BoundaryLine as distance
+ *
+ * @note distance is squared despite it still contains a sign to determine in-/outside!
  *
  * @param point of which to check the position
  * @param *LC LinkedCell structure
  *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
-{
-        Info FunctionInfo(__func__);
-  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
+ * @return >0 if outside, ==0 if on surface, <0 if inside (Note that distance can be at most LinkedCell::RADIUS.)
+ */
+double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
+{
+  Info FunctionInfo(__func__);
+  class BoundaryTriangleSet *result = FindClosestTriangleToVector(&Point, LC);
   Vector Center;
+  Vector helper;
+  Vector DistanceToCenter;
+  Vector Intersection;
+  double distance = 0.;
 
   if (result == NULL) {// is boundary point or only point in point cloud?
     Log() << Verbose(1) << Point << " is the only point in vicinity." << endl;
-    return false;
+    return LC->RADIUS;
+  } else {
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
   }
 
   result->GetCenter(&Center);
   Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
-  Center.SubtractVector(&Point);
-  Log() << Verbose(2) << "INFO: Vector from center to point to test is " << Center << "." << endl;
-  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
-    Log() << Verbose(1) << Point << " is an inner point." << endl;
-    return true;
+  DistanceToCenter.CopyVector(&Center);
+  DistanceToCenter.SubtractVector(&Point);
+  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl;
+
+  // check whether we are on boundary
+  if (fabs(DistanceToCenter.ScalarProduct(&result->NormalVector)) < MYEPSILON) {
+    // calculate whether inside of triangle
+    DistanceToCenter.CopyVector(&Point);
+    Center.CopyVector(&Point);
+    Center.SubtractVector(&result->NormalVector); // points towards MolCenter
+    DistanceToCenter.AddVector(&result->NormalVector); // points outside
+    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
+    if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
+      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
+      return 0.;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
+      return false;
+    }
   } else {
-    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
-    return false;
-  }
-}
-
-/** Checks whether the provided TesselPoint is within the tesselation structure.
- *
- * @param *Point of which to check the position
- * @param *LC Linked Cell structure
- *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const
-{
-        Info FunctionInfo(__func__);
-  return IsInnerPoint(*(Point->node), LC);
-}
+    // calculate smallest distance
+    distance = result->GetClosestPointInsideTriangle(&Point, &Intersection);
+    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
+    distance = Min(distance, (LC->RADIUS*LC->RADIUS));
+
+    // then check direction to boundary
+    if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
+      Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
+      return -distance;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl;
+      return +distance;
+    }
+  }
+};
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -3126 +3441 @@
  * @return set of the all points linked to the provided one
  */
-set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
-{
-        Info FunctionInfo(__func__);
-  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
+TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+        TesselPointSet *connectedPoints = new TesselPointSet;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -3170 +3485 @@
   }
 
-  if (connectedPoints->size() == 0) { // if have not found any points
+  if (connectedPoints->empty()) { // if have not found any points
     eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl;
     return NULL;
@@ -3191 +3506 @@
  * @return list of the all points linked to the provided one
  */
-list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
-  Vector center;
+  TesselPointList *connectedCircle = new TesselPointList;
   Vector PlaneNormal;
   Vector AngleZero;
   Vector OrthogonalVector;
   Vector helper;
+  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
+  TriangleList *triangles = NULL;
 
   if (SetOfNeighbours == NULL) {
@@ -3209 +3525 @@
 
   // calculate central point
-  for (set<TesselPoint*>::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
-    center.AddVector((*TesselRunner)->node);
-  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
-  //  << "; scale factor " << 1.0/connectedPoints.size();
-  center.Scale(1.0/SetOfNeighbours->size());
-  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
-
-  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
-  PlaneNormal.CopyVector(Point->node);
-  PlaneNormal.SubtractVector(&center);
+  triangles = FindTriangles(TrianglePoints);
+  if ((triangles != NULL) && (!triangles->empty())) {
+    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      PlaneNormal.AddVector(&(*Runner)->NormalVector);
+  } else {
+    eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl;
+    performCriticalExit();
+  }
+  PlaneNormal.Scale(1.0/triangles->size());
+  Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl;
   PlaneNormal.Normalize();
-  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
 
   // construct one orthogonal vector
@@ -3246 +3561 @@
 
   // go through all connected points and calculate angle
-  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
@@ -3262 +3577 @@
 }
 
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
+ * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
+ * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
+ * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
+ * triangle we are looking for.
+ *
+ * @param *SetOfNeighbours all points for which the angle should be calculated
+ * @param *Point of which get all connected points
+ * @param *Reference Reference vector for zero angle or NULL for no preference
+ * @return list of the all points linked to the provided one
+ */
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+{
+  Info FunctionInfo(__func__);
+  map<double, TesselPoint*> anglesOfPoints;
+  TesselPointList *connectedCircle = new TesselPointList;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+
+  if (SetOfNeighbours == NULL) {
+    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
+    delete(connectedCircle);
+    return NULL;
+  }
+
+  // check whether there's something to do
+  if (SetOfNeighbours->size() < 3) {
+    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+      connectedCircle->push_back(*TesselRunner);
+    return connectedCircle;
+  }
+
+  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
+  // calculate central point
+
+  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
+  TesselB++;
+  TesselC++;
+  TesselC++;
+  int counter = 0;
+  while (TesselC != SetOfNeighbours->end()) {
+    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
+    counter++;
+    TesselA++;
+    TesselB++;
+    TesselC++;
+    PlaneNormal.AddVector(&helper);
+  }
+  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
+  //  << "; scale factor " << counter;
+  PlaneNormal.Scale(1.0/(double)counter);
+//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
+//
+//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
+//  PlaneNormal.CopyVector(Point->node);
+//  PlaneNormal.SubtractVector(&center);
+//  PlaneNormal.Normalize();
+  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
+
+  // construct one orthogonal vector
+  if (Reference != NULL) {
+    AngleZero.CopyVector(Reference);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+  }
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
+    Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;
+    AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
+      performCriticalExit();
+    }
+  }
+  Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
+  if (AngleZero.NormSquared() > MYEPSILON)
+    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+  else
+    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+  Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
+
+  // go through all connected points and calculate angle
+  pair <map<double, TesselPoint*>::iterator, bool > InserterTest;
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+    helper.CopyVector((*listRunner)->node);
+    helper.SubtractVector(Point->node);
+    helper.ProjectOntoPlane(&PlaneNormal);
+    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
+    if (angle > M_PI) // the correction is of no use here (and not desired)
+      angle = 2.*M_PI - angle;
+    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    if (!InserterTest.second) {
+      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
+      performCriticalExit();
+    }
+  }
+
+  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
+    connectedCircle->push_back(AngleRunner->second);
+  }
+
+  return connectedCircle;
+}
+
 /** Gets all points connected to the provided point by triangulation lines, ordered such that we walk along a closed path.
  *
@@ -3268 +3695 @@
  * @return list of the all points linked to the provided one
  */
-list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+list< TesselPointList *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
+  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
@@ -3310 +3737 @@
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new list<TesselPoint*>;
+        connectedPath = new TesselPointList;
         triangle = NULL;
         CurrentLine = runner->second;
@@ -3384 +3811 @@
  * @return list of the closed paths
  */
-list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
-        Info FunctionInfo(__func__);
-  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point);
-  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
-  list<TesselPoint*> *newPath = NULL;
+list<TesselPointList *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+  list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
+  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
+  TesselPointList *newPath = NULL;
   int count = 0;
 
 
-  list<TesselPoint*>::iterator CircleRunner;
-  list<TesselPoint*>::iterator CircleStart;
-
-  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  TesselPointList::iterator CircleRunner;
+  TesselPointList::iterator CircleStart;
+
+  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
@@ -3406 +3833 @@
 
     // go through list, look for reappearance of starting Point and create list
-    list<TesselPoint*>::iterator Marker = CircleStart;
+    TesselPointList::iterator Marker = CircleStart;
     for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
         Log() << Verbose(1) << count+1 << ". closed path consists of: ";
-        newPath = new list<TesselPoint*>;
-        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        newPath = new TesselPointList;
+        TesselPointList::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
@@ -3446 +3873 @@
  * \return pointer to allocated list of triangles
  */
-set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
-{
-        Info FunctionInfo(__func__);
-  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+        Info FunctionInfo(__func__);
+        TriangleSet *connectedTriangles = new TriangleSet;
 
   if (Point == NULL) {
@@ -3498 +3925 @@
   }
 
-  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
 
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  map<class BoundaryTriangleSet *, int> Candidates;
+  TriangleMap Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
     }
   }
@@ -3516 +3943 @@
   count=0;
   NormalVector.Zero();
-  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
-    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
-    RemoveTesselationTriangle(Runner->first);
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl;
+    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->second);
     count++;
   }
   Log() << Verbose(1) << count << " triangles were removed." << endl;
 
-  list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
-  list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
-  map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
-  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+  TriangleMap::iterator NumberRunner = Candidates.begin();
+  TesselPointList::iterator StartNode, MiddleNode, EndNode;
   double angle;
   double smallestangle;
@@ -3541 +3968 @@
 
       // re-create all triangles by going through connected points list
-      list<class BoundaryLineSet *> NewLines;
+      LineList NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
@@ -3647 +4074 @@
       // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
       if (NewLines.size() > 1) {
-        list<class BoundaryLineSet *>::iterator Candidate;
+        LineList::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
-          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -3694 +4121 @@
  * Finds triangles belonging to the three provided points.
  *
- * @param *Points[3] list, is expected to contain three points
+ * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
  *
  * @return triangles which belong to the provided points, will be empty if there are none,
  *         will usually be one, in case of degeneration, there will be two
  */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
-{
-        Info FunctionInfo(__func__);
-  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
+TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+        TriangleList *result = new TriangleList;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
   class BoundaryPointSet *TrianglePoints[3];
+  size_t NoOfWildcards = 0;
 
   for (int i = 0; i < 3; i++) {
-    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      TrianglePoints[i] = FindPoint->second;
+    if (Points[i] == NULL) {
+      NoOfWildcards++;
+      TrianglePoints[i] = NULL;
     } else {
-      TrianglePoints[i] = NULL;
-    }
-  }
-
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (TrianglePoints[i] != NULL) {
-      for (int j = i+1; j < 3; j++) {
-        if (TrianglePoints[j] != NULL) {
-          for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
-              (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
-              FindLine++) {
-            for (FindTriangle = FindLine->second->triangles.begin();
-                FindTriangle != FindLine->second->triangles.end();
-                FindTriangle++) {
-              if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
-                result->push_back(FindTriangle->second);
+      PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
+      if (FindPoint != PointsOnBoundary.end()) {
+        TrianglePoints[i] = FindPoint->second;
+      } else {
+        TrianglePoints[i] = NULL;
+      }
+    }
+  }
+
+  switch (NoOfWildcards) {
+    case 0: // checks lines between the points in the Points for their adjacent triangles
+      for (int i = 0; i < 3; i++) {
+        if (TrianglePoints[i] != NULL) {
+          for (int j = i+1; j < 3; j++) {
+            if (TrianglePoints[j] != NULL) {
+              for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
+                  (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
+                  FindLine++) {
+                for (FindTriangle = FindLine->second->triangles.begin();
+                    FindTriangle != FindLine->second->triangles.end();
+                    FindTriangle++) {
+                  if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+                    result->push_back(FindTriangle->second);
+                  }
+                }
               }
+              // Is it sufficient to consider one of the triangle lines for this.
+              return result;
             }
           }
-          // Is it sufficient to consider one of the triangle lines for this.
-          return result;
         }
       }
-    }
+      break;
+    case 1: // copy all triangles of the respective line
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] == NULL)
+          break;
+      for (FindLine = TrianglePoints[(i+1)%3]->lines.find(TrianglePoints[(i+2)%3]->node->nr); // is a multimap!
+          (FindLine != TrianglePoints[(i+1)%3]->lines.end()) && (FindLine->first == TrianglePoints[(i+2)%3]->node->nr);
+          FindLine++) {
+        for (FindTriangle = FindLine->second->triangles.begin();
+            FindTriangle != FindLine->second->triangles.end();
+            FindTriangle++) {
+          if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+            result->push_back(FindTriangle->second);
+          }
+        }
+      }
+      break;
+    }
+    case 2: // copy all triangles of the respective point
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] != NULL)
+          break;
+      for (LineMap::const_iterator line = TrianglePoints[i]->lines.begin(); line != TrianglePoints[i]->lines.end(); line++)
+        for (TriangleMap::const_iterator triangle = line->second->triangles.begin(); triangle != line->second->triangles.end(); triangle++)
+          result->push_back(triangle->second);
+      result->sort();
+      result->unique();
+      break;
+    }
+    case 3: // copy all triangles
+    {
+      for (TriangleMap::const_iterator triangle = TrianglesOnBoundary.begin(); triangle != TrianglesOnBoundary.end(); triangle++)
+        result->push_back(triangle->second);
+      break;
+    }
+    default:
+      eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl;
+      performCriticalExit();
+      break;
   }
 
@@ -3779 +4257 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedLines()
+IndexToIndex * Tesselation::FindAllDegeneratedLines()
 {
         Info FunctionInfo(__func__);
         UniqueLines AllLines;
-  map<int, int> * DegeneratedLines = new map<int, int>;
+  IndexToIndex * DegeneratedLines = new IndexToIndex;
 
   // sanity check
@@ -3804 +4282 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
     const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
@@ -3823 +4301 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedTriangles()
-{
-        Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
-  map<int, int> * DegeneratedTriangles = new map<int, int>;
+IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
+{
+        Info FunctionInfo(__func__);
+  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
+  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
 
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
@@ -3833 +4311 @@
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
 
-  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
@@ -3854 +4332 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
@@ -3868 +4346 @@
 {
         Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
 
-  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
@@ -3961 +4439 @@
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
@@ -4128 +4606 @@
 
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
+  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
@@ -4138 +4617 @@
     for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++)
       for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
-        TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
-        if (TriangleInsertionTester.second)
-          Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {
+          TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
+          if (TriangleInsertionTester.second)
+            Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        } else {
+          Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl;
+        }
       }
     // check whether there are two that are parallel
@@ -4213 +4696 @@
     /// 4a. Gather all triangles of this polygon
     TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();
+
+    // check whether number is bigger than 2, otherwise it's just a simply degenerated one and nothing to do.
+    if (T->size() == 2) {
+      Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl;
+      delete(T);
+      continue;
+    }
+
+    // check whether number is even
+    // If this case occurs, we have to think about it!
+    // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
+    // connections to either polygon ...
+    if (T->size() % 2 != 0) {
+      eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl;
+      performCriticalExit();
+    }
 
     TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator
@@ -4259 +4758 @@
   }
 
-  map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
   Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;