Changes in src/tesselation.cpp [b1a6d8:27bd2f]

File:

• src/tesselation.cpp (modified) (175 diffs)

src/tesselation.cpp

@@ -7 +7 @@
 
 #include <fstream>
-#include <assert.h>
 
 #include "helpers.hpp"
@@ -15 +14 @@
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
-#include "triangleintersectionlist.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
@@ -26 +24 @@
  */
 BoundaryPointSet::BoundaryPointSet() :
-  LinesCount(0), value(0.), Nr(-1)
-{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (Log() << Verbose(1) << "Adding noname." << endl);
-}
-;
+    LinesCount(0),
+    value(0.),
+    Nr(-1)
+{
+        Info FunctionInfo(__func__);
+        Log() << Verbose(1) << "Adding noname." << endl;
+};
 
 /** Constructor of BoundaryPointSet with Tesselpoint.
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
-  LinesCount(0), node(Walker), value(0.), Nr(Walker->nr)
-{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (Log() << Verbose(1) << "Adding Node " << *Walker << endl);
-}
-;
+BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
+  LinesCount(0),
+  node(Walker),
+  value(0.),
+  Nr(Walker->nr)
+{
+        Info FunctionInfo(__func__);
+  Log() << Verbose(1) << "Adding Node " << *Walker << endl;
+};
 
 /** Destructor of BoundaryPointSet.
@@ -50 +51 @@
 BoundaryPointSet::~BoundaryPointSet()
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl;
   if (!lines.empty())
-    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);
+    eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl;
   node = NULL;
-}
-;
+};
 
 /** Add a line to the LineMap of this point.
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
-{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." << endl);
-  if (line->endpoints[0] == this) {
-    lines.insert(LinePair(line->endpoints[1]->Nr, line));
-  } else {
-    lines.insert(LinePair(line->endpoints[0]->Nr, line));
-  }
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+{
+        Info FunctionInfo(__func__);
+  Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "."
+      << endl;
+  if (line->endpoints[0] == this)
+    {
+      lines.insert(LinePair(line->endpoints[1]->Nr, line));
+    }
+  else
+    {
+      lines.insert(LinePair(line->endpoints[0]->Nr, line));
+    }
   LinesCount++;
-}
-;
+};
 
 /** output operator for BoundaryPointSet.
@@ -90 +93 @@
  */
 BoundaryLineSet::BoundaryLineSet() :
-  Nr(-1)
-{
-  Info FunctionInfo(__func__);
+    Nr(-1)
+{
+        Info FunctionInfo(__func__);
   for (int i = 0; i < 2; i++)
     endpoints[i] = NULL;
-}
-;
+};
 
 /** Constructor of BoundaryLineSet with two endpoints.
@@ -103 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
-{
-  Info FunctionInfo(__func__);
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
+{
+        Info FunctionInfo(__func__);
   // set number
   Nr = number;
@@ -116 +118 @@
   skipped = false;
   // clear triangles list
-  DoLog(0) && (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;
-  // clear triangles list
-  DoLog(0) && (Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl);
-}
-;
+  Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
+};
 
 /** Destructor for BoundaryLineSet.
@@ -149 +127 @@
 BoundaryLineSet::~BoundaryLineSet()
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   int Numbers[2];
 
@@ -180 +158 @@
         //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         if (endpoints[i] != NULL) {
-          delete (endpoints[i]);
+          delete(endpoints[i]);
           endpoints[i] = NULL;
         }
@@ -187 +165 @@
   }
   if (!triangles.empty())
-    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl);
-}
-;
+    eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+};
 
 /** Add triangle to TriangleMap of this boundary line.
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
-{
-  Info FunctionInfo(__func__);
-  DoLog(0) && (Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl);
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+{
+        Info FunctionInfo(__func__);
+  Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
   triangles.insert(TrianglePair(triangle->Nr, triangle));
-}
-;
+};
 
 /** Checks whether we have a common endpoint with given \a *line.
@@ -206 +182 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
-{
-  Info FunctionInfo(__func__);
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
+{
+        Info FunctionInfo(__func__);
   if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
     return true;
   else
     return false;
-}
-;
+};
 
 /** Checks whether the adjacent triangles of a baseline are convex or not.
@@ -222 +197 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion() const
-{
-  Info FunctionInfo(__func__);
+bool BoundaryLineSet::CheckConvexityCriterion()
+{
+        Info FunctionInfo(__func__);
   Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
   // get the two triangles
   if (triangles.size() != 2) {
-    DoeLog(0) && (eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);
+    eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;
     return true;
   }
@@ -236 +211 @@
   BaseLineCenter.CopyVector(endpoints[0]->node->node);
   BaseLineCenter.AddVector(endpoints[1]->node->node);
-  BaseLineCenter.Scale(1. / 2.);
+  BaseLineCenter.Scale(1./2.);
   BaseLine.CopyVector(endpoints[0]->node->node);
   BaseLine.SubtractVector(endpoints[1]->node->node);
@@ -244 +219 @@
   NormalCheck.Zero();
   double sign = -1.;
-  int i = 0;
+  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);
@@ -252 +227 @@
     sign = -sign;
     if (runner->second->NormalVector.NormSquared() > MYEPSILON)
-      BaseLineNormal.CopyVector(&runner->second->NormalVector); // yes, copy second on top of first
+      BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
     else {
-      DoeLog(0) && (eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
+      eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl;
     }
     node = runner->second->GetThirdEndpoint(this);
@@ -261 +236 @@
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
-      helper[i].MakeNormalVector(&BaseLine); // we want to compare the triangle's heights' angles!
+      helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
       //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
     } else {
-      DoeLog(1) && (eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);
+      eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl;
       return true;
     }
@@ -271 +246 @@
   //Log() << Verbose(0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
   if (NormalCheck.NormSquared() < MYEPSILON) {
-    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl);
+    Log() << Verbose(0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
     return true;
   }
@@ -277 +252 @@
   double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
   if ((angle - M_PI) > -MYEPSILON) {
-    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
+    Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl;
     return true;
   } else {
-    DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
+    Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl;
     return false;
   }
@@ -289 +264 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
-{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 2; i++)
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+        Info FunctionInfo(__func__);
+  for(int i=0;i<2;i++)
     if (point == endpoints[i])
       return true;
   return false;
-}
-;
+};
 
 /** Returns other endpoint of the line.
@@ -303 +277 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
-{
-  Info FunctionInfo(__func__);
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
+{
+        Info FunctionInfo(__func__);
   if (endpoints[0] == point)
     return endpoints[1];
@@ -312 +286 @@
   else
     return NULL;
-}
-;
+};
 
 /** output operator for BoundaryLineSet.
@@ -319 +292 @@
  * \param &a boundary line
  */
-ostream & operator <<(ostream &ost, const BoundaryLineSet &a)
+ostream & operator <<(ostream &ost, const  BoundaryLineSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]";
   return ost;
-}
-;
+};
 
 // ======================================== Triangles on Boundary =================================
@@ -333 +305 @@
   Nr(-1)
 {
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++) {
-    endpoints[i] = NULL;
-    lines[i] = NULL;
-  }
-}
-;
+        Info FunctionInfo(__func__);
+  for (int i = 0; i < 3; i++)
+    {
+      endpoints[i] = NULL;
+      lines[i] = NULL;
+    }
+};
 
 /** Constructor for BoundaryTriangleSet with three lines.
@@ -345 +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)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   // set number
   // set lines
@@ -360 +332 @@
     // for all three lines
     for (int j = 0; j < 2; j++) { // for both endpoints
-      OrderMap.insert(pair<int, class BoundaryPointSet *> (line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
+      OrderMap.insert(pair<int, class BoundaryPointSet *> (
+          line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
       // and we don't care whether insertion fails
     }
   // set endpoints
   int Counter = 0;
-  DoLog(0) && (Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl);
+  Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl;
   for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) {
     endpoints[Counter] = runner->second;
-    DoLog(0) && (Log() << Verbose(0) << " " << *endpoints[Counter] << endl);
+    Log() << Verbose(0) << " " << *endpoints[Counter] << endl;
     Counter++;
   }
   if (Counter < 3) {
-    DoeLog(0) && (eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl);
+    eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl;
     performCriticalExit();
   }
-}
-;
+};
 
 /** Destructor of BoundaryTriangleSet.
@@ -384 +356 @@
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   for (int i = 0; i < 3; i++) {
     if (lines[i] != NULL) {
@@ -391 +363 @@
       }
       if (lines[i]->triangles.empty()) {
-        //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-        delete (lines[i]);
-        lines[i] = NULL;
+          //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+          delete (lines[i]);
+          lines[i] = NULL;
       }
     }
   }
   //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl;
-}
-;
+};
 
 /** Calculates the normal vector for this triangle.
@@ -405 +376 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
-{
-  Info FunctionInfo(__func__);
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+{
+        Info FunctionInfo(__func__);
   // get normal vector
   NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);
@@ -414 +385 @@
   if (NormalVector.ScalarProduct(&OtherVector) > 0.)
     NormalVector.Scale(-1.);
-  DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl);
-}
-;
-
-/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
+  Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl;
+};
+
+/** 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
@@ -430 +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;
 
   if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) {
-    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
+    eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;
     return false;
   }
 
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
-
-  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
+  // Calculate cross point between one baseline and the line from the third endpoint to intersection
+  int i=0;
+  do {
+    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
+      helper.CopyVector(endpoints[(i+1)%3]->node->node);
+      helper.SubtractVector(endpoints[i%3]->node->node);
+    } else 
+      i++;
+    if (i>2)
+      break;
+  } while (CrossPoint.NormSquared() < MYEPSILON);
+  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
     return true;
-  } else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
-    return true;
-  } else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {
-    DoLog(1) && (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;
-  do {
-    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i % 3]->node->node, endpoints[(i + 1) % 3]->node->node, endpoints[(i + 2) % 3]->node->node, Intersection, &NormalVector)) {
-      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();
-      DoLog(1) && (Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl);
-      if ((s < -MYEPSILON) || ((s - 1.) > MYEPSILON)) {
-        DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
-        i = 4;
-        break;
-      }
-      i++;
-    } else
-      break;
-  } while (i < 3);
-  if (i == 3) {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
-    return true;
-  } else {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl);
+  } else { // outside!
+    Intersection->Zero();
     return false;
   }
-}
-;
-
-/** Finds the point on the triangle to the point \a *x.
- * We call Vector::GetIntersectionWithPlane() with \a * and the center of the triangle to receive an intersection point.
- * Then we check the in-plane part (the part projected down onto plane). We check whether it crosses one of the
- * boundary lines. If it does, we return this intersection as closest point, otherwise the projected point down.
- * 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
-  DoLog(1) && (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);
-
-  DoLog(2) && (Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl);
-  DoLog(2) && (Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl);
-  DoLog(2) && (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();
-    DoLog(2) && (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
-      DoLog(2) && (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);
-      }
-    }
-  }
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl);
-  return ShortestDistance;
-}
-;
+};
 
 /** Checks whether lines is any of the three boundary lines this triangle contains.
@@ -575 +440 @@
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
-{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++)
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
+{
+        Info FunctionInfo(__func__);
+  for(int i=0;i<3;i++)
     if (line == lines[i])
       return true;
   return false;
-}
-;
+};
 
 /** Checks whether point is any of the three endpoints this triangle contains.
@@ -589 +453 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
-{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+        Info FunctionInfo(__func__);
+  for(int i=0;i<3;i++)
     if (point == endpoints[i])
       return true;
   return false;
-}
-;
+};
 
 /** Checks whether point is any of the three endpoints this triangle contains.
@@ -603 +466 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
-{
-  Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
+{
+        Info FunctionInfo(__func__);
+  for(int i=0;i<3;i++)
     if (point == endpoints[i]->node)
       return true;
   return false;
-}
-;
+};
 
 /** Checks whether three given \a *Points coincide with triangle's endpoints.
@@ -617 +479 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
-{
-  Info FunctionInfo(__func__);
-  DoLog(1) && (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]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2])
-
-  ));
-}
-;
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
+{
+        Info FunctionInfo(__func__);
+  return (((endpoints[0] == Points[0])
+            || (endpoints[0] == Points[1])
+            || (endpoints[0] == Points[2])
+          ) && (
+            (endpoints[1] == Points[0])
+            || (endpoints[1] == Points[1])
+            || (endpoints[1] == Points[2])
+          ) && (
+            (endpoints[2] == Points[0])
+            || (endpoints[2] == Points[1])
+            || (endpoints[2] == Points[2])
+
+          ));
+};
 
 /** Checks whether three given \a *Points coincide with triangle's endpoints.
@@ -631 +501 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
-{
-  Info FunctionInfo(__func__);
-  return (((endpoints[0] == T->endpoints[0]) || (endpoints[0] == T->endpoints[1]) || (endpoints[0] == T->endpoints[2])) && ((endpoints[1] == T->endpoints[0]) || (endpoints[1] == T->endpoints[1]) || (endpoints[1] == T->endpoints[2])) && ((endpoints[2] == T->endpoints[0]) || (endpoints[2] == T->endpoints[1]) || (endpoints[2] == T->endpoints[2])
-
-  ));
-}
-;
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
+{
+        Info FunctionInfo(__func__);
+  return (((endpoints[0] == T->endpoints[0])
+            || (endpoints[0] == T->endpoints[1])
+            || (endpoints[0] == T->endpoints[2])
+          ) && (
+            (endpoints[1] == T->endpoints[0])
+            || (endpoints[1] == T->endpoints[1])
+            || (endpoints[1] == T->endpoints[2])
+          ) && (
+            (endpoints[2] == T->endpoints[0])
+            || (endpoints[2] == T->endpoints[1])
+            || (endpoints[2] == T->endpoints[2])
+
+          ));
+};
 
 /** Returns the endpoint which is not contained in the given \a *line.
@@ -644 +523 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
-{
-  Info FunctionInfo(__func__);
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
+{
+        Info FunctionInfo(__func__);
   // sanity check
   if (!ContainsBoundaryLine(line))
     return NULL;
-  for (int i = 0; i < 3; i++)
+  for(int i=0;i<3;i++)
     if (!line->ContainsBoundaryPoint(endpoints[i]))
       return endpoints[i];
   // actually, that' impossible :)
   return NULL;
-}
-;
+};
 
 /** Calculates the center point of the triangle.
@@ -662 +540 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector * const center) const
-{
-  Info FunctionInfo(__func__);
+void BoundaryTriangleSet::GetCenter(Vector *center)
+{
+        Info FunctionInfo(__func__);
   center->Zero();
-  for (int i = 0; i < 3; i++)
+  for(int i=0;i<3;i++)
     center->AddVector(endpoints[i]->node->node);
-  center->Scale(1. / 3.);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl);
+  center->Scale(1./3.);
 }
 
@@ -679 +556 @@
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
-  //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
-  //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
+//  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
+//      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
   return ost;
-}
-;
-
-// ======================================== Polygons on Boundary =================================
-
-/** Constructor for BoundaryPolygonSet.
- */
-BoundaryPolygonSet::BoundaryPolygonSet() :
-  Nr(-1)
+};
+
+// =========================================================== class TESSELPOINT ===========================================
+
+/** Constructor of class TesselPoint.
+ */
+TesselPoint::TesselPoint()
 {
   Info FunctionInfo(__func__);
-}
-;
-
-/** Destructor of BoundaryPolygonSet.
- * Just clears endpoints.
- * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
- */
-BoundaryPolygonSet::~BoundaryPolygonSet()
-{
-  Info FunctionInfo(__func__);
-  endpoints.clear();
-  DoLog(1) && (Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl);
-}
-;
-
-/** Calculates the normal vector for this triangle.
- * Is made unique by comparison with \a OtherVector to point in the other direction.
- * \param &OtherVector direction vector to make normal vector unique.
- * \return allocated vector in normal direction
- */
-Vector * BoundaryPolygonSet::GetNormalVector(const Vector &OtherVector) const
-{
-  Info FunctionInfo(__func__);
-  // get normal vector
-  Vector TemporaryNormal;
-  Vector *TotalNormal = new Vector;
-  PointSet::const_iterator Runner[3];
-  for (int i = 0; i < 3; i++) {
-    Runner[i] = endpoints.begin();
-    for (int j = 0; j < i; j++) { // go as much further
-      Runner[i]++;
-      if (Runner[i] == endpoints.end()) {
-        DoeLog(0) && (eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl);
-        performCriticalExit();
-      }
-    }
-  }
-  TotalNormal->Zero();
-  int counter = 0;
-  for (; Runner[2] != endpoints.end();) {
-    TemporaryNormal.MakeNormalVector((*Runner[0])->node->node, (*Runner[1])->node->node, (*Runner[2])->node->node);
-    for (int i = 0; i < 3; i++) // increase each of them
-      Runner[i]++;
-    TotalNormal->AddVector(&TemporaryNormal);
-  }
-  TotalNormal->Scale(1. / (double) counter);
-
-  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (TotalNormal->ScalarProduct(&OtherVector) > 0.)
-    TotalNormal->Scale(-1.);
-  DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl);
-
-  return TotalNormal;
-}
-;
-
-/** Calculates the center point of the triangle.
- * Is third of the sum of all endpoints.
- * \param *center central point on return.
- */
-void BoundaryPolygonSet::GetCenter(Vector * const center) const
-{
-  Info FunctionInfo(__func__);
-  center->Zero();
-  int counter = 0;
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    center->AddVector((*Runner)->node->node);
-    counter++;
-  }
-  center->Scale(1. / (double) counter);
-  DoLog(1) && (Log() << Verbose(1) << "Center is at " << *center << "." << endl);
-}
-
-/** Checks whether the polygons contains all three endpoints of the triangle.
- * \param *triangle triangle to test
- * \return true - triangle is contained polygon, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryTriangle(const BoundaryTriangleSet * const triangle) const
-{
-  Info FunctionInfo(__func__);
-  return ContainsPresentTupel(triangle->endpoints, 3);
-}
-;
-
-/** Checks whether the polygons contains both endpoints of the line.
- * \param *line line to test
- * \return true - line is of the triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
-{
-  Info FunctionInfo(__func__);
-  return ContainsPresentTupel(line->endpoints, 2);
-}
-;
-
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point point to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
-{
-  Info FunctionInfo(__func__);
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    DoLog(0) && (Log() << Verbose(0) << "Checking against " << **Runner << endl);
-    if (point == (*Runner)) {
-      DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
-      return true;
-    }
-  }
-  DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
-  return false;
-}
-;
-
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point TesselPoint to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsBoundaryPoint(const TesselPoint * const point) const
-{
-  Info FunctionInfo(__func__);
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
-    if (point == (*Runner)->node) {
-      DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
-      return true;
-    }
-  DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
-  return false;
-}
-;
-
-/** Checks whether given array of \a *Points coincide with polygons's endpoints.
- * \param **Points pointer to an array of BoundaryPointSet
- * \param dim dimension of array
- * \return true - set of points is contained in polygon, false - is not
- */
-bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const
-{
-  Info FunctionInfo(__func__);
-  int counter = 0;
-  DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
-  for (int i = 0; i < dim; i++) {
-    DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl);
-    if (ContainsBoundaryPoint(Points[i])) {
-      counter++;
-    }
-  }
-
-  if (counter == dim)
-    return true;
-  else
-    return false;
-}
-;
-
-/** Checks whether given PointList coincide with polygons's endpoints.
- * \param &endpoints PointList
- * \return true - set of points is contained in polygon, false - is not
- */
-bool BoundaryPolygonSet::ContainsPresentTupel(const PointSet &endpoints) const
-{
-  Info FunctionInfo(__func__);
-  size_t counter = 0;
-  DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << **Runner << endl);
-    if (ContainsBoundaryPoint(*Runner))
-      counter++;
-  }
-
-  if (counter == endpoints.size())
-    return true;
-  else
-    return false;
-}
-;
-
-/** Checks whether given set of \a *Points coincide with polygons's endpoints.
- * \param *P pointer to BoundaryPolygonSet
- * \return true - is the very triangle, false - is not
- */
-bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const
-{
-  return ContainsPresentTupel((const PointSet) P->endpoints);
-}
-;
-
-/** Gathers all the endpoints' triangles in a unique set.
- * \return set of all triangles
- */
-TriangleSet * BoundaryPolygonSet::GetAllContainedTrianglesFromEndpoints() const
-{
-  Info FunctionInfo(__func__);
-  pair<TriangleSet::iterator, bool> Tester;
-  TriangleSet *triangles = new TriangleSet;
-
-  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
-    for (LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)
-      for (TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) {
-        //Log() << Verbose(0) << " Testing triangle " << *(Sprinter->second) << endl;
-        if (ContainsBoundaryTriangle(Sprinter->second)) {
-          Tester = triangles->insert(Sprinter->second);
-          if (Tester.second)
-            DoLog(0) && (Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl);
-        }
-      }
-
-  DoLog(1) && (Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl);
-  return triangles;
-}
-;
-
-/** Fills the endpoints of this polygon from the triangles attached to \a *line.
- * \param *line lines with triangles attached
- * \return true - polygon contains endpoints, false - line was NULL
- */
-bool BoundaryPolygonSet::FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line)
-{
-  Info FunctionInfo(__func__);
-  pair<PointSet::iterator, bool> Tester;
-  if (line == NULL)
-    return false;
-  DoLog(1) && (Log() << Verbose(1) << "Filling polygon from line " << *line << endl);
-  for (TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {
-    for (int i = 0; i < 3; i++) {
-      Tester = endpoints.insert((Runner->second)->endpoints[i]);
-      if (Tester.second)
-        DoLog(1) && (Log() << Verbose(1) << "  Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl);
-    }
-  }
-
-  return true;
-}
-;
-
-/** output operator for BoundaryPolygonSet.
- * \param &ost output stream
- * \param &a boundary polygon
- */
-ostream &operator <<(ostream &ost, const BoundaryPolygonSet &a)
-{
-  ost << "[" << a.Nr << "|";
-  for (PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) {
-    ost << (*Runner)->node->Name;
-    Runner++;
-    if (Runner != a.endpoints.end())
-      ost << ",";
-  }
-  ost << "]";
-  return ost;
-}
-;
-
-// =========================================================== class TESSELPOINT ===========================================
-
-/** Constructor of class TesselPoint.
- */
-TesselPoint::TesselPoint()
-{
-  //Info FunctionInfo(__func__);
   node = NULL;
   nr = -1;
-  Name = NULL;
-}
-;
+  Name =  NULL;
+};
 
 /** Destructor for class TesselPoint.
@@ -963 +577 @@
 TesselPoint::~TesselPoint()
 {
-  //Info FunctionInfo(__func__);
-}
-;
+  Info FunctionInfo(__func__);
+};
 
 /** Prints LCNode to screen.
  */
-ostream & operator <<(ostream &ost, const TesselPoint &a)
+ostream & operator << (ostream &ost, const TesselPoint &a)
 {
   ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
   return ost;
-}
-;
+};
 
 /** Prints LCNode to screen.
  */
-ostream & TesselPoint::operator <<(ostream &ost)
-{
-  Info FunctionInfo(__func__);
+ostream & TesselPoint::operator << (ostream &ost)
+{
+        Info FunctionInfo(__func__);
   ost << "[" << (nr) << "|" << this << "]";
   return ost;
-}
-;
+};
+
 
 // =========================================================== class POINTCLOUD ============================================
@@ -992 +604 @@
 PointCloud::PointCloud()
 {
-  //Info FunctionInfo(__func__);
-}
-;
+        Info FunctionInfo(__func__);
+};
 
 /** Destructor for class PointCloud.
@@ -1000 +611 @@
 PointCloud::~PointCloud()
 {
-  //Info FunctionInfo(__func__);
-}
-;
+        Info FunctionInfo(__func__);
+};
 
 // ============================ CandidateForTesselation =============================
@@ -1008 +618 @@
 /** Constructor of class CandidateForTesselation.
  */
-CandidateForTesselation::CandidateForTesselation(BoundaryLineSet* line) :
-  BaseLine(line), ThirdPoint(NULL), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
-{
-  Info FunctionInfo(__func__);
-}
-;
+CandidateForTesselation::CandidateForTesselation (BoundaryLineSet* line) :
+  BaseLine(line),
+  ShortestAngle(2.*M_PI),
+  OtherShortestAngle(2.*M_PI)
+{
+        Info FunctionInfo(__func__);
+};
+
 
 /** Constructor of class CandidateForTesselation.
  */
-CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
-  BaseLine(line), ThirdPoint(point), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
-{
-  Info FunctionInfo(__func__);
+CandidateForTesselation::CandidateForTesselation (TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
+    BaseLine(line),
+    ShortestAngle(2.*M_PI),
+    OtherShortestAngle(2.*M_PI)
+{
+        Info FunctionInfo(__func__);
   OptCenter.CopyVector(&OptCandidateCenter);
   OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
-}
-;
+};
 
 /** Destructor for class CandidateForTesselation.
  */
-CandidateForTesselation::~CandidateForTesselation()
-{
-}
-;
-
-/** Checks validity of a given sphere of a candidate line.
- * Sphere must touch all candidates and the baseline endpoints and there must be no other atoms inside.
- * \param RADIUS radius of sphere
- * \param *LC LinkedCell structure with other atoms
- * \return true - sphere is valid, false - sphere contains other points
- */
-bool CandidateForTesselation::CheckValidity(const double RADIUS, const LinkedCell *LC) const
-{
-  Info FunctionInfo(__func__);
-
-  const double radiusSquared = RADIUS * RADIUS;
-  list<const Vector *> VectorList;
-  VectorList.push_back(&OptCenter);
-  //VectorList.push_back(&OtherOptCenter);  // don't check the other (wrong) center
-
-  if (!pointlist.empty())
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
-  else
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
-  // check baseline for OptCenter and OtherOptCenter being on sphere's surface
-  for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-    for (int i = 0; i < 2; i++) {
-      const double distance = fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->node) - radiusSquared);
-      if (distance > HULLEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
-        return false;
-      }
-    }
-  }
-
-  // check Candidates for OptCenter and OtherOptCenter being on sphere's surface
-  for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
-    const TesselPoint *Walker = *Runner;
-    for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-      const double distance = fabs((*VRunner)->DistanceSquared(Walker->node) - radiusSquared);
-      if (distance > HULLEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
-        return false;
-      } else {
-        DoLog(1) && (Log() << Verbose(1) << "Candidate " << *Walker << " is inside by " << distance << "." << endl);
-      }
-    }
-  }
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
-  bool flag = true;
-  for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-    // get all points inside the sphere
-    TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));
-
-    DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << OtherOptCenter << ":" << endl);
-    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&OtherOptCenter) << "." << endl);
-
-    // remove baseline's endpoints and candidates
-    for (int i = 0; i < 2; i++) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "." << endl);
-      ListofPoints->remove(BaseLine->endpoints[i]->node);
-    }
-    for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: removing candidate tesselpoint " << *(*Runner) << "." << endl);
-      ListofPoints->remove(*Runner);
-    }
-    if (!ListofPoints->empty()) {
-      DoeLog(1) && (eLog() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
-      flag = false;
-      DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl);
-      for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-        DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << endl);
-    }
-    delete (ListofPoints);
-
-    // check with animate_sphere.tcl VMD script
-    if (ThirdPoint != NULL) {
-      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl);
-    } else {
-      DoLog(1) && (Log() << Verbose(1) << "Check by: ... missing third point ..." << endl);
-      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl);
-    }
-  }
-  return flag;
-}
-;
+CandidateForTesselation::~CandidateForTesselation() {
+  BaseLine = NULL;
+};
 
 /** output operator for CandidateForTesselation.
@@ -1121 +649 @@
  * \param &a boundary line
  */
-ostream & operator <<(ostream &ost, const CandidateForTesselation &a)
+ostream & operator <<(ostream &ost, const  CandidateForTesselation &a)
 {
   ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->Name << "," << a.BaseLine->endpoints[1]->node->Name << "] with ";
@@ -1134 +662 @@
     for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++)
       ost << *(*Runner) << " ";
-    ost << " at angle " << (a.ShortestAngle) << ".";
+    ost << " at angle " << (a.ShortestAngle)<< ".";
   }
 
   return ost;
+};
+
+
+// =========================================================== class TESSELATION ===========================================
+
+/** Constructor of class Tesselation.
+ */
+Tesselation::Tesselation() :
+  PointsOnBoundaryCount(0),
+  LinesOnBoundaryCount(0),
+  TrianglesOnBoundaryCount(0),
+  LastTriangle(NULL),
+  TriangleFilesWritten(0),
+  InternalPointer(PointsOnBoundary.begin())
+{
+        Info FunctionInfo(__func__);
 }
 ;
 
-// =========================================================== class TESSELATION ===========================================
-
-/** Constructor of class Tesselation.
- */
-Tesselation::Tesselation() :
-  PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin())
-{
-  Info FunctionInfo(__func__);
-}
-;
-
 /** Destructor of class Tesselation.
  * We have to free all points, lines and triangles.
@@ -1157 +690 @@
 Tesselation::~Tesselation()
 {
-  Info FunctionInfo(__func__);
-  DoLog(0) && (Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl);
+        Info FunctionInfo(__func__);
+  Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl;
   for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
     if (runner->second != NULL) {
@@ -1164 +697 @@
       runner->second = NULL;
     } else
-      DoeLog(1) && (eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl);
-  }
-  DoLog(0) && (Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl);
+      eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl;
+  }
+  Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
 }
 ;
@@ -1172 +705 @@
 /** PointCloud implementation of GetCenter
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 Vector * Tesselation::GetCenter(ofstream *out) const
 {
-  Info FunctionInfo(__func__);
-  Vector *Center = new Vector(0., 0., 0.);
-  int num = 0;
+        Info FunctionInfo(__func__);
+  Vector *Center = new Vector(0.,0.,0.);
+  int num=0;
   for (GoToFirst(); (!IsEnd()); GoToNext()) {
     Center->AddVector(GetPoint()->node);
     num++;
   }
-  Center->Scale(1. / num);
+  Center->Scale(1./num);
   return Center;
-}
-;
+};
 
 /** PointCloud implementation of GoPoint
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 TesselPoint * Tesselation::GetPoint() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   return (InternalPointer->second->node);
-}
-;
+};
 
 /** PointCloud implementation of GetTerminalPoint.
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 TesselPoint * Tesselation::GetTerminalPoint() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   PointMap::const_iterator Runner = PointsOnBoundary.end();
   Runner--;
   return (Runner->second->node);
-}
-;
+};
 
 /** PointCloud implementation of GoToNext.
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 void Tesselation::GoToNext() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   if (InternalPointer != PointsOnBoundary.end())
     InternalPointer++;
-}
-;
+};
 
 /** PointCloud implementation of GoToPrevious.
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 void Tesselation::GoToPrevious() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   if (InternalPointer != PointsOnBoundary.begin())
     InternalPointer--;
-}
-;
+};
 
 /** PointCloud implementation of GoToFirst.
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 void Tesselation::GoToFirst() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   InternalPointer = PointsOnBoundary.begin();
-}
-;
+};
 
 /** PointCloud implementation of GoToLast.
@@ -1246 +773 @@
 void Tesselation::GoToLast() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   InternalPointer = PointsOnBoundary.end();
   InternalPointer--;
-}
-;
+};
 
 /** PointCloud implementation of IsEmpty.
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 bool Tesselation::IsEmpty() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   return (PointsOnBoundary.empty());
-}
-;
+};
 
 /** PointCloud implementation of IsLast.
  * Uses PointsOnBoundary and STL stuff.
- */
+ */   
 bool Tesselation::IsEnd() const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   return (InternalPointer == PointsOnBoundary.end());
-}
-;
+};
+
 
 /** Gueses first starting triangle of the convex envelope.
@@ -1277 +802 @@
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
  */
-void Tesselation::GuessStartingTriangle()
-{
-  Info FunctionInfo(__func__);
+void
+Tesselation::GuessStartingTriangle()
+{
+        Info FunctionInfo(__func__);
   // 4b. create a starting triangle
   // 4b1. create all distances
@@ -1289 +815 @@
 
   // with A chosen, take each pair B,C and sort
-  if (A != PointsOnBoundary.end()) {
-    B = A;
-    B++;
-    for (; B != PointsOnBoundary.end(); B++) {
-      C = B;
-      C++;
-      for (; C != PointsOnBoundary.end(); C++) {
-        tmp = A->second->node->node->DistanceSquared(B->second->node->node);
-        distance = tmp * tmp;
-        tmp = A->second->node->node->DistanceSquared(C->second->node->node);
-        distance += tmp * tmp;
-        tmp = B->second->node->node->DistanceSquared(C->second->node->node);
-        distance += tmp * tmp;
-        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
-      }
-    }
-  }
+  if (A != PointsOnBoundary.end())
+    {
+      B = A;
+      B++;
+      for (; B != PointsOnBoundary.end(); B++)
+        {
+          C = B;
+          C++;
+          for (; C != PointsOnBoundary.end(); C++)
+            {
+              tmp = A->second->node->node->DistanceSquared(B->second->node->node);
+              distance = tmp * tmp;
+              tmp = A->second->node->node->DistanceSquared(C->second->node->node);
+              distance += tmp * tmp;
+              tmp = B->second->node->node->DistanceSquared(C->second->node->node);
+              distance += tmp * tmp;
+              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
+            }
+        }
+    }
   //    // listing distances
   //    Log() << Verbose(1) << "Listing DistanceMMap:";
@@ -1315 +844 @@
   // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
   DistanceMultiMap::iterator baseline = DistanceMMap.begin();
-  for (; baseline != DistanceMMap.end(); baseline++) {
-    // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
-    // 2. next, we have to check whether all points reside on only one side of the triangle
-    // 3. construct plane vector
-    PlaneVector.MakeNormalVector(A->second->node->node, baseline->second.first->second->node->node, baseline->second.second->second->node->node);
-    DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl);
-    // 4. loop over all points
-    double sign = 0.;
-    PointMap::iterator checker = PointsOnBoundary.begin();
-    for (; checker != PointsOnBoundary.end(); checker++) {
-      // (neglecting A,B,C)
-      if ((checker == A) || (checker == baseline->second.first) || (checker == baseline->second.second))
-        continue;
-      // 4a. project onto plane vector
-      TrialVector.CopyVector(checker->second->node->node);
-      TrialVector.SubtractVector(A->second->node->node);
-      distance = TrialVector.ScalarProduct(&PlaneVector);
-      if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
-        continue;
-      DoLog(2) && (Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl);
-      tmp = distance / fabs(distance);
-      // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
-      if ((sign != 0) && (tmp != sign)) {
-        // 4c. If so, break 4. loop and continue with next candidate in 1. loop
-        DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name << "," << baseline->second.second->second->node->Name << " leaves " << checker->second->node->Name << " outside the convex hull." << endl);
-        break;
-      } else { // note the sign for later
-        DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name << "," << baseline->second.second->second->node->Name << " leave " << checker->second->node->Name << " inside the convex hull." << endl);
-        sign = tmp;
-      }
-      // 4d. Check whether the point is inside the triangle (check distance to each node
-      tmp = checker->second->node->node->DistanceSquared(A->second->node->node);
-      int innerpoint = 0;
-      if ((tmp < A->second->node->node->DistanceSquared(baseline->second.first->second->node->node)) && (tmp < A->second->node->node->DistanceSquared(baseline->second.second->second->node->node)))
-        innerpoint++;
-      tmp = checker->second->node->node->DistanceSquared(baseline->second.first->second->node->node);
-      if ((tmp < baseline->second.first->second->node->node->DistanceSquared(A->second->node->node)) && (tmp < baseline->second.first->second->node->node->DistanceSquared(baseline->second.second->second->node->node)))
-        innerpoint++;
-      tmp = checker->second->node->node->DistanceSquared(baseline->second.second->second->node->node);
-      if ((tmp < baseline->second.second->second->node->node->DistanceSquared(baseline->second.first->second->node->node)) && (tmp < baseline->second.second->second->node->node->DistanceSquared(A->second->node->node)))
-        innerpoint++;
-      // 4e. If so, break 4. loop and continue with next candidate in 1. loop
-      if (innerpoint == 3)
-        break;
-    }
-    // 5. come this far, all on same side? Then break 1. loop and construct triangle
-    if (checker == PointsOnBoundary.end()) {
-      DoLog(2) && (Log() << Verbose(2) << "Looks like we have a candidate!" << endl);
-      break;
-    }
-  }
-  if (baseline != DistanceMMap.end()) {
-    BPS[0] = baseline->second.first->second;
-    BPS[1] = baseline->second.second->second;
-    BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-    BPS[0] = A->second;
-    BPS[1] = baseline->second.second->second;
-    BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-    BPS[0] = baseline->second.first->second;
-    BPS[1] = A->second;
-    BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-
-    // 4b3. insert created triangle
-    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-    TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-    TrianglesOnBoundaryCount++;
-    for (int i = 0; i < NDIM; i++) {
-      LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
-      LinesOnBoundaryCount++;
-    }
-
-    DoLog(1) && (Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl);
-  } else {
-    DoeLog(0) && (eLog() << Verbose(0) << "No starting triangle found." << endl);
-  }
+  for (; baseline != DistanceMMap.end(); baseline++)
+    {
+      // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+      // 2. next, we have to check whether all points reside on only one side of the triangle
+      // 3. construct plane vector
+      PlaneVector.MakeNormalVector(A->second->node->node,
+          baseline->second.first->second->node->node,
+          baseline->second.second->second->node->node);
+      Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl;
+      // 4. loop over all points
+      double sign = 0.;
+      PointMap::iterator checker = PointsOnBoundary.begin();
+      for (; checker != PointsOnBoundary.end(); checker++)
+        {
+          // (neglecting A,B,C)
+          if ((checker == A) || (checker == baseline->second.first) || (checker
+              == baseline->second.second))
+            continue;
+          // 4a. project onto plane vector
+          TrialVector.CopyVector(checker->second->node->node);
+          TrialVector.SubtractVector(A->second->node->node);
+          distance = TrialVector.ScalarProduct(&PlaneVector);
+          if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
+            continue;
+          Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl;
+          tmp = distance / fabs(distance);
+          // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
+          if ((sign != 0) && (tmp != sign))
+            {
+              // 4c. If so, break 4. loop and continue with next candidate in 1. loop
+              Log() << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leaves "
+                  << checker->second->node->Name << " outside the convex hull."
+                  << endl;
+              break;
+            }
+          else
+            { // note the sign for later
+              Log() << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leave "
+                  << checker->second->node->Name << " inside the convex hull."
+                  << endl;
+              sign = tmp;
+            }
+          // 4d. Check whether the point is inside the triangle (check distance to each node
+          tmp = checker->second->node->node->DistanceSquared(A->second->node->node);
+          int innerpoint = 0;
+          if ((tmp < A->second->node->node->DistanceSquared(
+              baseline->second.first->second->node->node)) && (tmp
+              < A->second->node->node->DistanceSquared(
+                  baseline->second.second->second->node->node)))
+            innerpoint++;
+          tmp = checker->second->node->node->DistanceSquared(
+              baseline->second.first->second->node->node);
+          if ((tmp < baseline->second.first->second->node->node->DistanceSquared(
+              A->second->node->node)) && (tmp
+              < baseline->second.first->second->node->node->DistanceSquared(
+                  baseline->second.second->second->node->node)))
+            innerpoint++;
+          tmp = checker->second->node->node->DistanceSquared(
+              baseline->second.second->second->node->node);
+          if ((tmp < baseline->second.second->second->node->node->DistanceSquared(
+              baseline->second.first->second->node->node)) && (tmp
+              < baseline->second.second->second->node->node->DistanceSquared(
+                  A->second->node->node)))
+            innerpoint++;
+          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+          if (innerpoint == 3)
+            break;
+        }
+      // 5. come this far, all on same side? Then break 1. loop and construct triangle
+      if (checker == PointsOnBoundary.end())
+        {
+          Log() << Verbose(2) << "Looks like we have a candidate!" << endl;
+          break;
+        }
+    }
+  if (baseline != DistanceMMap.end())
+    {
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = A->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = A->second;
+      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+
+      // 4b3. insert created triangle
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+      TrianglesOnBoundaryCount++;
+      for (int i = 0; i < NDIM; i++)
+        {
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
+          LinesOnBoundaryCount++;
+        }
+
+      Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+    }
+  else
+    {
+      eLog() << Verbose(0) << "No starting triangle found." << endl;
+    }
 }
 ;
@@ -1408 +971 @@
 void Tesselation::TesselateOnBoundary(const PointCloud * const cloud)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   bool flag;
   PointMap::iterator winner;
@@ -1427 +990 @@
         // get peak point with respect to this base line's only triangle
         BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
-        DoLog(0) && (Log() << Verbose(0) << "Current baseline is between " << *(baseline->second) << "." << endl);
+        Log() << Verbose(0) << "Current baseline is between " << *(baseline->second) << "." << endl;
         for (int i = 0; i < 3; i++)
           if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
             peak = BTS->endpoints[i];
-        DoLog(1) && (Log() << Verbose(1) << " and has peak " << *peak << "." << endl);
+        Log() << Verbose(1) << " and has peak " << *peak << "." << endl;
 
         // prepare some auxiliary vectors
@@ -1446 +1009 @@
           CenterVector.AddVector(BTS->endpoints[i]->node->node);
         CenterVector.Scale(1. / 3.);
-        DoLog(2) && (Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl);
+        Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl;
 
         // normal vector of triangle
@@ -1453 +1016 @@
         BTS->GetNormalVector(NormalVector);
         NormalVector.CopyVector(&BTS->NormalVector);
-        DoLog(2) && (Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl);
+        Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl;
 
         // vector in propagation direction (out of triangle)
@@ -1463 +1026 @@
         if (PropagationVector.ScalarProduct(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
           PropagationVector.Scale(-1.);
-        DoLog(2) && (Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl);
+        Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl;
         winner = PointsOnBoundary.end();
 
@@ -1469 +1032 @@
         for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
           if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
-            DoLog(1) && (Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl);
+            Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl;
 
             // first check direction, so that triangles don't intersect
@@ -1476 +1039 @@
             VirtualNormalVector.ProjectOntoPlane(&NormalVector);
             TempAngle = VirtualNormalVector.Angle(&PropagationVector);
-            DoLog(2) && (Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl);
-            if (TempAngle > (M_PI / 2.)) { // no bends bigger than Pi/2 (90 degrees)
-              DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl);
+            Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
+            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
+              Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
               continue;
             } else
-              DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl);
+              Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
 
             // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
@@ -1487 +1050 @@
             LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
             if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) {
-              DoLog(2) && (Log() << Verbose(2) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl);
+              Log() << Verbose(2) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl;
               continue;
             }
             if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) {
-              DoLog(2) && (Log() << Verbose(2) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl);
+              Log() << Verbose(2) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl;
               continue;
             }
@@ -1497 +1060 @@
             // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
             if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
-              DoLog(4) && (Log() << Verbose(4) << "Current target is peak!" << endl);
+              Log() << Verbose(4) << "Current target is peak!" << endl;
               continue;
             }
@@ -1508 +1071 @@
             helper.ProjectOntoPlane(&TempVector);
             if (fabs(helper.NormSquared()) < MYEPSILON) {
-              DoLog(2) && (Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl);
+              Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl;
               continue;
             }
@@ -1518 +1081 @@
             TempVector.AddVector(baseline->second->endpoints[1]->node->node);
             TempVector.AddVector(target->second->node->node);
-            TempVector.Scale(1. / 3.);
+            TempVector.Scale(1./3.);
             TempVector.SubtractVector(Center);
             // make it always point outward
@@ -1525 +1088 @@
             // calculate angle
             TempAngle = NormalVector.Angle(&VirtualNormalVector);
-            DoLog(2) && (Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl);
+            Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
             if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
               SmallestAngle = TempAngle;
               winner = target;
-              DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);
+              Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
             } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
               // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
@@ -1547 +1110 @@
                 SmallestAngle = TempAngle;
                 winner = target;
-                DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl);
+                Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
               } else
-                DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl);
+                Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
             } else
-              DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);
+              Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
           }
         } // end of loop over all boundary points
@@ -1557 +1120 @@
         // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
         if (winner != PointsOnBoundary.end()) {
-          DoLog(0) && (Log() << Verbose(0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl);
+          Log() << Verbose(0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;
           // create the lins of not yet present
           BLS[0] = baseline->second;
@@ -1587 +1150 @@
           TrianglesOnBoundaryCount++;
         } else {
-          DoeLog(2) && (eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl);
+          eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
         }
 
         // 5d. If the set of lines is not yet empty, go to 5. and continue
       } else
-        DoLog(0) && (Log() << Verbose(0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl);
+        Log() << Verbose(0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl;
   } while (flag);
 
   // exit
-  delete (Center);
-}
-;
+  delete(Center);
+};
 
 /** Inserts all points outside of the tesselated surface into it by adding new triangles.
@@ -1608 +1170 @@
 bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   Vector Intersection, Normal;
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  TriangleList *triangles = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
@@ -1618 +1180 @@
   cloud->GoToFirst();
   BoundaryPoints = new LinkedCell(this, 5.);
-  while (!cloud->IsEnd()) { // we only have to go once through all points, as boundary can become only bigger
+  while (!cloud->IsEnd()) {  // we only have to go once through all points, as boundary can become only bigger
     if (AddFlag) {
-      delete (BoundaryPoints);
+      delete(BoundaryPoints);
       BoundaryPoints = new LinkedCell(this, 5.);
       AddFlag = false;
     }
     Walker = cloud->GetPoint();
-    DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);
+    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))) {
-      DoLog(0) && (Log() << Verbose(0) << "No triangles found, probably a tesselation point itself." << endl);
+      Log() << Verbose(0) << "No triangles found, probably a tesselation point itself." << endl;
       cloud->GoToNext();
       continue;
     } else {
     }
-    DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);
+    Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl;
     // get the intersection point
     if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) {
-      DoLog(0) && (Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl);
+      Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl;
       // we have the intersection, check whether in- or outside of boundary
       if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) {
         // inside, next!
-        DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
+        Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;
       } else {
         // outside!
-        DoLog(0) && (Log() << Verbose(0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl);
+        Log() << Verbose(0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl;
         class BoundaryLineSet *OldLines[3], *NewLines[3];
         class BoundaryPointSet *OldPoints[3], *NewPoint;
         // store the three old lines and old points
-        for (int i = 0; i < 3; i++) {
+        for (int i=0;i<3;i++) {
           OldLines[i] = BTS->lines[i];
           OldPoints[i] = BTS->endpoints[i];
@@ -1655 +1217 @@
         Normal.CopyVector(&BTS->NormalVector);
         // add Walker to boundary points
-        DoLog(0) && (Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl);
+        Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl;
         AddFlag = true;
-        if (AddBoundaryPoint(Walker, 0))
+        if (AddBoundaryPoint(Walker,0))
           NewPoint = BPS[0];
         else
           continue;
         // remove triangle
-        DoLog(0) && (Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl);
+        Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl;
         TrianglesOnBoundary.erase(BTS->Nr);
-        delete (BTS);
+        delete(BTS);
         // create three new boundary lines
-        for (int i = 0; i < 3; i++) {
+        for (int i=0;i<3;i++) {
           BPS[0] = NewPoint;
           BPS[1] = OldPoints[i];
           NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-          DoLog(1) && (Log() << Verbose(1) << "Creating new line " << *NewLines[i] << "." << endl);
+          Log() << Verbose(1) << "Creating new line " << *NewLines[i] << "." << endl;
           LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
           LinesOnBoundaryCount++;
         }
         // create three new triangle with new point
-        for (int i = 0; i < 3; i++) { // find all baselines
+        for (int i=0;i<3;i++) { // find all baselines
           BLS[0] = OldLines[i];
           int n = 1;
-          for (int j = 0; j < 3; j++) {
+          for (int j=0;j<3;j++) {
             if (NewLines[j]->IsConnectedTo(BLS[0])) {
-              if (n > 2) {
-                DoeLog(2) && (eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl);
+              if (n>2) {
+                eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl;
                 return false;
               } else
@@ -1692 +1254 @@
           BTS->GetNormalVector(Normal);
           Normal.Scale(-1.);
-          DoLog(0) && (Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl);
+          Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl;
           TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
           TrianglesOnBoundaryCount++;
@@ -1698 +1260 @@
       }
     } else { // something is wrong with FindClosestTriangleToPoint!
-      DoeLog(1) && (eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl);
+      eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl;
       return false;
     }
@@ -1705 +1267 @@
 
   // exit
-  delete (Center);
+  delete(Center);
   return true;
-}
-;
+};
 
 /** Adds a point to the tesselation::PointsOnBoundary list.
@@ -1717 +1278 @@
 bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   PointTestPair InsertUnique;
   BPS[n] = new class BoundaryPointSet(Walker);
@@ -1725 +1286 @@
     return true;
   } else {
-    delete (BPS[n]);
+    delete(BPS[n]);
     BPS[n] = InsertUnique.first->second;
     return false;
@@ -1739 +1300 @@
 void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   PointTestPair InsertUnique;
   TPS[n] = new class BoundaryPointSet(Candidate);
@@ -1747 +1308 @@
   } else {
     delete TPS[n];
-    DoLog(0) && (Log() << Verbose(0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl);
+    Log() << Verbose(0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
     TPS[n] = (InsertUnique.first)->second;
   }
@@ -1760 +1321 @@
 void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr);
   if (FindPoint != PointsOnBoundary.end())
@@ -1766 +1327 @@
   else
     TPS[n] = NULL;
-}
-;
+};
 
 /** Function tries to add line from current Points in BPS to BoundaryLineSet.
  * If successful it raises the line count and inserts the new line into the BLS,
  * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
- * @param *OptCenter desired OptCenter if there are more than one candidate line
- * @param *candidate third point of the triangle to be, for checking between multiple open line candidates
  * @param *a first endpoint
  * @param *b second endpoint
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
-{
+void Tesselation::AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) {
   bool insertNewLine = true;
-  LineMap::iterator FindLine = a->lines.find(b->node->nr);
-  BoundaryLineSet *WinningLine = NULL;
-  if (FindLine != a->lines.end()) {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl);
-
-    pair<LineMap::iterator, LineMap::iterator> FindPair;
+
+  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);
-
-    for (FindLine = FindPair.first; (FindLine != FindPair.second) && (insertNewLine); FindLine++) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl);
+    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++) {
       // If there is a line with less than two attached triangles, we don't need a new line.
-      if (FindLine->second->triangles.size() == 1) {
-        CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
-        if (!Finder->second->pointlist.empty())
-          DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl);
-        else
-          DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate." << endl);
-        // get open line
-        for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
-          if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(&Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches
-            DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << "." << endl);
-            insertNewLine = false;
-            WinningLine = FindLine->second;
-            break;
-          } else {
-            DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << "." << endl);
-          }
-        }
+      if (FindLine->second->triangles.size() < 2) {
+        insertNewLine = false;
+        Log() << Verbose(0) << "Using existing line " << *FindLine->second << endl;
+
+        BPS[0] = FindLine->second->endpoints[0];
+        BPS[1] = FindLine->second->endpoints[1];
+        BLS[n] = FindLine->second;
+
+        // remove existing line from OpenLines
+        CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
+        delete(CandidateLine->second);
+        OpenLines.erase(CandidateLine);
+
+        break;
       }
     }
@@ -1814 +1366 @@
 
   if (insertNewLine) {
-    AddNewTesselationTriangleLine(a, b, n);
-  } else {
-    AddExistingTesselationTriangleLine(WinningLine, n);
+    AlwaysAddTesselationTriangleLine(a, b, n);
   }
 }
@@ -1829 +1379 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
-{
-  Info FunctionInfo(__func__);
-  DoLog(0) && (Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl);
+void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
+{
+        Info FunctionInfo(__func__);
+  Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
-  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); // this also adds the line to the local maps
+  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
   // add line to global map
   LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
@@ -1842 +1392 @@
   // also add to open lines
   CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
-  OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
-}
-;
-
-/** Uses an existing line for a new triangle.
- * Sets Tesselation::BLS[\a n] and removes the lines from Tesselation::OpenLines.
- * \param *FindLine the line to add
- * \param n index of the line to set in Tesselation::BLS
- */
-void Tesselation::AddExistingTesselationTriangleLine(class BoundaryLineSet *Line, int n)
-{
-  Info FunctionInfo(__func__);
-  DoLog(0) && (Log() << Verbose(0) << "Using existing line " << *Line << endl);
-
-  // set endpoints and line
-  BPS[0] = Line->endpoints[0];
-  BPS[1] = Line->endpoints[1];
-  BLS[n] = Line;
-  // remove existing line from OpenLines
-  CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
-  if (CandidateLine != OpenLines.end()) {
-    DoLog(1) && (Log() << Verbose(1) << " Removing line from OpenLines." << endl);
-    delete (CandidateLine->second);
-    OpenLines.erase(CandidateLine);
-  } else {
-    DoeLog(1) && (eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
-  }
-}
-;
+  OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
+};
 
 /** Function adds triangle to global list.
@@ -1877 +1400 @@
 void Tesselation::AddTesselationTriangle()
 {
-  Info FunctionInfo(__func__);
-  DoLog(1) && (Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl);
+        Info FunctionInfo(__func__);
+  Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
 
   // add triangle to global map
@@ -1888 +1411 @@
 
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
-}
-;
+};
 
 /** Function adds triangle to global list.
@@ -1897 +1419 @@
 void Tesselation::AddTesselationTriangle(const int nr)
 {
-  Info FunctionInfo(__func__);
-  DoLog(0) && (Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl);
+        Info FunctionInfo(__func__);
+  Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl;
 
   // add triangle to global map
@@ -1907 +1429 @@
 
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
-}
-;
+};
 
 /** Removes a triangle from the tesselation.
@@ -1917 +1438 @@
 void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   if (triangle == NULL)
     return;
   for (int i = 0; i < 3; i++) {
     if (triangle->lines[i] != NULL) {
-      DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl);
+      Log() << Verbose(0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl;
       triangle->lines[i]->triangles.erase(triangle->Nr);
       if (triangle->lines[i]->triangles.empty()) {
-        DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl);
-        RemoveTesselationLine(triangle->lines[i]);
+          Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
+          RemoveTesselationLine(triangle->lines[i]);
       } else {
-        DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: ");
-        OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL));
-        for (TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
-          DoLog(0) && (Log() << Verbose(0) << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t");
-        DoLog(0) && (Log() << Verbose(0) << endl);
-        //        for (int j=0;j<2;j++) {
-        //          Log() << Verbose(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
-        //          for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
-        //            Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
-        //          Log() << Verbose(0) << endl;
-        //        }
+        Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: ";
+        for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
+          Log() << Verbose(0) << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t";
+        Log() << Verbose(0) << endl;
+//        for (int j=0;j<2;j++) {
+//          Log() << Verbose(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
+//          for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
+//            Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
+//          Log() << Verbose(0) << endl;
+//        }
       }
-      triangle->lines[i] = NULL; // free'd or not: disconnect
+      triangle->lines[i] = NULL;  // free'd or not: disconnect
     } else
-      DoeLog(1) && (eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl);
+      eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl;
   }
 
   if (TrianglesOnBoundary.erase(triangle->Nr))
-    DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl);
-  delete (triangle);
-}
-;
+    Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl;
+  delete(triangle);
+};
 
 /** Removes a line from the tesselation.
@@ -1957 +1476 @@
 void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   int Numbers[2];
 
@@ -1978 +1497 @@
         for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
           if ((*Runner).second == line) {
-            DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);
+            Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;
             line->endpoints[i]->lines.erase(Runner);
             break;
@@ -1984 +1503 @@
       } else { // there's just a single line left
         if (line->endpoints[i]->lines.erase(line->Nr))
-          DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);
+          Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;
       }
       if (line->endpoints[i]->lines.empty()) {
-        DoLog(0) && (Log() << Verbose(0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl);
+        Log() << Verbose(0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         RemoveTesselationPoint(line->endpoints[i]);
       } else {
-        DoLog(0) && (Log() << Verbose(0) << *line->endpoints[i] << " has still lines it's attached to: ");
-        for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
-          DoLog(0) && (Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t");
-        DoLog(0) && (Log() << Verbose(0) << endl);
+        Log() << Verbose(0) << *line->endpoints[i] << " has still lines it's attached to: ";
+        for(LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
+          Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
+        Log() << Verbose(0) << endl;
       }
-      line->endpoints[i] = NULL; // free'd or not: disconnect
+      line->endpoints[i] = NULL;  // free'd or not: disconnect
     } else
-      DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl);
+      eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl;
   }
   if (!line->triangles.empty())
-    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl);
+    eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl;
 
   if (LinesOnBoundary.erase(line->Nr))
-    DoLog(0) && (Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl);
-  delete (line);
-}
-;
+    Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl;
+  delete(line);
+};
 
 /** Removes a point from the tesselation.
@@ -2015 +1533 @@
 void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   if (point == NULL)
     return;
   if (PointsOnBoundary.erase(point->Nr))
-    DoLog(0) && (Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl);
-  delete (point);
-}
-;
-
-/** Checks validity of a given sphere of a candidate line.
- * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
- * We check CandidateForTesselation::OtherOptCenter
- * \param &CandidateLine contains other degenerated candidates which we have to subtract as well
- * \param RADIUS radius of sphere
- * \param *LC LinkedCell structure with other atoms
- * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
- */
-bool Tesselation::CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) const
-{
-  Info FunctionInfo(__func__);
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
-  bool flag = true;
-
-  DoLog(1) && (Log() << Verbose(1) << "Check by: draw sphere {" << CandidateLine.OtherOptCenter.x[0] << " " << CandidateLine.OtherOptCenter.x[1] << " " << CandidateLine.OtherOptCenter.x[2] << "} radius " << RADIUS << " resolution 30" << endl);
-  // get all points inside the sphere
-  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);
-
-  DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
-  for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl);
-
-  // remove triangles's endpoints
-  for (int i = 0; i < 2; i++)
-    ListofPoints->remove(CandidateLine.BaseLine->endpoints[i]->node);
-
-  // remove other candidates
-  for (TesselPointList::const_iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); ++Runner)
-    ListofPoints->remove(*Runner);
-
-  // check for other points
-  if (!ListofPoints->empty()) {
-    DoLog(1) && (Log() << Verbose(1) << "CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
-    flag = false;
-    DoLog(1) && (Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
-    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl);
-  }
-  delete (ListofPoints);
-
-  return flag;
-}
-;
+    Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl;
+  delete(point);
+};
 
 /** Checks whether the triangle consisting of the three points is already present.
@@ -2080 +1552 @@
 int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   int adjacentTriangleCount = 0;
   class BoundaryPointSet *Points[3];
@@ -2102 +1574 @@
           for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
             TriangleMap *triangles = &FindLine->second->triangles;
-            DoLog(1) && (Log() << Verbose(1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl);
+            Log() << Verbose(1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl;
             for (TriangleMap::const_iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
               if (FindTriangle->second->IsPresentTupel(Points)) {
@@ -2108 +1580 @@
               }
             }
-            DoLog(1) && (Log() << Verbose(1) << "end." << endl);
+            Log() << Verbose(1) << "end." << endl;
           }
           // Only one of the triangle lines must be considered for the triangle count.
@@ -2118 +1590 @@
   }
 
-  DoLog(0) && (Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl);
+  Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
   return adjacentTriangleCount;
-}
-;
+};
 
 /** Checks whether the triangle consisting of the three points is already present.
@@ -2133 +1604 @@
 class BoundaryTriangleSet * Tesselation::GetPresentTriangle(TesselPoint *Candidates[3])
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   class BoundaryTriangleSet *triangle = NULL;
   class BoundaryPointSet *Points[3];
@@ -2171 +1642 @@
 
   return triangle;
-}
-;
+};
+
 
 /** Finds the starting triangle for FindNonConvexBorder().
@@ -2181 +1652 @@
  * \param RADIUS radius of virtual rolling sphere
  * \param *LC LinkedCell structure with neighbouring TesselPoint's
- * \return true - a starting triangle has been created, false - no valid triple of points found
- */
-bool Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC)
-{
-  Info FunctionInfo(__func__);
+ */
+void Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC)
+{
+        Info FunctionInfo(__func__);
   int i = 0;
   TesselPoint* MaxPoint[NDIM];
   TesselPoint* Temporary;
   double maxCoordinate[NDIM];
-  BoundaryLineSet *BaseLine = NULL;
+  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++) {
@@ -2207 +1674 @@
 
   // 1. searching topmost point with respect to each axis
-  for (int i = 0; i < NDIM; i++) { // each axis
-    LC->n[i] = LC->N[i] - 1; // current axis is topmost cell
-    for (LC->n[(i + 1) % NDIM] = 0; LC->n[(i + 1) % NDIM] < LC->N[(i + 1) % NDIM]; LC->n[(i + 1) % NDIM]++)
-      for (LC->n[(i + 2) % NDIM] = 0; LC->n[(i + 2) % NDIM] < LC->N[(i + 2) % NDIM]; LC->n[(i + 2) % NDIM]++) {
-        const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
+  for (int i=0;i<NDIM;i++) { // each axis
+    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
+    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
+      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
+        const LinkedNodes *List = LC->GetCurrentCell();
         //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
         if (List != NULL) {
-          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+          for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) {
             if ((*Runner)->node->x[i] > maxCoordinate[i]) {
-              DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl);
+              Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
               maxCoordinate[i] = (*Runner)->node->x[i];
               MaxPoint[i] = (*Runner);
@@ -2222 +1689 @@
           }
         } else {
-          DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
+          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
         }
       }
   }
 
-  DoLog(1) && (Log() << Verbose(1) << "Found maximum coordinates: ");
-  for (int i = 0; i < NDIM; i++)
-    DoLog(0) && (Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t");
-  DoLog(0) && (Log() << Verbose(0) << endl);
+  Log() << Verbose(1) << "Found maximum coordinates: ";
+  for (int i=0;i<NDIM;i++)
+    Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t";
+  Log() << Verbose(0) << endl;
 
   BTS = NULL;
-  for (int k = 0; k < NDIM; k++) {
-    NormalVector.Zero();
-    NormalVector.x[k] = 1.;
-    BaseLine = new BoundaryLineSet();
-    BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
-    DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl);
+  for (int k=0;k<NDIM;k++) {
+    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;
 
     double ShortestAngle;
     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.
 
-    Temporary = NULL;
-    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?
-      delete BaseLine;
+    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();
-    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)
+    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);
+    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)
 
     // 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
-    DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl);
-    DoLog(0) && (Log() << Verbose(0) << "Found second point is at " << *BaseLine->endpoints[1]->node << ".\n");
+    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
+    Log() << Verbose(0) << "Found second point is at " << *BaseLine.endpoints[1]->node << ".\n";
 
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    CandidateForTesselation OptCandidates(BaseLine);
-    FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
-    DoLog(0) && (Log() << Verbose(0) << "List of third Points is:" << endl);
+    CandidateForTesselation OptCandidates(&BaseLine);
+    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++) {
-      DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl);
-    }
-    if (!OptCandidates.pointlist.empty()) {
-      BTS = NULL;
-      AddCandidatePolygon(OptCandidates, RADIUS, LC);
-    } else {
-      delete BaseLine;
-      continue;
-    }
-
-    if (BTS != NULL) { // we have created one starting triangle
-      delete BaseLine;
+        Log() << Verbose(0) << " " << *(*it) << endl;
+    }
+
+    BTS = NULL;
+    AddCandidateTriangle(OptCandidates);
+//    delete(BaseLine.endpoints[0]);
+//    delete(BaseLine.endpoints[1]);
+
+    if (BTS != NULL) // we have created one starting triangle
       break;
-    } else {
+    else {
       // remove all candidates from the list and then the list itself
       OptCandidates.pointlist.clear();
     }
-    delete BaseLine;
-  }
-
-  return (BTS != NULL);
-}
-;
+  }
+};
 
 /** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates.
@@ -2378 +1828 @@
 //            if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
 //              // rotated the wrong way!
-//              DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
+//              eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
 //            }
 //
@@ -2435 +1885 @@
 //          }
 //        } else {
-//          DoeLog(2) && (eLog()<< Verbose(2) << "Baseline is connected to two triangles already?" << endl);
+//          eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl;
 //        }
 //      } else {
@@ -2442 +1892 @@
 //    }
 //  } else {
-//    DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl);
+//    eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl;
 //  }
 //
@@ -2456 +1906 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
-{
-  Info FunctionInfo(__func__);
+bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
+{
+        Info FunctionInfo(__func__);
+  bool result = true;
+
   Vector CircleCenter;
   Vector CirclePlaneNormal;
-  Vector RelativeSphereCenter;
+  Vector OldSphereCenter;
   Vector SearchDirection;
   Vector helper;
-  BoundaryPointSet *ThirdPoint = NULL;
+  TesselPoint *ThirdNode = NULL;
   LineMap::iterator testline;
   double radius, CircleRadius;
 
-  for (int i = 0; i < 3; i++)
-    if ((T.endpoints[i] != CandidateLine.BaseLine->endpoints[0]) && (T.endpoints[i] != CandidateLine.BaseLine->endpoints[1])) {
-      ThirdPoint = T.endpoints[i];
-      break;
-    }
-  DoLog(0) && (Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdPoint " << *ThirdPoint << " of triangle " << T << "." << endl);
-
-  CandidateLine.T = &T;
+  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))
+      ThirdNode = T.endpoints[i]->node;
 
   // construct center of circle
@@ -2488 +1936 @@
   // calculate squared radius of circle
   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.;
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
-    DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
-
-    DoLog(1) && (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(ThirdPoint->node->node);
+    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);
+    helper.SubtractVector(ThirdNode->node);
     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
-    DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
-    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+    SearchDirection.Normalize();
+    Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
-      DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
+      eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
     }
 
     // add third point
-    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdPoint, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
 
   } else {
-    DoLog(0) && (Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl);
+    Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl;
   }
 
   if (CandidateLine.pointlist.empty()) {
-    DoeLog(2) && (eLog() << Verbose(2) << "Could not find a suitable candidate." << endl);
+    eLog() << Verbose(2) << "Could not find a suitable candidate." << endl;
     return false;
   }
-  DoLog(0) && (Log() << Verbose(0) << "Third Points are: " << endl);
+  Log() << Verbose(0) << "Third Points are: " << endl;
   for (TesselPointList::iterator it = CandidateLine.pointlist.begin(); it != CandidateLine.pointlist.end(); ++it) {
-    DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl);
+    Log() << Verbose(0) << " " << *(*it) << endl;
   }
 
   return true;
-}
-;
-
-/** Walks through Tesselation::OpenLines() and finds candidates for newly created ones.
- * \param *&LCList atoms in LinkedCell list
- * \param RADIUS radius of the virtual sphere
- * \return true - for all open lines without candidates so far, a candidate has been found,
- *         false - at least one open line without candidate still
- */
-bool Tesselation::FindCandidatesforOpenLines(const double RADIUS, const LinkedCell *&LCList)
-{
-  bool TesselationFailFlag = true;
-  CandidateForTesselation *baseline = NULL;
-  BoundaryTriangleSet *T = NULL;
-
-  for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) {
-    baseline = Runner->second;
-    if (baseline->pointlist.empty()) {
-      assert((baseline->BaseLine->triangles.size() == 1) && ("Open line without exactly one attached triangle"));
-      T = (((baseline->BaseLine->triangles.begin()))->second);
-      DoLog(1) && (Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl);
-      TesselationFailFlag = TesselationFailFlag && FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
-    }
-  }
-  return TesselationFailFlag;
-}
-;
+
+//  BoundaryLineSet *BaseRay = CandidateLine.BaseLine;
+//  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
+//    Log() << Verbose(0) << "Third point candidate is " << *(*it)->point
+//    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
+//    Log() << Verbose(0) << "Baseline is " << *BaseRay << endl;
+//
+//    // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+//    TesselPoint *PointCandidates[3];
+//    PointCandidates[0] = (*it)->point;
+//    PointCandidates[1] = BaseRay->endpoints[0]->node;
+//    PointCandidates[2] = BaseRay->endpoints[1]->node;
+//    int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates);
+//
+//    BTS = NULL;
+//    // check for present edges and whether we reach better candidates from them
+//    //if (HasOtherBaselineBetterCandidate(BaseRay, (*it)->point, ShortestAngle, RADIUS, LC) ) {
+//    if (0) {
+//      result = false;
+//      break;
+//    } else {
+//      // If there is no triangle, add it regularly.
+//      if (existentTrianglesCount == 0) {
+//        AddTesselationPoint((*it)->point, 0);
+//        AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
+//        AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
+//
+//        if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
+//          CandidateLine.point = (*it)->point;
+//          CandidateLine.OptCenter.CopyVector(&((*it)->OptCenter));
+//          CandidateLine.OtherOptCenter.CopyVector(&((*it)->OtherOptCenter));
+//          CandidateLine.ShortestAngle = ShortestAngle;
+//        } else {
+////          eLog() << Verbose(1) << "This triangle consisting of ";
+////          Log() << Verbose(0) << *(*it)->point << ", ";
+////          Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
+////          Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
+////          Log() << Verbose(0) << "exists and is not added, as it 0x80000000006fc150(does not seem helpful!" << endl;
+//          result = false;
+//        }
+//      } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
+//          AddTesselationPoint((*it)->point, 0);
+//          AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
+//          AddTesselationPoint(BaseRay->endpoints[1]->node, 2);
+//
+//          // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
+//          // i.e. at least one of the three lines must be present with TriangleCount <= 1
+//          if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS) || CandidateLine.BaseLine->skipped) {
+//            CandidateLine.point = (*it)->point;
+//            CandidateLine.OptCenter.CopyVector(&(*it)->OptCenter);
+//            CandidateLine.OtherOptCenter.CopyVector(&(*it)->OtherOptCenter);
+//            CandidateLine.ShortestAngle = ShortestAngle+2.*M_PI;
+//
+//          } else {
+////            eLog() << Verbose(1) << "This triangle consisting of " << *(*it)->point << ", " << *BaseRay->endpoints[0]->node << " and " << *BaseRay->endpoints[1]->node << " " << "exists and is not added, as it does not seem helpful!" << endl;
+//            result = false;
+//          }
+//      } else {
+////        Log() << Verbose(1) << "This triangle consisting of ";
+////        Log() << Verbose(0) << *(*it)->point << ", ";
+////        Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and ";
+////        Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " ";
+////        Log() << Verbose(0) << "is invalid!" << endl;
+//        result = false;
+//      }
+//    }
+//
+//    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
+//    BaseRay = BLS[0];
+//    if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) {
+//      eLog() << Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl;
+//      exit(255);
+//    }
+//
+//  }
+//
+//  // remove all candidates from the list and then the list itself
+//  class CandidateForTesselation *remover = NULL;
+//  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
+//    remover = *it;
+//    delete(remover);
+//  }
+//  delete(OptCandidates);
+  return result;
+};
 
 /** Adds the present line and candidate point from \a &CandidateLine to the Tesselation.
  * \param CandidateLine triangle to add
- * \param RADIUS Radius of sphere
- * \param *LC LinkedCell structure
- * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in
- * AddTesselationLine() in AddCandidateTriangle()
- */
-void Tesselation::AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC)
-{
-  Info FunctionInfo(__func__);
+ * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in AddTesselationLine()
+ */
+void Tesselation::AddCandidateTriangle(CandidateForTesselation CandidateLine)
+{
+        Info FunctionInfo(__func__);
   Vector Center;
   TesselPoint * const TurningPoint = CandidateLine.BaseLine->endpoints[0]->node;
-  TesselPointList::iterator Runner;
-  TesselPointList::iterator Sprinter;
 
   // 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);
-
-  DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl);
-  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
-    DoLog(0) && (Log() << Verbose(0) << " " << **TesselRunner << endl);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
 
   // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
-  Runner = connectedClosestPoints->begin();
-  Sprinter = Runner;
+  TesselPointList::iterator Runner = connectedClosestPoints->begin();
+  TesselPointList::iterator Sprinter = Runner;
   Sprinter++;
-  while (Sprinter != connectedClosestPoints->end()) {
-    DoLog(0) && (Log() << Verbose(0) << "Current Runner is " << *(*Runner) << " and sprinter is " << *(*Sprinter) << "." << endl);
-
+  while(Sprinter != connectedClosestPoints->end()) {
+    // add the points
     AddTesselationPoint(TurningPoint, 0);
-    AddTesselationPoint(*Runner, 1);
-    AddTesselationPoint(*Sprinter, 2);
-
-    AddCandidateTriangle(CandidateLine, Opt);
-
+    AddTesselationPoint((*Runner), 1);
+    AddTesselationPoint((*Sprinter), 2);
+
+    Center.CopyVector(&CandidateLine.OptCenter);
+    // add the lines
+    AddTesselationLine(TPS[0], TPS[1], 0);
+    AddTesselationLine(TPS[0], TPS[2], 1);
+    AddTesselationLine(TPS[1], TPS[2], 2);
+
+    // add the triangles
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    AddTesselationTriangle();
+    Center.Scale(-1.);
+    BTS->GetNormalVector(Center);
+
+    Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << "." << endl;
     Runner = Sprinter;
     Sprinter++;
-    if (Sprinter != connectedClosestPoints->end()) {
-      // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
-      FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OptCenter); // Assume BTS contains last triangle
-      DoLog(0) && (Log() << Verbose(0) << " There are still more triangles to add." << endl);
-    }
-    // pick candidates for other open lines as well
-    FindCandidatesforOpenLines(RADIUS, LC);
-
-    // check whether we add a degenerate or a normal triangle
-    if (CheckDegeneracy(CandidateLine, RADIUS, LC)) {
-      // add normal and degenerate triangles
-      DoLog(1) && (Log() << Verbose(1) << "Triangle of endpoints " << *TPS[0] << "," << *TPS[1] << " and " << *TPS[2] << " is degenerated, adding both sides." << endl);
-      AddCandidateTriangle(CandidateLine, OtherOpt);
-
-      if (Sprinter != connectedClosestPoints->end()) {
-        // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
-        FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OtherOptCenter);
-      }
-      // pick candidates for other open lines as well
-      FindCandidatesforOpenLines(RADIUS, LC);
-    }
-  }
-  delete (connectedClosestPoints);
-};
-
-/** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate.
- * \param *Sprinter next candidate to which internal open lines are set
- * \param *OptCenter OptCenter for this candidate
- */
-void Tesselation::FindDegeneratedCandidatesforOpenLines(TesselPoint * const Sprinter, const Vector * const OptCenter)
-{
-  Info FunctionInfo(__func__);
-
-  pair<LineMap::iterator, LineMap::iterator> FindPair = TPS[0]->lines.equal_range(TPS[2]->node->nr);
-  for (LineMap::const_iterator FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl);
-    // If there is a line with less than two attached triangles, we don't need a new line.
-    if (FindLine->second->triangles.size() == 1) {
-      CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
-      if (!Finder->second->pointlist.empty())
-        DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl);
-      else {
-        DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter) << endl);
-        Finder->second->T = BTS;  // is last triangle
-        Finder->second->pointlist.push_back(Sprinter);
-        Finder->second->ShortestAngle = 0.;
-        Finder->second->OptCenter.CopyVector(OptCenter);
-      }
-    }
-  }
-};
-
-/** If a given \a *triangle is degenerated, this adds both sides.
- * i.e. the triangle with same BoundaryPointSet's but NormalVector in opposite direction.
- * Note that endpoints are stored in Tesselation::TPS
- * \param CandidateLine CanddiateForTesselation structure for the desired BoundaryLine
- * \param RADIUS radius of sphere
- * \param *LC pointer to LinkedCell structure
- */
-void Tesselation::AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC)
-{
-  Info FunctionInfo(__func__);
-  Vector Center;
-  CandidateMap::const_iterator CandidateCheck = OpenLines.end();
-  BoundaryTriangleSet *triangle = NULL;
-
-  /// 1. Create or pick the lines for the first triangle
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for first triangle ..." << endl);
-  for (int i = 0; i < 3; i++) {
-    BLS[i] = NULL;
-    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
-    AddTesselationLine(&CandidateLine.OptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
-  }
-
-  /// 2. create the first triangle and NormalVector and so on
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding first triangle with center at " << CandidateLine.OptCenter << " ..." << endl);
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  AddTesselationTriangle();
-
-  // create normal vector
-  BTS->GetCenter(&Center);
-  Center.SubtractVector(&CandidateLine.OptCenter);
-  BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
-  BTS->GetNormalVector(Center);
-  // give some verbose output about the whole procedure
-  if (CandidateLine.T != NULL)
-    DoLog(0) && (Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
-  else
-    DoLog(0) && (Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
-  triangle = BTS;
-
-  /// 3. Gather candidates for each new line
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding candidates to new lines ..." << endl);
-  for (int i = 0; i < 3; i++) {
-    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
-    CandidateCheck = OpenLines.find(BLS[i]);
-    if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
-      if (CandidateCheck->second->T == NULL)
-        CandidateCheck->second->T = triangle;
-      FindNextSuitableTriangle(*(CandidateCheck->second), *CandidateCheck->second->T, RADIUS, LC);
-    }
-  }
-
-  /// 4. Create or pick the lines for the second triangle
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for second triangle ..." << endl);
-  for (int i = 0; i < 3; i++) {
-    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
-    AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
-  }
-
-  /// 5. create the second triangle and NormalVector and so on
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangle with center at " << CandidateLine.OtherOptCenter << " ..." << endl);
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  AddTesselationTriangle();
-
-  BTS->SphereCenter.CopyVector(&CandidateLine.OtherOptCenter);
-  // create normal vector in other direction
-  BTS->GetNormalVector(&triangle->NormalVector);
-  BTS->NormalVector.Scale(-1.);
-  // give some verbose output about the whole procedure
-  if (CandidateLine.T != NULL)
-    DoLog(0) && (Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
-  else
-    DoLog(0) && (Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
-
-  /// 6. Adding triangle to new lines
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangles to new lines ..." << endl);
-  for (int i = 0; i < 3; i++) {
-    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
-    CandidateCheck = OpenLines.find(BLS[i]);
-    if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
-      if (CandidateCheck->second->T == NULL)
-        CandidateCheck->second->T = BTS;
-    }
-  }
-}
-;
-
-/** Adds a triangle to the Tesselation structure from three given TesselPoint's.
- * Note that endpoints are in Tesselation::TPS.
- * \param CandidateLine CandidateForTesselation structure contains other information
- * \param type which opt center to add (i.e. which side) and thus which NormalVector to take
- */
-void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine, enum centers type)
-{
-  Info FunctionInfo(__func__);
-  Vector Center;
-  Vector *OptCenter = (type == Opt) ? &CandidateLine.OptCenter : &CandidateLine.OtherOptCenter;
-
-  // add the lines
-  AddTesselationLine(OptCenter, TPS[2], TPS[0], TPS[1], 0);
-  AddTesselationLine(OptCenter, TPS[1], TPS[0], TPS[2], 1);
-  AddTesselationLine(OptCenter, TPS[0], TPS[1], TPS[2], 2);
-
-  // add the triangles
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  AddTesselationTriangle();
-
-  // create normal vector
-  BTS->GetCenter(&Center);
-  Center.SubtractVector(OptCenter);
-  BTS->SphereCenter.CopyVector(OptCenter);
-  BTS->GetNormalVector(Center);
-
-  // give some verbose output about the whole procedure
-  if (CandidateLine.T != NULL)
-    DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
-  else
-    DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
-}
-;
+  }
+};
 
 /** Checks whether the quadragon of the two triangles connect to \a *Base is convex.
@@ -2778 +2123 @@
 class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   class BoundaryPointSet *Spot = NULL;
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint;
 
-  int m = 0;
-  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
+  int m=0;
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
       if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
         BPS[m++] = runner->second->endpoints[j];
-  OtherBase = new class BoundaryLineSet(BPS, -1);
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Other base line is " << *OtherBase << "." << endl);
+  OtherBase = new class BoundaryLineSet(BPS,-1);
+
+  Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl;
+  Log() << Verbose(1) << "INFO: Other base line is " << *OtherBase << "." << endl;
 
   // get the closest point on each line to the other line
@@ -2797 +2142 @@
 
   // delete the temporary other base line
-  delete (OtherBase);
+  delete(OtherBase);
 
   // get the distance vector from Base line to OtherBase line
@@ -2804 +2149 @@
   BaseLine.CopyVector(Base->endpoints[1]->node->node);
   BaseLine.SubtractVector(Base->endpoints[0]->node->node);
-  for (int i = 0; i < 2; i++) {
+  for (int i=0;i<2;i++) {
     DistanceToIntersection[i].CopyVector(ClosestPoint);
     DistanceToIntersection[i].SubtractVector(Base->endpoints[i]->node->node);
     distance[i] = BaseLine.ScalarProduct(&DistanceToIntersection[i]);
   }
-  delete (ClosestPoint);
-  if ((distance[0] * distance[1]) > 0) { // have same sign?
-    DoLog(1) && (Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl);
+  delete(ClosestPoint);
+  if ((distance[0] * distance[1]) > 0)  { // have same sign?
+    Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1]  << ". " << *Base << "' rectangle is concave." << endl;
     if (distance[0] < distance[1]) {
       Spot = Base->endpoints[0];
@@ -2818 +2163 @@
     }
     return Spot;
-  } else { // different sign, i.e. we are in between
-    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl);
+  } else {  // different sign, i.e. we are in between
+    Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;
     return NULL;
   }
 
-}
-;
+};
 
 void Tesselation::PrintAllBoundaryPoints(ofstream *out) const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   // print all lines
-  DoLog(0) && (Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl);
-  for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin(); PointRunner != PointsOnBoundary.end(); PointRunner++)
-    DoLog(0) && (Log() << Verbose(0) << *(PointRunner->second) << endl);
-}
-;
+  Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl;
+  for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++)
+    Log() << Verbose(0) << *(PointRunner->second) << endl;
+};
 
 void Tesselation::PrintAllBoundaryLines(ofstream *out) const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   // print all lines
-  DoLog(0) && (Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl);
+  Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl;
   for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
-    DoLog(0) && (Log() << Verbose(0) << *(LineRunner->second) << endl);
-}
-;
+    Log() << Verbose(0) << *(LineRunner->second) << endl;
+};
 
 void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   // print all triangles
-  DoLog(0) && (Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl);
+  Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl;
   for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
-    DoLog(0) && (Log() << Verbose(0) << *(TriangleRunner->second) << endl);
-}
-;
+    Log() << Verbose(0) << *(TriangleRunner->second) << endl;
+};
 
 /** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher".
@@ -2863 +2204 @@
 double Tesselation::PickFarthestofTwoBaselines(class BoundaryLineSet *Base)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint[2];
   double volume;
 
-  int m = 0;
-  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
+  int m=0;
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
       if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
         BPS[m++] = runner->second->endpoints[j];
-  OtherBase = new class BoundaryLineSet(BPS, -1);
-
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl);
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Other base line is " << *OtherBase << "." << endl);
+  OtherBase = new class BoundaryLineSet(BPS,-1);
+
+  Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl;
+  Log() << Verbose(0) << "INFO: Other base line is " << *OtherBase << "." << endl;
 
   // get the closest point on each line to the other line
@@ -2891 +2232 @@
 
   // delete the temporary other base line and the closest points
-  delete (ClosestPoint[0]);
-  delete (ClosestPoint[1]);
-  delete (OtherBase);
+  delete(ClosestPoint[0]);
+  delete(ClosestPoint[1]);
+  delete(OtherBase);
 
   if (Distance.NormSquared() < MYEPSILON) { // check for intersection
-    DoLog(0) && (Log() << Verbose(0) << "REJECT: Both lines have an intersection: Nothing to do." << endl);
+    Log() << Verbose(0) << "REJECT: Both lines have an intersection: Nothing to do." << endl;
     return false;
   } else { // check for sign against BaseLineNormal
@@ -2902 +2243 @@
     BaseLineNormal.Zero();
     if (Base->triangles.size() < 2) {
-      DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
+      eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl;
       return 0.;
     }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
+      Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;
       BaseLineNormal.AddVector(&(runner->second->NormalVector));
     }
-    BaseLineNormal.Scale(1. / 2.);
+    BaseLineNormal.Scale(1./2.);
 
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
-      DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl);
+      Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;
       // calculate volume summand as a general tetraeder
       return volume;
-    } else { // Base higher than OtherBase -> do nothing
-      DoLog(0) && (Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl);
+    } else {  // Base higher than OtherBase -> do nothing
+      Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl;
       return 0.;
     }
   }
-}
-;
+};
 
 /** For a given baseline and its two connected triangles, flips the baseline.
@@ -2932 +2272 @@
 class BoundaryLineSet * Tesselation::FlipBaseline(class BoundaryLineSet *Base)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
   Vector BaseLineNormal;
   int OldTriangleNrs[2], OldBaseLineNr;
-  int i, m;
+  int i,m;
 
   // calculate NormalVector for later use
   BaseLineNormal.Zero();
   if (Base->triangles.size() < 2) {
-    DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
+    eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl;
     return NULL;
   }
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
+    Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;
     BaseLineNormal.AddVector(&(runner->second->NormalVector));
   }
-  BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
+  BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
 
   // get the two triangles
   // gather four endpoints and four lines
-  for (int j = 0; j < 4; j++)
+  for (int j=0;j<4;j++)
     OldLines[j] = NULL;
-  for (int j = 0; j < 2; j++)
+  for (int j=0;j<2;j++)
     OldPoints[j] = NULL;
-  i = 0;
-  m = 0;
-  DoLog(0) && (Log() << Verbose(0) << "The four old lines are: ");
-  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
+  i=0;
+  m=0;
+  Log() << Verbose(0) << "The four old lines are: ";
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
       if (runner->second->lines[j] != Base) { // pick not the central baseline
         OldLines[i++] = runner->second->lines[j];
-        DoLog(0) && (Log() << Verbose(0) << *runner->second->lines[j] << "\t");
+        Log() << Verbose(0) << *runner->second->lines[j] << "\t";
       }
-  DoLog(0) && (Log() << Verbose(0) << endl);
-  DoLog(0) && (Log() << Verbose(0) << "The two old points are: ");
-  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
+  Log() << Verbose(0) << endl;
+  Log() << Verbose(0) << "The two old points are: ";
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j=0;j<3;j++) // all of their endpoints and baselines
       if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints
         OldPoints[m++] = runner->second->endpoints[j];
-        DoLog(0) && (Log() << Verbose(0) << *runner->second->endpoints[j] << "\t");
+        Log() << Verbose(0) << *runner->second->endpoints[j] << "\t";
       }
-  DoLog(0) && (Log() << Verbose(0) << endl);
+  Log() << Verbose(0) << endl;
 
   // check whether everything is in place to create new lines and triangles
-  if (i < 4) {
-    DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
+  if (i<4) {
+    eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl;
     return NULL;
   }
-  for (int j = 0; j < 4; j++)
+  for (int j=0;j<4;j++)
     if (OldLines[j] == NULL) {
-      DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
+      eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl;
       return NULL;
     }
-  for (int j = 0; j < 2; j++)
+  for (int j=0;j<2;j++)
     if (OldPoints[j] == NULL) {
-      DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl);
+      eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl;
       return NULL;
     }
 
   // remove triangles and baseline removes itself
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting baseline " << *Base << " from global list." << endl);
+  Log() << Verbose(0) << "INFO: Deleting baseline " << *Base << " from global list." << endl;
   OldBaseLineNr = Base->Nr;
-  m = 0;
-  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-    DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl);
+  m=0;
+  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+    Log() << Verbose(0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;
     OldTriangleNrs[m++] = runner->second->Nr;
     RemoveTesselationTriangle(runner->second);
@@ -3007 +2347 @@
   NewLine = new class BoundaryLineSet(BPS, OldBaseLineNr);
   LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
-  DoLog(0) && (Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl);
+  Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl;
 
   // construct new triangles with flipped baseline
-  i = -1;
+  i=-1;
   if (OldLines[0]->IsConnectedTo(OldLines[2]))
-    i = 2;
+    i=2;
   if (OldLines[0]->IsConnectedTo(OldLines[3]))
-    i = 3;
-  if (i != -1) {
+    i=3;
+  if (i!=-1) {
     BLS[0] = OldLines[0];
     BLS[1] = OldLines[i];
@@ -3022 +2362 @@
     BTS->GetNormalVector(BaseLineNormal);
     AddTesselationTriangle(OldTriangleNrs[0]);
-    DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);
-
-    BLS[0] = (i == 2 ? OldLines[3] : OldLines[2]);
+    Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl;
+
+    BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
     BLS[1] = OldLines[1];
     BLS[2] = NewLine;
@@ -3030 +2370 @@
     BTS->GetNormalVector(BaseLineNormal);
     AddTesselationTriangle(OldTriangleNrs[1]);
-    DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);
+    Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl;
   } else {
-    DoeLog(0) && (eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl);
+    eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl;
     return NULL;
   }
 
   return NewLine;
-}
-;
+};
+
 
 /** Finds the second point of starting triangle.
@@ -3050 +2390 @@
 void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   Vector AngleCheck;
   class TesselPoint* Candidate = NULL;
@@ -3059 +2399 @@
   int Nupper[NDIM];
 
-  if (LC->SetIndexToNode(a)) { // get cell for the starting point
-    for (int i = 0; i < NDIM; i++) // store indices of this cell
+  if (LC->SetIndexToNode(a)) {  // get cell for the starting point
+    for(int i=0;i<NDIM;i++) // store indices of this cell
       N[i] = LC->n[i];
   } else {
-    DoeLog(1) && (eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl);
+    eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl;
     return;
   }
   // then go through the current and all neighbouring cells and check the contained points for possible candidates
-  for (int i = 0; i < NDIM; i++) {
-    Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
-    Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
-  }
-  DoLog(0) && (Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl);
+  for (int i=0;i<NDIM;i++) {
+    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+  }
+  Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :"
+    << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << 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 LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
+        const LinkedNodes *List = LC->GetCurrentCell();
         //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
         if (List != NULL) {
-          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             Candidate = (*Runner);
             // check if we only have one unique point yet ...
@@ -3105 +2446 @@
               norm = aCandidate.Norm();
               // second point shall have smallest angle with respect to Oben vector
-              if (norm < RADIUS * 2.) {
+              if (norm < RADIUS*2.) {
                 angle = AngleCheck.Angle(&Oben);
                 if (angle < Storage[0]) {
                   //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                  DoLog(1) && (Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n");
+                  Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
                   OptCandidate = Candidate;
                   Storage[0] = angle;
@@ -3124 +2465 @@
           }
         } else {
-          DoLog(0) && (Log() << Verbose(0) << "Linked cell list is empty." << endl);
+          Log() << Verbose(0) << "Linked cell list is empty." << endl;
         }
       }
-}
-;
+};
+
 
 /** This recursive function finds a third point, to form a triangle with two given ones.
@@ -3156 +2497 @@
  * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
  * @param CandidateLine CandidateForTesselation with the current base line and list of candidates and ShortestAngle
- * @param ThirdPoint third point to avoid in search
+ * @param ThirdNode third point to avoid in search
  * @param RADIUS radius of sphere
  * @param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
-{
-  Info FunctionInfo(__func__);
-  Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
+void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class TesselPoint  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const
+{
+        Info FunctionInfo(__func__);
+  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 NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
-  Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
-  Vector NewNormalVector; // normal vector of the Candidate's triangle
+  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+  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];
   TesselPoint *Candidate = NULL;
 
-  DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl);
-
-  // copy old center
-  CandidateLine.OldCenter.CopyVector(&OldSphereCenter);
-  CandidateLine.ThirdPoint = ThirdPoint;
-  CandidateLine.pointlist.clear();
+  Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
 
   // construct center of circle
@@ -3195 +2528 @@
   CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
 
-  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
-  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
-
-  // calculate squared radius TesselPoint *ThirdPoint,f circle
-  radius = CirclePlaneNormal.NormSquared() / 4.;
-  if (radius < RADIUS * RADIUS) {
-    CircleRadius = RADIUS * RADIUS - radius;
+  // calculate squared radius TesselPoint *ThirdNode,f circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
-    DoLog(1) && (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) {
-      DoeLog(1) && (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) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl);
+      //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 
       // check SearchDirection
-      DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
-      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
-        DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);
+      //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;
       }
 
       // get cell for the starting point
       if (LC->SetIndexToVector(&CircleCenter)) {
-        for (int i = 0; i < NDIM; i++) // store indices of this cell
-          N[i] = LC->n[i];
+        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;
       } else {
-        DoeLog(1) && (eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl);
+        eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
         return;
       }
       // then go through the current and all neighbouring cells and check the contained points for possible candidates
       //Log() << Verbose(1) << "LC Intervals:";
-      for (int i = 0; i < NDIM; i++) {
-        Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
-        Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
+      for (int i=0;i<NDIM;i++) {
+        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
         //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
       }
@@ -3240 +2570 @@
         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 LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
+            const LinkedNodes *List = LC->GetCurrentCell();
             //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
             if (List != NULL) {
-              for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+              for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
                 Candidate = (*Runner);
 
                 // check for three unique points
-                DoLog(2) && (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)) {
-
-                  // find center on the plane
-                  GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
-                  DoLog(1) && (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)) {
-                    DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl);
-                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
-                    DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
-                    DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
-                    DoLog(1) && (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) {
-                        // construct both new centers
-                        NewSphereCenter.CopyVector(&NewPlaneCenter);
-                        OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
-                        helper.CopyVector(&NewNormalVector);
-                        helper.Scale(sqrt(RADIUS * RADIUS - radius));
-                        DoLog(2) && (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);
-                        DoLog(2) && (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);
-                        DoLog(2) && (Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl);
-
-                        alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                        Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                        if ((ThirdPoint != NULL) && (Candidate == ThirdPoint->node)) { // in that case only the other circlecenter is valid
-                          if (OldSphereCenter.DistanceSquared(&NewSphereCenter) < OldSphereCenter.DistanceSquared(&OtherNewSphereCenter))
-                            alpha = Otheralpha;
-                        } else
-                          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
-                        if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
-                          if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                            CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
-                            CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
-                          } else {
-                            CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
-                            CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
-                          }
-                          // if there is an equal candidate, add it to the list without clearing the list
-                          if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
-                            CandidateLine.pointlist.push_back(Candidate);
-                            DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl);
-                          } else {
-                            // remove all candidates from the list and then the list itself
-                            CandidateLine.pointlist.clear();
-                            CandidateLine.pointlist.push_back(Candidate);
-                            DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl);
-                          }
-                          CandidateLine.ShortestAngle = alpha;
-                          DoLog(0) && (Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << 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) ){
+
+                  // 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(&NewSphereCenter);
+                    if (radius < RADIUS*RADIUS) {
+                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                      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;
+
+                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      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
+                      if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
+                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
+                          CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
+                          CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
                         } else {
-                          if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
-                            DoLog(1) && (Log() << Verbose(1) << "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl);
-                          } else {
-                            DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl);
-                          }
+                          CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
+                          CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
                         }
+                        // if there is an equal candidate, add it to the list without clearing the list
+                        if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
+                          CandidateLine.pointlist.push_back(Candidate);
+                          Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with "
+                            << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                        } else {
+                          // remove all candidates from the list and then the list itself
+                          CandidateLine.pointlist.clear();
+                          CandidateLine.pointlist.push_back(Candidate);
+                          Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with "
+                            << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                        }
+                        CandidateLine.ShortestAngle = alpha;
+                        Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl;
                       } else {
-                        DoLog(1) && (Log() << Verbose(1) << "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius) << endl);
+                        if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
+                          Log() << Verbose(1) << "REJECT: Old candidate " << *(Candidate) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl;
+                        } else {
+                          Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                        }
                       }
+
                     } else {
-                      DoLog(1) && (Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl);
+                      Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
                     }
                   } else {
-                    DoLog(1) && (Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl);
+                    Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
                   }
                 } else {
-                  if (ThirdPoint != NULL) {
-                    DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << "." << endl);
+                  if (ThirdNode != NULL) {
+                    Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
                   } else {
-                    DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl);
+                    Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl;
                   }
                 }
@@ -3333 +2655 @@
           }
     } else {
-      DoeLog(1) && (eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
+      eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
     }
   } else {
-    if (ThirdPoint != NULL)
-      DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!" << endl);
+    if (ThirdNode != NULL)
+      Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
     else
-      DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl);
-  }
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl);
+      Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl;
+  }
+
+  Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl;
   if (CandidateLine.pointlist.size() > 1) {
     CandidateLine.pointlist.unique();
     CandidateLine.pointlist.sort(); //SortCandidates);
   }
-
-  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
-    DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl);
-    performCriticalExit();
-  }
-}
-;
+};
 
 /** Finds the endpoint two lines are sharing.
@@ -3362 +2678 @@
 class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   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
-    for (int j = 0; j < 2; j++) { // for both endpoints
-      OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
-      if (!OrderTest.second) { // if insertion fails, we have common endpoint
-        node = OrderTest.first->second;
-        DoLog(1) && (Log() << Verbose(1) << "Common endpoint of lines " << *line1 << " and " << *line2 << " is: " << *node << "." << endl);
-        j = 2;
-        i = 2;
-        break;
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
+        if (!OrderTest.second)
+          { // if insertion fails, we have common endpoint
+            node = OrderTest.first->second;
+            Log() << Verbose(1) << "Common endpoint of lines " << *line1
+                << " and " << *line2 << " is: " << *node << "." << endl;
+            j = 2;
+            i = 2;
+            break;
+          }
       }
-    }
   return node;
-}
-;
-
-/** 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()) {
-    DoeLog(1) && (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];
-  DoLog(1) && (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 LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
-        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedCell::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));
-              DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
-            }
-          }
-        } else {
-          DoeLog(1) && (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()) {
-    DoeLog(1) && (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) {
-    DoeLog(1) && (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
-  DoLog(1) && (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;
-          DoLog(1) && (Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl);
-        } else {
-          DoLog(1) && (Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl);
-        }
-      } else {
-        DoLog(1) && (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) {
-    DoLog(0) && (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) {
-    DoeLog(1) && (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
-  DoLog(1) && (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;
-          DoLog(1) && (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);
-          DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl);
-        } else { // line is worse
-          DoLog(1) && (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);
+  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;
+  }
+  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;
         }
-      } 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()) {
-          DoeLog(0) && (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;
-          DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl);
-        } else {
-          DoLog(2) && (Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl);
-        }
+        //Log() << Verbose(1) << "List of triangle points:" << endl;
+        //Log() << Verbose(2) << *trianglePoints[i] << endl;
       }
-    }
-  }
-  delete (points);
-
-  // check whether closest line is "too close" :), then it's inside
-  if (ClosestLines.empty()) {
-    DoLog(0) && (Log() << Verbose(0) << "Is the only point, no one else is closeby." << 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;
-  }
-  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;
-}
-;
+  } else
+    return triangles;
+};
 
 /** Finds closest triangle to a point.
@@ -3593 +2789 @@
  * \param *out output stream for debugging
  * \param *x Vector to look from
- * \param &distance contains found distance on return
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
-{
-  Info FunctionInfo(__func__);
+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;
-      DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl);
-    } else {
-      DoLog(1) && (Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl);
-    }
-  }
-  delete (triangles);
-
+  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
-{
-  Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point, this, LC);
-
-  return Intersections.IsInside();
-}
-;
-
-/** 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
- */
-double Tesselation::GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const
-{
-  Info FunctionInfo(__func__);
+ * @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 (triangle == NULL) {// is boundary point or only point in point cloud?
-    DoLog(1) && (Log() << Verbose(1) << "No triangle given!" << endl);
-    return -1.;
+
+  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;
+  }
+
+  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;
   } else {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl);
-  }
-
-  triangle->GetCenter(&Center);
-  DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
-  DistanceToCenter.CopyVector(&Center);
-  DistanceToCenter.SubtractVector(&Point);
-  DoLog(2) && (Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl);
-
-  // check whether we are on boundary
-  if (fabs(DistanceToCenter.ScalarProduct(&triangle->NormalVector)) < MYEPSILON) {
-    // calculate whether inside of triangle
-    DistanceToCenter.CopyVector(&Point);
-    Center.CopyVector(&Point);
-    Center.SubtractVector(&triangle->NormalVector); // points towards MolCenter
-    DistanceToCenter.AddVector(&triangle->NormalVector); // points outside
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl);
-    if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
-      DoLog(1) && (Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl);
-      return 0.;
-    } else {
-      DoLog(1) && (Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl);
-      return false;
-    }
-  } else {
-    // calculate smallest distance
-    distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection);
-    DoLog(1) && (Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl);
-
-    // then check direction to boundary
-    if (DistanceToCenter.ScalarProduct(&triangle->NormalVector) > MYEPSILON) {
-      DoLog(1) && (Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl);
-      return -distance;
-    } else {
-      DoLog(1) && (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;
   }
 }
-;
-
-/** Calculates minimum distance from \a&Point to a tesselated surface.
- * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
- * \param &Point point to calculate distance from
- * \param *LC needed for finding closest points fast
- * \return distance squared to closest point on surface
- */
-double Tesselation::GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const
-{
-  Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point, this, LC);
-
-  return Intersections.GetSmallestDistance();
+
+/** 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);
 }
-;
-
-/** Calculates minimum distance from \a&Point to a tesselated surface.
- * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
- * \param &Point point to calculate distance from
- * \param *LC needed for finding closest points fast
- * \return distance squared to closest point on surface
- */
-BoundaryTriangleSet * Tesselation::GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const
-{
-  Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point, this, LC);
-
-  return Intersections.GetClosestTriangle();
-}
-;
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -3753 +2871 @@
  * @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;
   bool takePoint = false;
+
   // find the respective boundary point
   PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
@@ -3765 +2884 @@
     ReferencePoint = PointRunner->second;
   } else {
-    DoeLog(2) && (eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
+    eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
     ReferencePoint = NULL;
   }
@@ -3771 +2890 @@
   // little trick so that we look just through lines connect to the BoundaryPoint
   // OR fall-back to look through all lines if there is no such BoundaryPoint
-  const LineMap *Lines;
-  ;
+  const LineMap *Lines;;
   if (ReferencePoint != NULL)
     Lines = &(ReferencePoint->lines);
@@ -3779 +2897 @@
   LineMap::const_iterator findLines = Lines->begin();
   while (findLines != Lines->end()) {
-    takePoint = false;
-
-    if (findLines->second->endpoints[0]->Nr == Point->nr) {
-      takePoint = true;
-      current = findLines->second->endpoints[1]->node;
-    } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
-      takePoint = true;
-      current = findLines->second->endpoints[0]->node;
-    }
-
-    if (takePoint) {
-      DoLog(1) && (Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl);
-      connectedPoints->insert(current);
-    }
-
-    findLines++;
-  }
-
-  if (connectedPoints->empty()) { // if have not found any points
-    DoeLog(1) && (eLog() << Verbose(1) << "We have not found any connected points to " << *Point << "." << endl);
+   takePoint = false;
+
+   if (findLines->second->endpoints[0]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[1]->node;
+   } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[0]->node;
+   }
+
+   if (takePoint) {
+     Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
+     connectedPoints->insert(current);
+   }
+
+   findLines++;
+  }
+
+  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;
   }
 
   return connectedPoints;
-}
-;
+};
+
 
 /** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
@@ -3818 +2936 @@
  * @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
-{
-  Info FunctionInfo(__func__);
+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) {
-    DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);
-    delete (connectedCircle);
+    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
+    delete(connectedCircle);
     return NULL;
   }
 
   // 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 {
-    DoeLog(0) && (eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);
-    performCriticalExit();
-  }
-  PlaneNormal.Scale(1.0 / triangles->size());
-  DoLog(1) && (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
@@ -3855 +2973 @@
     AngleZero.ProjectOntoPlane(&PlaneNormal);
   }
-  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
-    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
+  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) {
-      DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
+      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
       performCriticalExit();
     }
   }
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);
+  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);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);
+  Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
 
   // 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);
     helper.ProjectOntoPlane(&PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
-    DoLog(0) && (Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl);
-    anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
-  }
-
-  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 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) {
-    DoeLog(2) && (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;
-  }
-
-  DoLog(1) && (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);
-    DoLog(0) && (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();
-  DoLog(1) && (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)) {
-    DoLog(1) && (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) {
-      DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
-      performCriticalExit();
-    }
-  }
-  DoLog(1) && (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);
-  DoLog(1) && (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;
-    DoLog(0) && (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) {
-      DoeLog(0) && (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++) {
+    Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
+    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+  }
+
+  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
     connectedCircle->push_back(AngleRunner->second);
   }
@@ -4006 +3013 @@
  * @return list of the all points linked to the provided one
  */
-ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
-  Info FunctionInfo(__func__);
+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;
@@ -4022 +3029 @@
   class BoundaryLineSet *CurrentLine = NULL;
   class BoundaryLineSet *StartLine = NULL;
+
   // find the respective boundary point
   PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
@@ -4027 +3035 @@
     ReferencePoint = PointRunner->second;
   } else {
-    DoeLog(1) && (eLog() << Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
+    eLog() << Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
     return NULL;
   }
 
-  map<class BoundaryLineSet *, bool> TouchedLine;
-  map<class BoundaryTriangleSet *, bool> TouchedTriangle;
-  map<class BoundaryLineSet *, bool>::iterator LineRunner;
-  map<class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
+  map <class BoundaryLineSet *, bool> TouchedLine;
+  map <class BoundaryTriangleSet *, bool> TouchedTriangle;
+  map <class BoundaryLineSet *, bool>::iterator LineRunner;
+  map <class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
   for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
-    TouchedLine.insert(pair<class BoundaryLineSet *, bool> (Runner->second, false));
+    TouchedLine.insert( pair <class BoundaryLineSet *, bool>(Runner->second, false) );
     for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++)
-      TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool> (Sprinter->second, false));
+      TouchedTriangle.insert( pair <class BoundaryTriangleSet *, bool>(Sprinter->second, false) );
   }
   if (!ReferencePoint->lines.empty()) {
@@ -4044 +3052 @@
       LineRunner = TouchedLine.find(runner->second);
       if (LineRunner == TouchedLine.end()) {
-        DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl);
+        eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl;
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new TesselPointList;
+        connectedPath = new list<TesselPoint*>;
         triangle = NULL;
         CurrentLine = runner->second;
         StartLine = CurrentLine;
         CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
-        DoLog(1) && (Log() << Verbose(1) << "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl);
+        Log() << Verbose(1)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;
         do {
           // push current one
-          DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl);
+          Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
           connectedPath->push_back(CurrentPoint->node);
 
           // find next triangle
           for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
-            DoLog(1) && (Log() << Verbose(1) << "INFO: Inspecting triangle " << *Runner->second << "." << endl);
+            Log() << Verbose(1) << "INFO: Inspecting triangle " << *Runner->second << "." << endl;
             if ((Runner->second != triangle)) { // look for first triangle not equal to old one
               triangle = Runner->second;
@@ -4067 +3075 @@
                 if (!TriangleRunner->second) {
                   TriangleRunner->second = true;
-                  DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl);
+                  Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl;
                   break;
                 } else {
-                  DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl);
+                  Log() << Verbose(1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl;
                   triangle = NULL;
                 }
               } else {
-                DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl);
+                eLog() << Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl;
                 triangle = NULL;
               }
@@ -4082 +3090 @@
             break;
           // find next line
-          for (int i = 0; i < 3; i++) {
+          for (int i=0;i<3;i++) {
             if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
               CurrentLine = triangle->lines[i];
-              DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl);
+              Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl;
               break;
             }
@@ -4091 +3099 @@
           LineRunner = TouchedLine.find(CurrentLine);
           if (LineRunner == TouchedLine.end())
-            DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl);
+            eLog() << Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl;
           else
             LineRunner->second = true;
@@ -4099 +3107 @@
         } while (CurrentLine != StartLine);
         // last point is missing, as it's on start line
-        DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl);
+        Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
         if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
           connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
@@ -4105 +3113 @@
         ListOfPaths->push_back(connectedPath);
       } else {
-        DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl);
+        Log() << Verbose(1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl;
       }
     }
   } else {
-    DoeLog(1) && (eLog() << Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl);
+    eLog() << Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl;
   }
 
@@ -4121 +3129 @@
  * @return list of the closed paths
  */
-ListOfTesselPointList * 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;
 
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Current path is " << connectedPath << "." << endl);
+    Log() << Verbose(1) << "INFO: Current path is " << connectedPath << "." << endl;
 
     // go through list, look for reappearance of starting Point and count
     CircleStart = connectedPath->begin();
+
     // 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
-        DoLog(1) && (Log() << Verbose(1) << count + 1 << ". closed path consists of: ");
-        newPath = new TesselPointList;
-        TesselPointList::iterator CircleSprinter = Marker;
+        Log() << Verbose(1) << count+1 << ". closed path consists of: ";
+        newPath = new list<TesselPoint*>;
+        list<TesselPoint*>::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
-          DoLog(0) && (Log() << Verbose(0) << (**CircleSprinter) << " <-> ");
+          Log() << Verbose(0) << (**CircleSprinter) << " <-> ";
         }
-        DoLog(0) && (Log() << Verbose(0) << ".." << endl);
+        Log() << Verbose(0) << ".." << endl;
         count++;
         Marker = CircleRunner;
@@ -4160 +3171 @@
     }
   }
-  DoLog(1) && (Log() << Verbose(1) << "INFO: " << count << " closed additional path(s) have been created." << endl);
+  Log() << Verbose(1) << "INFO: " << count << " closed additional path(s) have been created." << endl;
 
   // delete list of paths
@@ -4166 +3177 @@
     connectedPath = *(ListofPaths->begin());
     ListofPaths->remove(connectedPath);
-    delete (connectedPath);
-  }
-  delete (ListofPaths);
+    delete(connectedPath);
+  }
+  delete(ListofPaths);
 
   // exit
   return ListofClosedPaths;
-}
-;
+};
+
 
 /** Gets all belonging triangles for a given BoundaryPointSet.
@@ -4180 +3191 @@
  * \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) {
-    DoeLog(1) && (eLog() << Verbose(1) << "Point given is NULL." << endl);
+    eLog() << Verbose(1) << "Point given is NULL." << endl;
   } else {
     // go through its lines and insert all triangles
     for (LineMap::const_iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
       for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
-        connectedTriangles->insert(TriangleRunner->second);
-      }
+      connectedTriangles->insert(TriangleRunner->second);
+    }
   }
 
   return connectedTriangles;
-}
-;
+};
+
 
 /** Removes a boundary point from the envelope while keeping it closed.
@@ -4210 +3221 @@
  * \return volume added to the volume inside the tesselated surface by the removal
  */
-double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point)
-{
+double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point) {
   class BoundaryLineSet *line = NULL;
   class BoundaryTriangleSet *triangle = NULL;
@@ -4219 +3229 @@
 
   if (point == NULL) {
-    DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl);
+    eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl;
     return 0.;
   } else
-    DoLog(0) && (Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl);
+    Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl;
 
   // copy old location for the volume
@@ -4229 +3239 @@
   // get list of connected points
   if (point->lines.empty()) {
-    DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl);
+    eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl;
     return 0.;
   }
 
-  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;
+    count+=LineRunner->second->triangles.size();
+  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) );
     }
   }
 
   // remove all triangles
-  count = 0;
+  count=0;
   NormalVector.Zero();
-  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    DoLog(1) && (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++;
   }
-  DoLog(1) && (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;
+  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;
   double angle;
   double smallestangle;
   Vector Point, Reference, OrthogonalVector;
-  if (count > 2) { // less than three triangles, then nothing will be created
+  if (count > 2) {  // less than three triangles, then nothing will be created
     class TesselPoint *TriangleCandidates[3];
     count = 0;
-    for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
+    for ( ; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) {  // go through all closed paths
       if (ListAdvance != ListOfClosedPaths->end())
         ListAdvance++;
 
       connectedPath = *ListRunner;
+
       // re-create all triangles by going through connected points list
-      LineList NewLines;
-      for (; !connectedPath->empty();) {
+      list<class BoundaryLineSet *> NewLines;
+      for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
         EndNode = connectedPath->end();
         smallestangle = 0.;
         for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
-          DoLog(1) && (Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);
+          Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
           // construct vectors to next and previous neighbour
           StartNode = MiddleNode;
@@ -4302 +3313 @@
           angle = GetAngle(Point, Reference, OrthogonalVector);
           //if (angle < M_PI)  // no wrong-sided triangles, please?
-          if (fabs(angle - M_PI) < fabs(smallestangle - M_PI)) { // get straightest angle (i.e. construct those triangles with smallest area first)
-            smallestangle = angle;
-            EndNode = MiddleNode;
-          }
+            if(fabs(angle - M_PI) < fabs(smallestangle - M_PI)) {  // get straightest angle (i.e. construct those triangles with smallest area first)
+              smallestangle = angle;
+              EndNode = MiddleNode;
+            }
         }
         MiddleNode = EndNode;
         if (MiddleNode == connectedPath->end()) {
-          DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl);
+          eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl;
           performCriticalExit();
         }
@@ -4319 +3330 @@
         if (EndNode == connectedPath->end())
           EndNode = connectedPath->begin();
-        DoLog(2) && (Log() << Verbose(2) << "INFO: StartNode is " << **StartNode << "." << endl);
-        DoLog(2) && (Log() << Verbose(2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);
-        DoLog(2) && (Log() << Verbose(2) << "INFO: EndNode is " << **EndNode << "." << endl);
-        DoLog(1) && (Log() << Verbose(1) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl);
+        Log() << Verbose(2) << "INFO: StartNode is " << **StartNode << "." << endl;
+        Log() << Verbose(2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+        Log() << Verbose(2) << "INFO: EndNode is " << **EndNode << "." << endl;
+        Log() << Verbose(1) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl;
         TriangleCandidates[0] = *StartNode;
         TriangleCandidates[1] = *MiddleNode;
@@ -4328 +3339 @@
         triangle = GetPresentTriangle(TriangleCandidates);
         if (triangle != NULL) {
-          DoeLog(0) && (eLog() << Verbose(0) << "New triangle already present, skipping!" << endl);
+          eLog() << Verbose(0) << "New triangle already present, skipping!" << endl;
           StartNode++;
           MiddleNode++;
@@ -4340 +3351 @@
           continue;
         }
-        DoLog(3) && (Log() << Verbose(3) << "Adding new triangle points." << endl);
+        Log() << Verbose(3) << "Adding new triangle points."<< endl;
         AddTesselationPoint(*StartNode, 0);
         AddTesselationPoint(*MiddleNode, 1);
         AddTesselationPoint(*EndNode, 2);
-        DoLog(3) && (Log() << Verbose(3) << "Adding new triangle lines." << endl);
-        AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
-        AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+        Log() << Verbose(3) << "Adding new triangle lines."<< endl;
+        AddTesselationLine(TPS[0], TPS[1], 0);
+        AddTesselationLine(TPS[0], TPS[2], 1);
         NewLines.push_back(BLS[1]);
-        AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
+        AddTesselationLine(TPS[1], TPS[2], 2);
         BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
         BTS->GetNormalVector(NormalVector);
@@ -4359 +3370 @@
         // prepare nodes for next triangle
         StartNode = EndNode;
-        DoLog(2) && (Log() << Verbose(2) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl);
+        Log() << Verbose(2) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl;
         connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
         if (connectedPath->size() == 2) { // we are done
@@ -4366 +3377 @@
           break;
         } else if (connectedPath->size() < 2) { // something's gone wrong!
-          DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl);
+          eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl;
           performCriticalExit();
         } else {
@@ -4381 +3392 @@
       // 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) {
@@ -4395 +3406 @@
           if (maxgain != 0) {
             volume += maxgain;
-            DoLog(1) && (Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl);
+            Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl;
             OtherBase = FlipBaseline(*Candidate);
             NewLines.erase(Candidate);
@@ -4404 +3415 @@
 
       ListOfClosedPaths->remove(connectedPath);
-      delete (connectedPath);
-    }
-    DoLog(0) && (Log() << Verbose(0) << count << " triangles were created." << endl);
+      delete(connectedPath);
+    }
+    Log() << Verbose(0) << count << " triangles were created." << endl;
   } else {
     while (!ListOfClosedPaths->empty()) {
@@ -4412 +3423 @@
       connectedPath = *ListRunner;
       ListOfClosedPaths->remove(connectedPath);
-      delete (connectedPath);
-    }
-    DoLog(0) && (Log() << Verbose(0) << "No need to create any triangles." << endl);
-  }
-  delete (ListOfClosedPaths);
-
-  DoLog(0) && (Log() << Verbose(0) << "Removed volume is " << volume << "." << endl);
+      delete(connectedPath);
+    }
+    Log() << Verbose(0) << "No need to create any triangles." << endl;
+  }
+  delete(ListOfClosedPaths);
+
+  Log() << Verbose(0) << "Removed volume is " << volume << "." << endl;
 
   return volume;
-}
-;
+};
+
+
 
 /**
  * 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:
-      DoeLog(0) && (eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl);
-      performCriticalExit();
-      break;
+    }
   }
 
   return result;
 }
-
-struct BoundaryLineSetCompare
-{
-  bool operator()(const BoundaryLineSet * const a, const BoundaryLineSet * const b)
-  {
-    int lowerNra = -1;
-    int lowerNrb = -1;
-
-    if (a->endpoints[0] < a->endpoints[1])
-      lowerNra = 0;
-    else
-      lowerNra = 1;
-
-    if (b->endpoints[0] < b->endpoints[1])
-      lowerNrb = 0;
-    else
-      lowerNrb = 1;
-
-    if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])
-      return true;
-    else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
-      return false;
-    else { // both lower-numbered endpoints are the same ...
-      if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2])
-        return true;
-      else if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2])
-        return false;
-    }
-    return false;
-  }
-  ;
-};
-
-#define UniqueLines set < class BoundaryLineSet *, BoundaryLineSetCompare>
 
 /**
@@ -4560 +3493 @@
  *         in the list, once as key and once as value
  */
-IndexToIndex * Tesselation::FindAllDegeneratedLines()
-{
-  Info FunctionInfo(__func__);
-  UniqueLines AllLines;
-  IndexToIndex * DegeneratedLines = new IndexToIndex;
+map<int, int> * Tesselation::FindAllDegeneratedLines()
+{
+        Info FunctionInfo(__func__);
+  map<int, class BoundaryLineSet *> AllLines;
+  map<int, int> * DegeneratedLines = new map<int, int>;
 
   // sanity check
   if (LinesOnBoundary.empty()) {
-    DoeLog(2) && (eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.");
+    eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.";
     return DegeneratedLines;
   }
+
   LineMap::iterator LineRunner1;
-  pair<UniqueLines::iterator, bool> tester;
+  pair<LineMap::iterator, bool> tester;
   for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
-    tester = AllLines.insert(LineRunner1->second);
-    if (!tester.second) { // found degenerated line
-      DegeneratedLines->insert(pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr));
-      DegeneratedLines->insert(pair<int, int> ((*tester.first)->Nr, LineRunner1->second->Nr));
+    tester = AllLines.insert( pair<int,BoundaryLineSet *> (LineRunner1->second->endpoints[0]->Nr, LineRunner1->second) );
+    if ((!tester.second) && (tester.first->second->endpoints[1]->Nr == LineRunner1->second->endpoints[1]->Nr)) { // found degenerated line
+      DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, tester.first->second->Nr) );
+      DegeneratedLines->insert ( pair<int, int> (tester.first->second->Nr, LineRunner1->second->Nr) );
     }
   }
@@ -4583 +3517 @@
   AllLines.clear();
 
-  DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl);
-  IndexToIndex::iterator it;
-  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
-    const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
-    const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second);
-    if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end())
-      DoLog(0) && (Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl);
-    else
-      DoeLog(1) && (eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl);
-  }
+  Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
+  map<int,int>::iterator it;
+  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++)
+      Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
 
   return DegeneratedLines;
@@ -4603 +3531 @@
  *         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;
   LineMap::iterator Liner;
   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);
@@ -4622 +3551 @@
     for (TriangleRunner1 = line1->triangles.begin(); TriangleRunner1 != line1->triangles.end(); ++TriangleRunner1) {
       for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) {
-        if ((TriangleRunner1->second != TriangleRunner2->second) && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
-          DegeneratedTriangles->insert(pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr));
-          DegeneratedTriangles->insert(pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr));
+        if ((TriangleRunner1->second != TriangleRunner2->second)
+          && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr) );
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr) );
         }
       }
     }
   }
-  delete (DegeneratedLines);
-
-  DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl);
-  IndexToIndex::iterator it;
+  delete(DegeneratedLines);
+
+  Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
+  map<int,int>::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
-    DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);
+      Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
 
   return DegeneratedTriangles;
@@ -4645 +3575 @@
 void Tesselation::RemoveDegeneratedTriangles()
 {
-  Info FunctionInfo(__func__);
-  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
+        Info FunctionInfo(__func__);
+  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
-  int count = 0;
-
-  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner) {
+  int count  = 0;
+
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+    TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
+  ) {
     finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
     if (finder != TrianglesOnBoundary.end())
@@ -4668 +3600 @@
         trianglesShareLine = trianglesShareLine || triangle->lines[i] == partnerTriangle->lines[j];
 
-    if (trianglesShareLine && (triangle->endpoints[1]->LinesCount > 2) && (triangle->endpoints[2]->LinesCount > 2) && (triangle->endpoints[0]->LinesCount > 2)) {
+    if (trianglesShareLine
+      && (triangle->endpoints[1]->LinesCount > 2)
+      && (triangle->endpoints[2]->LinesCount > 2)
+      && (triangle->endpoints[0]->LinesCount > 2)
+    ) {
       // check whether we have to fix lines
       BoundaryTriangleSet *Othertriangle = NULL;
@@ -4688 +3624 @@
             // the line of triangle receives the degenerated ones
             triangle->lines[i]->triangles.erase(Othertriangle->Nr);
-            triangle->lines[i]->triangles.insert(TrianglePair(partnerTriangle->Nr, partnerTriangle));
-            for (int k = 0; k < 3; k++)
+            triangle->lines[i]->triangles.insert( TrianglePair( partnerTriangle->Nr, partnerTriangle) );
+            for (int k=0;k<3;k++)
               if (triangle->lines[i] == Othertriangle->lines[k]) {
                 Othertriangle->lines[k] = partnerTriangle->lines[j];
@@ -4695 +3631 @@
               }
             // the line of partnerTriangle receives the non-degenerated ones
-            partnerTriangle->lines[j]->triangles.erase(partnerTriangle->Nr);
-            partnerTriangle->lines[j]->triangles.insert(TrianglePair(Othertriangle->Nr, Othertriangle));
+            partnerTriangle->lines[j]->triangles.erase( partnerTriangle->Nr);
+            partnerTriangle->lines[j]->triangles.insert( TrianglePair( Othertriangle->Nr, Othertriangle) );
             partnerTriangle->lines[j] = triangle->lines[i];
           }
@@ -4702 +3638 @@
       // erase the pair
       count += (int) DegeneratedTriangles->erase(triangle->Nr);
-      DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl);
+      Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;
       RemoveTesselationTriangle(triangle);
       count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
-      DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl);
+      Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;
       RemoveTesselationTriangle(partnerTriangle);
     } else {
-      DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle << " and its partner " << *partnerTriangle << " because it is essential for at" << " least one of the endpoints to be kept in the tesselation structure." << endl);
-    }
-  }
-  delete (DegeneratedTriangles);
+      Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle
+        << " and its partner " << *partnerTriangle << " because it is essential for at"
+        << " least one of the endpoints to be kept in the tesselation structure." << endl;
+    }
+  }
+  delete(DegeneratedTriangles);
   if (count > 0)
     LastTriangle = NULL;
 
-  DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl);
+  Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl;
 }
 
@@ -4728 +3666 @@
 void Tesselation::AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   // 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;
@@ -4741 +3679 @@
     NearestBoundaryPoint = PointRunner->second;
   } else {
-    DoeLog(1) && (eLog() << Verbose(1) << "I cannot find the boundary point." << endl);
+    eLog() << Verbose(1) << "I cannot find the boundary point." << endl;
     return;
   }
-  DoLog(0) && (Log() << Verbose(0) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl);
+  Log() << Verbose(0) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl;
 
   // go through its lines and find the best one to split
@@ -4759 +3697 @@
     CenterToPoint.SubtractVector(point->node);
     angle = CenterToPoint.Angle(&BaseLine);
-    if (fabs(angle - M_PI / 2.) < fabs(BestAngle - M_PI / 2.)) {
+    if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
       BestAngle = angle;
       BestLine = Runner->second;
@@ -4769 +3707 @@
   BestLine->triangles.erase(TempTriangle->Nr);
   int nr = -1;
-  for (int i = 0; i < 3; i++) {
+  for (int i=0;i<3; i++) {
     if (TempTriangle->lines[i] == BestLine) {
       nr = i;
@@ -4777 +3715 @@
 
   // create new triangle to connect point (connects automatically with the missing spot of the chosen line)
-  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);
+  Log() << Verbose(2) << "Adding new triangle points."<< endl;
   AddTesselationPoint((BestLine->endpoints[0]->node), 0);
   AddTesselationPoint((BestLine->endpoints[1]->node), 1);
   AddTesselationPoint(point, 2);
-  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);
-  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
-  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
-  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
+  Log() << Verbose(2) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   BTS->GetNormalVector(TempTriangle->NormalVector);
   BTS->NormalVector.Scale(-1.);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl);
+  Log() << Verbose(1) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl;
   AddTesselationTriangle();
 
   // create other side of this triangle and close both new sides of the first created triangle
-  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);
+  Log() << Verbose(2) << "Adding new triangle points."<< endl;
   AddTesselationPoint((BestLine->endpoints[0]->node), 0);
   AddTesselationPoint((BestLine->endpoints[1]->node), 1);
   AddTesselationPoint(point, 2);
-  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);
-  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
-  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
-  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
+  Log() << Verbose(2) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   BTS->GetNormalVector(TempTriangle->NormalVector);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl);
+  Log() << Verbose(1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;
   AddTesselationTriangle();
 
   // add removed triangle to the last open line of the second triangle
-  for (int i = 0; i < 3; i++) { // look for the same line as BestLine (only it's its degenerated companion)
+  for (int i=0;i<3;i++) { // look for the same line as BestLine (only it's its degenerated companion)
     if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
-      if (BestLine == BTS->lines[i]) {
-        DoeLog(0) && (eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl);
+      if (BestLine == BTS->lines[i]){
+        eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl;
         performCriticalExit();
       }
-      BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle));
+      BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) );
       TempTriangle->lines[nr] = BTS->lines[i];
       break;
     }
   }
-}
-;
+};
 
 /** Writes the envelope to file.
@@ -4827 +3764 @@
 void Tesselation::Output(const char *filename, const PointCloud * const cloud)
 {
-  Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
   ofstream *tempstream = NULL;
   string NameofTempFile;
@@ -4833 +3770 @@
 
   if (LastTriangle != NULL) {
-    sprintf(NumberName, "-%04d-%s_%s_%s", (int) TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name);
+    sprintf(NumberName, "-%04d-%s_%s_%s", (int)TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name);
     if (DoTecplotOutput) {
       string NameofTempFile(filename);
       NameofTempFile.append(NumberName);
-      for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+      NameofTempFile.erase(npos, 1);
       NameofTempFile.append(TecplotSuffix);
-      DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");
+      Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
       tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
       WriteTecplotFile(tempstream, this, cloud, TriangleFilesWritten);
       tempstream->close();
       tempstream->flush();
-      delete (tempstream);
+      delete(tempstream);
     }
 
@@ -4851 +3788 @@
       string NameofTempFile(filename);
       NameofTempFile.append(NumberName);
-      for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+      NameofTempFile.erase(npos, 1);
       NameofTempFile.append(Raster3DSuffix);
-      DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");
+      Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
       tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
       WriteRaster3dFile(tempstream, this, cloud);
@@ -4860 +3797 @@
       tempstream->close();
       tempstream->flush();
-      delete (tempstream);
+      delete(tempstream);
     }
   }
   if (DoTecplotOutput || DoRaster3DOutput)
     TriangleFilesWritten++;
-}
-;
-
-struct BoundaryPolygonSetCompare
-{
-  bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const
-  {
-    if (s1->endpoints.size() < s2->endpoints.size())
-      return true;
-    else if (s1->endpoints.size() > s2->endpoints.size())
-      return false;
-    else { // equality of number of endpoints
-      PointSet::const_iterator Walker1 = s1->endpoints.begin();
-      PointSet::const_iterator Walker2 = s2->endpoints.begin();
-      while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
-        if ((*Walker1)->Nr < (*Walker2)->Nr)
-          return true;
-        else if ((*Walker1)->Nr > (*Walker2)->Nr)
-          return false;
-        Walker1++;
-        Walker2++;
-      }
-      return false;
-    }
-  }
-};
-
-#define UniquePolygonSet set < BoundaryPolygonSet *, BoundaryPolygonSetCompare>
-
-/** Finds all degenerated polygons and calls ReTesselateDegeneratedPolygon()/
- * \return number of polygons found
- */
-int Tesselation::CorrectAllDegeneratedPolygons()
-{
-  Info FunctionInfo(__func__);
-  /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
-  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
-  set<BoundaryPointSet *> EndpointCandidateList;
-  pair<set<BoundaryPointSet *>::iterator, bool> InsertionTester;
-  pair<map<int, Vector *>::iterator, bool> TriangleInsertionTester;
-  for (PointMap::const_iterator Runner = PointsOnBoundary.begin(); Runner != PointsOnBoundary.end(); Runner++) {
-    DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Runner->second << "." << endl);
-    map<int, Vector *> TriangleVectors;
-    // gather all NormalVectors
-    DoLog(1) && (Log() << Verbose(1) << "Gathering triangles ..." << endl);
-    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)
-            DoLog(1) && (Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl);
-        } else {
-          DoLog(1) && (Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl);
-        }
-      }
-    // check whether there are two that are parallel
-    DoLog(1) && (Log() << Verbose(1) << "Finding two parallel triangles ..." << endl);
-    for (map<int, Vector *>::iterator VectorWalker = TriangleVectors.begin(); VectorWalker != TriangleVectors.end(); VectorWalker++)
-      for (map<int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++)
-        if (VectorWalker != VectorRunner) { // skip equals
-          const double SCP = VectorWalker->second->ScalarProduct(VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles
-          DoLog(1) && (Log() << Verbose(1) << "Checking " << *VectorWalker->second << " against " << *VectorRunner->second << ": " << SCP << endl);
-          if (fabs(SCP + 1.) < ParallelEpsilon) {
-            InsertionTester = EndpointCandidateList.insert((Runner->second));
-            if (InsertionTester.second)
-              DoLog(0) && (Log() << Verbose(0) << " Adding " << *Runner->second << " to endpoint candidate list." << endl);
-            // and break out of both loops
-            VectorWalker = TriangleVectors.end();
-            VectorRunner = TriangleVectors.end();
-            break;
-          }
-        }
-  }
-  delete (DegeneratedTriangles);
-  /// 3. Find connected endpoint candidates and put them into a polygon
-  UniquePolygonSet ListofDegeneratedPolygons;
-  BoundaryPointSet *Walker = NULL;
-  BoundaryPointSet *OtherWalker = NULL;
-  BoundaryPolygonSet *Current = NULL;
-  stack<BoundaryPointSet*> ToCheckConnecteds;
-  while (!EndpointCandidateList.empty()) {
-    Walker = *(EndpointCandidateList.begin());
-    if (Current == NULL) { // create a new polygon with current candidate
-      DoLog(0) && (Log() << Verbose(0) << "Starting new polygon set at point " << *Walker << endl);
-      Current = new BoundaryPolygonSet;
-      Current->endpoints.insert(Walker);
-      EndpointCandidateList.erase(Walker);
-      ToCheckConnecteds.push(Walker);
-    }
-
-    // go through to-check stack
-    while (!ToCheckConnecteds.empty()) {
-      Walker = ToCheckConnecteds.top(); // fetch ...
-      ToCheckConnecteds.pop(); // ... and remove
-      for (LineMap::const_iterator LineWalker = Walker->lines.begin(); LineWalker != Walker->lines.end(); LineWalker++) {
-        OtherWalker = (LineWalker->second)->GetOtherEndpoint(Walker);
-        DoLog(1) && (Log() << Verbose(1) << "Checking " << *OtherWalker << endl);
-        set<BoundaryPointSet *>::iterator Finder = EndpointCandidateList.find(OtherWalker);
-        if (Finder != EndpointCandidateList.end()) { // found a connected partner
-          DoLog(1) && (Log() << Verbose(1) << " Adding to polygon." << endl);
-          Current->endpoints.insert(OtherWalker);
-          EndpointCandidateList.erase(Finder); // remove from candidates
-          ToCheckConnecteds.push(OtherWalker); // but check its partners too
-        } else {
-          DoLog(1) && (Log() << Verbose(1) << " is not connected to " << *Walker << endl);
-        }
-      }
-    }
-
-    DoLog(0) && (Log() << Verbose(0) << "Final polygon is " << *Current << endl);
-    ListofDegeneratedPolygons.insert(Current);
-    Current = NULL;
-  }
-
-  const int counter = ListofDegeneratedPolygons.size();
-
-  DoLog(0) && (Log() << Verbose(0) << "The following " << counter << " degenerated polygons have been found: " << endl);
-  for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++)
-    DoLog(0) && (Log() << Verbose(0) << " " << **PolygonRunner << endl);
-
-  /// 4. Go through all these degenerated polygons
-  for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++) {
-    stack<int> TriangleNrs;
-    Vector NormalVector;
-    /// 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) {
-      DoLog(1) && (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) {
-      DoeLog(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
-    /// 4a. Get NormalVector for one side (this is "front")
-    NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
-    DoLog(1) && (Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl);
-    TriangleWalker++;
-    TriangleSet::iterator TriangleSprinter = TriangleWalker; // is the inner advanced iterator
-    /// 4b. Remove all triangles whose NormalVector is in opposite direction (i.e. "back")
-    BoundaryTriangleSet *triangle = NULL;
-    while (TriangleSprinter != T->end()) {
-      TriangleWalker = TriangleSprinter;
-      triangle = *TriangleWalker;
-      TriangleSprinter++;
-      DoLog(1) && (Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl);
-      if (triangle->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list
-        DoLog(1) && (Log() << Verbose(1) << " Removing ... " << endl);
-        TriangleNrs.push(triangle->Nr);
-        T->erase(TriangleWalker);
-        RemoveTesselationTriangle(triangle);
-      } else
-        DoLog(1) && (Log() << Verbose(1) << " Keeping ... " << endl);
-    }
-    /// 4c. Copy all "front" triangles but with inverse NormalVector
-    TriangleWalker = T->begin();
-    while (TriangleWalker != T->end()) { // go through all front triangles
-      DoLog(1) && (Log() << Verbose(1) << " Re-creating triangle " << **TriangleWalker << " with NormalVector " << (*TriangleWalker)->NormalVector << endl);
-      for (int i = 0; i < 3; i++)
-        AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i);
-      AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
-      AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
-      AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
-      if (TriangleNrs.empty())
-        DoeLog(0) && (eLog() << Verbose(0) << "No more free triangle numbers!" << endl);
-      BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ...
-      AddTesselationTriangle(); // ... and add
-      TriangleNrs.pop();
-      BTS->NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
-      BTS->NormalVector.Scale(-1.);
-      TriangleWalker++;
-    }
-    if (!TriangleNrs.empty()) {
-      DoeLog(0) && (eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl);
-    }
-    delete (T); // remove the triangleset
-  }
-  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
-  DoLog(0) && (Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl);
-  IndexToIndex::iterator it;
-  for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
-    DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);
-  delete (SimplyDegeneratedTriangles);
-  /// 5. exit
-  UniquePolygonSet::iterator PolygonRunner;
-  while (!ListofDegeneratedPolygons.empty()) {
-    PolygonRunner = ListofDegeneratedPolygons.begin();
-    delete (*PolygonRunner);
-    ListofDegeneratedPolygons.erase(PolygonRunner);
-  }
-
-  return counter;
-}
-;
+};