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

File:

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

molecuilder/src/tesselation.cpp

@@ -35 +35 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
   LinesCount(0),
   node(Walker),
@@ -61 +61 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -105 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *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;
@@ -193 +171 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
         Info FunctionInfo(__func__);
@@ -204 +182 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -219 +197 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion() const
+bool BoundaryLineSet::CheckConvexityCriterion()
 {
         Info FunctionInfo(__func__);
@@ -243 +221 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for(TriangleMap::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);
@@ -286 +264 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -299 +277 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -339 +317 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
   Nr(number)
 {
@@ -398 +376 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
 {
         Info FunctionInfo(__func__);
@@ -410 +388 @@
 };
 
-/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
+/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
  * 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.
+ * This 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
@@ -422 +400 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
-{
-  Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
+{
+        Info FunctionInfo(__func__);
   Vector CrossPoint;
   Vector helper;
@@ -433 +411 @@
   }
 
-  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;
@@ -453 +417 @@
       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;
-      }
+    } else 
       i++;
-    } else 
+    if (i>2)
       break;
-  } while (i<3);
+  } while (CrossPoint.NormSquared() < MYEPSILON);
   if (i==3) {
-    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
+    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
     return true;
-  } else {
-    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
+  } else { // outside!
+    Intersection->Zero();
     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;
 };
 
@@ -562 +440 @@
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -575 +453 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -588 +466 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
 {
         Info FunctionInfo(__func__);
@@ -601 +479 @@
  * \return true - is the very triangle, false - is not
  */
-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;
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
+{
+        Info FunctionInfo(__func__);
   return (((endpoints[0] == Points[0])
             || (endpoints[0] == Points[1])
@@ -624 +501 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
 {
         Info FunctionInfo(__func__);
@@ -646 +523 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -663 +540 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector * const center) const
+void BoundaryTriangleSet::GetCenter(Vector *center)
 {
         Info FunctionInfo(__func__);
@@ -670 +547 @@
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
-  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -1546 +1422 @@
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  TriangleList *triangles = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
@@ -1561 +1437 @@
     Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
@@ -1711 +1587 @@
   bool insertNewLine = true;
 
-  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;
-
+  if (a->lines.find(b->node->nr) != a->lines.end()) {
+    LineMap::iterator FindLine = a->lines.find(b->node->nr);
     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++) {
@@ -2040 +1915 @@
   double maxCoordinate[NDIM];
   BoundaryLineSet BaseLine;
+  Vector Oben;
   Vector helper;
   Vector Chord;
   Vector SearchDirection;
-  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();
+
+  Oben.Zero();
 
   for (i = 0; i < 3; i++) {
@@ -2083 +1955 @@
   BTS = NULL;
   for (int k=0;k<NDIM;k++) {
-    NormalVector.Zero();
-    NormalVector.x[k] = 1.;
+    Oben.Zero();
+    Oben.x[k] = 1.;
     BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
     Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
@@ -2091 +1963 @@
     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, NormalVector, 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, Oben, 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);
 
-    // 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();
+    helper.CopyVector(BaseLine.endpoints[0]->node->node);
+    helper.SubtractVector(BaseLine.endpoints[1]->node->node);
+    helper.Normalize();
+    Oben.ProjectOntoPlane(&helper);
+    Oben.Normalize();
+    helper.VectorProduct(&Oben);
+    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.);
-
-    NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
-    NormalVector.Normalize();
-    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-
-    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)
+    helper.CopyVector(&Oben);
+    helper.Scale(CircleRadius);
+    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
@@ -2126 +1993 @@
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
     CandidateForTesselation OptCandidates(&BaseLine);
-    FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
+    FindThirdPointForTesselation(Oben, SearchDirection, helper, 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++) {
@@ -2300 +2167 @@
   Vector CircleCenter;
   Vector CirclePlaneNormal;
-  Vector RelativeSphereCenter;
+  Vector OldSphereCenter;
   Vector SearchDirection;
   Vector helper;
@@ -2307 +2174 @@
   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
@@ -2326 +2191 @@
   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;
 
-    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl;
-
-    // construct SearchDirection and an "outward pointer"
-    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
-    helper.CopyVector(&CircleCenter);
+    // 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);
     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(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
@@ -2349 +2219 @@
 
     // add third point
-    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
 
   } else {
@@ -2462 +2332 @@
 
   // fill the set of neighbours
-  TesselPointSet SetOfNeighbours;
+  Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);
+  Center.SubtractVector(TurningPoint->node);
+  set<TesselPoint*> 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, CandidateLine.BaseLine->endpoints[1]->node->node);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
 
   // 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;
@@ -2481 +2350 @@
     AddTesselationPoint((*Sprinter), 2);
 
+    Center.CopyVector(&CandidateLine.OptCenter);
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -2489 +2359 @@
     BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     AddTesselationTriangle();
-    BTS->GetCenter(&Center);
-    Center.SubtractVector(&CandidateLine.OptCenter);
-    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+    Center.Scale(-1.);
     BTS->GetNormalVector(Center);
 
@@ -2497 +2365 @@
     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);
@@ -2901 +2766 @@
   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];
@@ -2921 +2783 @@
   CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
 
-  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
-  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
-
   // calculate squared radius TesselPoint *ThirdNode,f circle
-  radius = CirclePlaneNormal.NormSquared()/4.;
-  if (radius < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius;
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
     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(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(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(radius - CircleRadius) < HULLEPSILON) {
-      Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl;
+      //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 
       // check SearchDirection
-      Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      //Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
         eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
       }
@@ -2973 +2832 @@
 
                 // check for three unique points
-                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl;
+                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << *(Candidate->node) << "." << endl;
                 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){
 
-                  // 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)
+                  // 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)
                   ) {
+                    helper.CopyVector(&NewNormalVector);
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    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;
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
                     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 NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl;
+                      Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " 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;
 
@@ -3009 +2861 @@
                       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
@@ -3041 +2892 @@
                         }
                       }
+
                     } else {
                       Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
@@ -3084 +2936 @@
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  PointMap OrderMap;
-  PointTestPair OrderTest;
+  map<int, class BoundaryPointSet *> OrderMap;
+  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
@@ -3105 +2957 @@
 };
 
-/** 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 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 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 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(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) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
-  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;
+  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;
 };
 
@@ -3318 +3046 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
 {
         Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
-  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
-  TriangleList candidates;
+  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
   Vector Center;
-  Vector helper;
-
-  if ((triangles == NULL) || (triangles->empty()))
+
+  if (triangles == NULL)
     return NULL;
 
-  // 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;
+  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;
+      }
     }
   }
   delete(triangles);
-
   return result;
 };
 
-/** 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!
+/** Checks whether the provided Vector is within the tesselation structure.
  *
  * @param point of which to check the position
  * @param *LC LinkedCell structure
  *
- * @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);
+ * @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);
   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 LC->RADIUS;
-  } else {
-    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
+    return false;
   }
 
   result->GetCenter(&Center);
   Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
-  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;
-    }
+  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;
   } else {
-    // 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;
-    }
-  }
-};
+    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);
+}
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -3441 +3126 @@
  * @return set of the all points linked to the provided one
  */
-TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
-{
-        Info FunctionInfo(__func__);
-        TesselPointSet *connectedPoints = new TesselPointSet;
+set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -3485 +3170 @@
   }
 
-  if (connectedPoints->empty()) { // if have not found any points
+  if (connectedPoints->size() == 0) { // if have not found any points
     eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl;
     return NULL;
@@ -3506 +3191 @@
  * @return list of the all points linked to the provided one
  */
-TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  TesselPointList *connectedCircle = new TesselPointList;
+  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  Vector center;
   Vector PlaneNormal;
   Vector AngleZero;
   Vector OrthogonalVector;
   Vector helper;
-  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
-  TriangleList *triangles = NULL;
 
   if (SetOfNeighbours == NULL) {
@@ -3525 +3209 @@
 
   // calculate central point
-  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;
+  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);
   PlaneNormal.Normalize();
+  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
 
   // construct one orthogonal vector
@@ -3561 +3246 @@
 
   // go through all connected points and calculate angle
-  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
@@ -3577 +3262 @@
 }
 
-/** 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.
  *
@@ -3695 +3268 @@
  * @return list of the all points linked to the provided one
  */
-list< TesselPointList *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
-  TesselPointList *connectedPath = NULL;
+  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
@@ -3737 +3310 @@
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new TesselPointList;
+        connectedPath = new list<TesselPoint*>;
         triangle = NULL;
         CurrentLine = runner->second;
@@ -3811 +3384 @@
  * @return list of the closed paths
  */
-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;
+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;
   int count = 0;
 
 
-  TesselPointList::iterator CircleRunner;
-  TesselPointList::iterator CircleStart;
-
-  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  list<TesselPoint*>::iterator CircleRunner;
+  list<TesselPoint*>::iterator CircleStart;
+
+  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
@@ -3833 +3406 @@
 
     // go through list, look for reappearance of starting Point and create list
-    TesselPointList::iterator Marker = CircleStart;
+    list<TesselPoint*>::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 TesselPointList;
-        TesselPointList::iterator CircleSprinter = Marker;
+        newPath = new list<TesselPoint*>;
+        list<TesselPoint*>::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
@@ -3873 +3446 @@
  * \return pointer to allocated list of triangles
  */
-TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
-{
-        Info FunctionInfo(__func__);
-        TriangleSet *connectedTriangles = new TriangleSet;
+set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+        Info FunctionInfo(__func__);
+  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
 
   if (Point == NULL) {
@@ -3925 +3498 @@
   }
 
-  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  TesselPointList *connectedPath = NULL;
+  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  list<TesselPoint*> *connectedPath = NULL;
 
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  TriangleMap Candidates;
+  map<class BoundaryTriangleSet *, int> 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( TrianglePair (triangle->Nr, triangle) );
+      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
     }
   }
@@ -3943 +3516 @@
   count=0;
   NormalVector.Zero();
-  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);
+  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);
     count++;
   }
   Log() << Verbose(1) << count << " triangles were removed." << endl;
 
-  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
-  list<TesselPointList *>::iterator ListRunner = ListAdvance;
-  TriangleMap::iterator NumberRunner = Candidates.begin();
-  TesselPointList::iterator StartNode, MiddleNode, EndNode;
+  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;
   double angle;
   double smallestangle;
@@ -3968 +3541 @@
 
       // re-create all triangles by going through connected points list
-      LineList NewLines;
+      list<class BoundaryLineSet *> NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
@@ -4074 +3647 @@
       // 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) {
-        LineList::iterator Candidate;
+        list<class BoundaryLineSet *>::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
-          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -4121 +3694 @@
  * Finds triangles belonging to the three provided points.
  *
- * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
+ * @param *Points[3] list, is expected to contain three points
  *
  * @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
  */
-TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
-{
-        Info FunctionInfo(__func__);
-        TriangleList *result = new TriangleList;
+list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
   class BoundaryPointSet *TrianglePoints[3];
-  size_t NoOfWildcards = 0;
 
   for (int i = 0; i < 3; i++) {
-    if (Points[i] == NULL) {
-      NoOfWildcards++;
+    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      TrianglePoints[i] = FindPoint->second;
+    } else {
       TrianglePoints[i] = NULL;
-    } else {
-      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);
-                  }
-                }
+    }
+  }
+
+  // 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;
+    }
   }
 
@@ -4257 +3779 @@
  *         in the list, once as key and once as value
  */
-IndexToIndex * Tesselation::FindAllDegeneratedLines()
+map<int, int> * Tesselation::FindAllDegeneratedLines()
 {
         Info FunctionInfo(__func__);
         UniqueLines AllLines;
-  IndexToIndex * DegeneratedLines = new IndexToIndex;
+  map<int, int> * DegeneratedLines = new map<int, int>;
 
   // sanity check
@@ -4282 +3804 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
-  IndexToIndex::iterator it;
+  map<int,int>::iterator it;
   for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
     const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
@@ -4301 +3823 @@
  *         in the list, once as key and once as value
  */
-IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
-{
-        Info FunctionInfo(__func__);
-  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
-  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
+map<int, int> * Tesselation::FindAllDegeneratedTriangles()
+{
+        Info FunctionInfo(__func__);
+  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
+  map<int, int> * DegeneratedTriangles = new map<int, int>;
 
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
@@ -4311 +3833 @@
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
 
-  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
@@ -4332 +3854 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
-  IndexToIndex::iterator it;
+  map<int,int>::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
@@ -4346 +3868 @@
 {
         Info FunctionInfo(__func__);
-  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
 
-  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
@@ -4439 +3961 @@
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
@@ -4606 +4128 @@
 
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
-  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
@@ -4617 +4138 @@
     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++) {
-        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;
-        }
+        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;
       }
     // check whether there are two that are parallel
@@ -4696 +4213 @@
     /// 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
@@ -4758 +4259 @@
   }
 
-  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
   Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
-  IndexToIndex::iterator it;
+  map<int,int>::iterator it;
   for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;