Changeset b8d1aeb for src/tesselation.cpp

Timestamp:

Feb 2, 2010, 11:38:06 AM (16 years ago)

Author:

Tillmann Crueger <crueger@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

7ba324

Parents:

9fe36b (diff), 2ededc2 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'MenuRefactoring' into QT4Refactoring

Conflicts:

molecuilder/src/Makefile.am
molecuilder/src/builder.cpp
molecuilder/src/unittests/Makefile.am

File:

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

src/tesselation.cpp

@@ -9 +9 @@
 
 #include "helpers.hpp"
+#include "info.hpp"
 #include "linkedcell.hpp"
 #include "log.hpp"
@@ -22 +23 @@
 /** Constructor of BoundaryPointSet.
  */
-BoundaryPointSet::BoundaryPointSet()
-{
-  LinesCount = 0;
-  Nr = -1;
-  value = 0.;
+BoundaryPointSet::BoundaryPointSet() :
+    LinesCount(0),
+    value(0.),
+    Nr(-1)
+{
+        Info FunctionInfo(__func__);
+        Log() << Verbose(1) << "Adding noname." << endl;
 };
 
@@ -32 +35 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker)
-{
-  node = Walker;
-  LinesCount = 0;
-  Nr = Walker->nr;
-  value = 0.;
+BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
+  LinesCount(0),
+  node(Walker),
+  value(0.),
+  Nr(Walker->nr)
+{
+        Info FunctionInfo(__func__);
+  Log() << Verbose(1) << "Adding Node " << *Walker << endl;
 };
 
@@ -46 +51 @@
 BoundaryPointSet::~BoundaryPointSet()
 {
-  //Log() << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+        Info FunctionInfo(__func__);
+  //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl;
   if (!lines.empty())
     eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl;
@@ -55 +61 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
-{
-  Log() << Verbose(6) << "Adding " << *this << " to line " << *line << "."
+void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
+{
+        Info FunctionInfo(__func__);
+  Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "."
       << endl;
   if (line->endpoints[0] == this)
@@ -85 +92 @@
 /** Constructor of BoundaryLineSet.
  */
-BoundaryLineSet::BoundaryLineSet()
-{
+BoundaryLineSet::BoundaryLineSet() :
+    Nr(-1)
+{
+        Info FunctionInfo(__func__);
   for (int i = 0; i < 2; i++)
     endpoints[i] = NULL;
-  Nr = -1;
 };
 
@@ -97 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
-{
+BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
+{
+        Info FunctionInfo(__func__);
   // set number
   Nr = number;
@@ -106 +115 @@
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
+  // set skipped to false
+  skipped = false;
   // clear triangles list
-  Log() << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+  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
+  Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
 };
 
@@ -116 +149 @@
 BoundaryLineSet::~BoundaryLineSet()
 {
+        Info FunctionInfo(__func__);
   int Numbers[2];
 
@@ -134 +168 @@
         for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
           if ((*Runner).second == this) {
-            //Log() << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+            //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
             endpoints[i]->lines.erase(Runner);
             break;
@@ -140 +174 @@
       } else { // there's just a single line left
         if (endpoints[i]->lines.erase(Nr)) {
-          //Log() << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+          //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
         }
       }
       if (endpoints[i]->lines.empty()) {
-        //Log() << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+        //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         if (endpoints[i] != NULL) {
           delete(endpoints[i]);
@@ -159 +193 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
-{
-  Log() << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
+void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
+{
+        Info FunctionInfo(__func__);
+  Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
   triangles.insert(TrianglePair(triangle->Nr, triangle));
 };
@@ -169 +204 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
-{
+bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
+{
+        Info FunctionInfo(__func__);
   if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
     return true;
@@ -183 +219 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion()
-{
+bool BoundaryLineSet::CheckConvexityCriterion() const
+{
+        Info FunctionInfo(__func__);
   Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
   // get the two triangles
   if (triangles.size() != 2) {
-    eLog() << Verbose(1) << "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;
   }
   // check normal vectors
   // have a normal vector on the base line pointing outwards
-  //Log() << Verbose(3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
+  //Log() << Verbose(0) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
   BaseLineCenter.CopyVector(endpoints[0]->node->node);
   BaseLineCenter.AddVector(endpoints[1]->node->node);
@@ -199 +236 @@
   BaseLine.CopyVector(endpoints[0]->node->node);
   BaseLine.SubtractVector(endpoints[1]->node->node);
-  //Log() << Verbose(3) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
+  //Log() << Verbose(0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
 
   BaseLineNormal.Zero();
@@ -206 +243 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
-    //Log() << Verbose(3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
+  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+    //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     NormalCheck.AddVector(&runner->second->NormalVector);
     NormalCheck.Scale(sign);
@@ -214 +251 @@
       BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
     else {
-      Log() << Verbose(1) << "CRITICAL: Triangle " << *runner->second << " has zero normal vector!" << endl;
-      exit(255);
+      eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl;
     }
     node = runner->second->GetThirdEndpoint(this);
     if (node != NULL) {
-      //Log() << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
+      //Log() << Verbose(0) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
       helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
-      //Log() << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
+      //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
     } else {
-      //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;
     }
   }
-  //Log() << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
+  //Log() << Verbose(0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
   if (NormalCheck.NormSquared() < MYEPSILON) {
-    Log() << Verbose(3) << "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;
   }
@@ -238 +274 @@
   double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
   if ((angle - M_PI) > -MYEPSILON) {
-    Log() << Verbose(3) << "ACCEPT: Angle is greater than pi: convex." << endl;
+    Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl;
     return true;
   } else {
-    Log() << Verbose(3) << "REJECT: Angle is less than pi: concave." << endl;
+    Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl;
     return false;
   }
@@ -250 +286 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
-{
+bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+{
+        Info FunctionInfo(__func__);
   for(int i=0;i<2;i++)
     if (point == endpoints[i])
@@ -262 +299 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
-{
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
+{
+        Info FunctionInfo(__func__);
   if (endpoints[0] == point)
     return endpoints[1];
@@ -286 +324 @@
 /** Constructor for BoundaryTriangleSet.
  */
-BoundaryTriangleSet::BoundaryTriangleSet()
-{
+BoundaryTriangleSet::BoundaryTriangleSet() :
+  Nr(-1)
+{
+        Info FunctionInfo(__func__);
   for (int i = 0; i < 3; i++)
     {
@@ -293 +333 @@
       lines[i] = NULL;
     }
-  Nr = -1;
 };
 
@@ -300 +339 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
-{
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
+  Nr(number)
+{
+        Info FunctionInfo(__func__);
   // set number
-  Nr = number;
   // set lines
-  Log() << Verbose(5) << "New triangle " << Nr << ":" << endl;
-  for (int i = 0; i < 3; i++)
-    {
-      lines[i] = line[i];
-      lines[i]->AddTriangle(this);
-    }
+  for (int i = 0; i < 3; i++) {
+    lines[i] = line[i];
+    lines[i]->AddTriangle(this);
+  }
   // get ascending order of endpoints
-  map<int, class BoundaryPointSet *> OrderMap;
+  PointMap OrderMap;
   for (int i = 0; i < 3; i++)
     // 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]));
-        // and we don't care whether insertion fails
-      }
+    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]));
+      // and we don't care whether insertion fails
+    }
   // set endpoints
   int Counter = 0;
-  Log() << Verbose(6) << " with end points ";
-  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
-      != OrderMap.end(); runner++)
-    {
-      endpoints[Counter] = runner->second;
-      Log() << Verbose(0) << " " << *endpoints[Counter];
-      Counter++;
-    }
-  if (Counter < 3)
-    {
-      eLog() << Verbose(0) << "ERROR! We have a triangle with only two distinct endpoints!" << endl;
-      performCriticalExit();
-    }
-  Log() << Verbose(0) << "." << 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;
+    Log() << Verbose(0) << " " << *endpoints[Counter] << endl;
+    Counter++;
+  }
+  if (Counter < 3) {
+    eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl;
+    performCriticalExit();
+  }
 };
 
@@ -345 +378 @@
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
+        Info FunctionInfo(__func__);
   for (int i = 0; i < 3; i++) {
     if (lines[i] != NULL) {
       if (lines[i]->triangles.erase(Nr)) {
-        //Log() << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+        //Log() << Verbose(0) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
       }
       if (lines[i]->triangles.empty()) {
-          //Log() << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+          //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
           delete (lines[i]);
           lines[i] = NULL;
@@ -357 +391 @@
     }
   }
-  //Log() << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+  //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl;
 };
 
@@ -364 +398 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
-{
+void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
+{
+        Info FunctionInfo(__func__);
   // get normal vector
   NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);
@@ -372 +407 @@
   if (NormalVector.ScalarProduct(&OtherVector) > 0.)
     NormalVector.Scale(-1.);
-};
-
-/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+  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.
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
  * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
@@ -386 +422 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
-{
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
+{
+  Info FunctionInfo(__func__);
   Vector CrossPoint;
   Vector helper;
 
   if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) {
-    Log() << 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;
   }
 
+  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
+  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
+
+  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl;
+    return true;
+  }
   // Calculate cross point between one baseline and the line from the third endpoint to intersection
   int i=0;
@@ -402 +453 @@
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
+      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+      const double s = CrossPoint.ScalarProduct(&helper)/helper.NormSquared();
+      Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl;
+      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
+        Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl;
+        i=4;
+        break;
+      }
+      i++;
     } else 
-      i++;
-    if (i>2)
       break;
-  } while (CrossPoint.NormSquared() < MYEPSILON);
+  } while (i<3);
   if (i==3) {
-    eLog() << Verbose(1) << "Could not find any cross points, something's utterly wrong here!" << endl;
-    exit(255);
-  }
-  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-
-  // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
     return true;
-  } else { // outside!
-    Intersection->Zero();
+  } else {
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
     return false;
   }
+};
+
+/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
+ * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+ * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
+ * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
+ * the first two basepoints) or not.
+ * \param *x point
+ * \param *ClosestPoint desired closest point inside triangle to \a *x, is absolute vector
+ * \return Distance squared between \a *x and closest point inside triangle
+ */
+double BoundaryTriangleSet::GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const
+{
+  Info FunctionInfo(__func__);
+  Vector Direction;
+
+  // 1. get intersection with plane
+  Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
+  GetCenter(&Direction);
+  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
+    ClosestPoint->CopyVector(x);
+  }
+
+  // 2. Calculate in plane part of line (x, intersection)
+  Vector InPlane;
+  InPlane.CopyVector(x);
+  InPlane.SubtractVector(ClosestPoint);  // points from plane intersection to straight-down point
+  InPlane.ProjectOntoPlane(&NormalVector);
+  InPlane.AddVector(ClosestPoint);
+
+  Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl;
+  Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl;
+
+  // Calculate cross point between one baseline and the desired point such that distance is shortest
+  double ShortestDistance = -1.;
+  bool InsideFlag = false;
+  Vector CrossDirection[3];
+  Vector CrossPoint[3];
+  Vector helper;
+  for (int i=0;i<3;i++) {
+    // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
+    Direction.CopyVector(endpoints[(i+1)%3]->node->node);
+    Direction.SubtractVector(endpoints[i%3]->node->node);
+    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
+    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
+    CrossDirection[i].CopyVector(&CrossPoint[i]);
+    CrossDirection[i].SubtractVector(&InPlane);
+    CrossPoint[i].SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+    const double s = CrossPoint[i].ScalarProduct(&Direction)/Direction.NormSquared();
+    Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl;
+    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
+      CrossPoint[i].AddVector(endpoints[i%3]->node->node);  // make cross point absolute again
+      Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i%3]->node->node << " and " << *endpoints[(i+1)%3]->node->node << "." << endl;
+      const double distance = CrossPoint[i].DistanceSquared(x);
+      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
+        ShortestDistance = distance;
+        ClosestPoint->CopyVector(&CrossPoint[i]);
+      }
+    } else
+      CrossPoint[i].Zero();
+  }
+  InsideFlag = true;
+  for (int i=0;i<3;i++) {
+    const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+1)%3]);
+    const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+2)%3]);;
+    if ((sign > -MYEPSILON) && (othersign > -MYEPSILON))  // have different sign
+      InsideFlag = false;
+  }
+  if (InsideFlag) {
+    ClosestPoint->CopyVector(&InPlane);
+    ShortestDistance = InPlane.DistanceSquared(x);
+  } else {  // also check endnodes
+    for (int i=0;i<3;i++) {
+      const double distance = x->DistanceSquared(endpoints[i]->node->node);
+      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
+        ShortestDistance = distance;
+        ClosestPoint->CopyVector(endpoints[i]->node->node);
+      }
+    }
+  }
+  Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl;
+  return ShortestDistance;
 };
 
@@ -426 +562 @@
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
-{
+bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+{
+        Info FunctionInfo(__func__);
   for(int i=0;i<3;i++)
     if (line == lines[i])
@@ -438 +575 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
-{
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+{
+        Info FunctionInfo(__func__);
   for(int i=0;i<3;i++)
     if (point == endpoints[i])
@@ -450 +588 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
-{
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+{
+        Info FunctionInfo(__func__);
   for(int i=0;i<3;i++)
     if (point == endpoints[i]->node)
@@ -462 +601 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
-{
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+        Log() << Verbose(1) << "INFO: Checking " << Points[0] << ","  << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl;
   return (((endpoints[0] == Points[0])
             || (endpoints[0] == Points[1])
@@ -483 +624 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
-{
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
+{
+        Info FunctionInfo(__func__);
   return (((endpoints[0] == T->endpoints[0])
             || (endpoints[0] == T->endpoints[1])
@@ -504 +646 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
-{
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
+{
+        Info FunctionInfo(__func__);
   // sanity check
   if (!ContainsBoundaryLine(line))
@@ -520 +663 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector *center)
-{
+void BoundaryTriangleSet::GetCenter(Vector * const center) const
+{
+        Info FunctionInfo(__func__);
   center->Zero();
   for(int i=0;i<3;i++)
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
+  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -534 +679 @@
 ostream &operator <<(ostream &ost, const BoundaryTriangleSet &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 << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
+  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 << "]";
   return ost;
 };
 
+// ======================================== Polygons on Boundary =================================
+
+/** Constructor for BoundaryPolygonSet.
+ */
+BoundaryPolygonSet::BoundaryPolygonSet() :
+  Nr(-1)
+{
+  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();
+  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()) {
+        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.);
+  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);
+  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++) {
+    Log() << Verbose(0) << "Checking against " << **Runner << endl;
+    if (point == (*Runner)) {
+      Log() << Verbose(0) << " Contained." << endl;
+      return true;
+    }
+  }
+  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) {
+      Log() << Verbose(0) << " Contained." << endl;
+      return true;
+    }
+  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;
+  Log() << Verbose(1) << "Polygon is " << *this << endl;
+  for(int i=0;i<dim;i++) {
+    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;
+  Log() << Verbose(1) << "Polygon is " << *this << endl;
+  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+    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)
+            Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl;
+        }
+      }
+
+  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;
+  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)
+        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 ===========================================
 
@@ -545 +939 @@
 TesselPoint::TesselPoint()
 {
+  //Info FunctionInfo(__func__);
   node = NULL;
   nr = -1;
@@ -554 +949 @@
 TesselPoint::~TesselPoint()
 {
+  //Info FunctionInfo(__func__);
 };
 
@@ -568 +964 @@
 ostream & TesselPoint::operator << (ostream &ost)
 {
-  ost << "[" << (Name) << "|" << this << "]";
+        Info FunctionInfo(__func__);
+  ost << "[" << (nr) << "|" << this << "]";
   return ost;
 };
@@ -579 +976 @@
 PointCloud::PointCloud()
 {
-
+        //Info FunctionInfo(__func__);
 };
 
@@ -586 +983 @@
 PointCloud::~PointCloud()
 {
-
+        //Info FunctionInfo(__func__);
 };
 
@@ -593 +990 @@
 /** Constructor of class CandidateForTesselation.
  */
-CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) {
-  point = candidate;
-  BaseLine = line;
+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, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
+    BaseLine(line),
+    ShortestAngle(2.*M_PI),
+    OtherShortestAngle(2.*M_PI)
+{
+        Info FunctionInfo(__func__);
   OptCenter.CopyVector(&OptCandidateCenter);
   OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
@@ -603 +1014 @@
  */
 CandidateForTesselation::~CandidateForTesselation() {
-  point = NULL;
   BaseLine = NULL;
 };
 
+/** output operator for CandidateForTesselation.
+ * \param &ost output stream
+ * \param &a boundary line
+ */
+ostream & operator <<(ostream &ost, const  CandidateForTesselation &a)
+{
+  ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->Name << "," << a.BaseLine->endpoints[1]->node->Name << "] with ";
+  if (a.pointlist.empty())
+    ost << "no candidate.";
+  else {
+    ost << "candidate";
+    if (a.pointlist.size() != 1)
+      ost << "s ";
+    else
+      ost << " ";
+    for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++)
+      ost << *(*Runner) << " ";
+    ost << " at angle " << (a.ShortestAngle)<< ".";
+  }
+
+  return ost;
+};
+
+
 // =========================================================== class TESSELATION ===========================================
 
 /** Constructor of class Tesselation.
  */
-Tesselation::Tesselation()
-{
-  PointsOnBoundaryCount = 0;
-  LinesOnBoundaryCount = 0;
-  TrianglesOnBoundaryCount = 0;
-  InternalPointer = PointsOnBoundary.begin();
-  LastTriangle = NULL;
-  TriangleFilesWritten = 0;
+Tesselation::Tesselation() :
+  PointsOnBoundaryCount(0),
+  LinesOnBoundaryCount(0),
+  TrianglesOnBoundaryCount(0),
+  LastTriangle(NULL),
+  TriangleFilesWritten(0),
+  InternalPointer(PointsOnBoundary.begin())
+{
+        Info FunctionInfo(__func__);
 }
 ;
@@ -627 +1062 @@
 Tesselation::~Tesselation()
 {
-  Log() << Verbose(1) << "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) {
@@ -635 +1071 @@
       eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl;
   }
-  Log() << Verbose(1) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
+  Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
 }
 ;
@@ -644 +1080 @@
 Vector * Tesselation::GetCenter(ofstream *out) const
 {
+        Info FunctionInfo(__func__);
   Vector *Center = new Vector(0.,0.,0.);
   int num=0;
@@ -659 +1096 @@
 TesselPoint * Tesselation::GetPoint() const
 {
+        Info FunctionInfo(__func__);
   return (InternalPointer->second->node);
 };
@@ -667 +1105 @@
 TesselPoint * Tesselation::GetTerminalPoint() const
 {
+        Info FunctionInfo(__func__);
   PointMap::const_iterator Runner = PointsOnBoundary.end();
   Runner--;
@@ -677 +1116 @@
 void Tesselation::GoToNext() const
 {
+        Info FunctionInfo(__func__);
   if (InternalPointer != PointsOnBoundary.end())
     InternalPointer++;
@@ -686 +1126 @@
 void Tesselation::GoToPrevious() const
 {
+        Info FunctionInfo(__func__);
   if (InternalPointer != PointsOnBoundary.begin())
     InternalPointer--;
@@ -695 +1136 @@
 void Tesselation::GoToFirst() const
 {
+        Info FunctionInfo(__func__);
   InternalPointer = PointsOnBoundary.begin();
 };
@@ -703 +1145 @@
 void Tesselation::GoToLast() const
 {
+        Info FunctionInfo(__func__);
   InternalPointer = PointsOnBoundary.end();
   InternalPointer--;
@@ -712 +1155 @@
 bool Tesselation::IsEmpty() const
 {
+        Info FunctionInfo(__func__);
   return (PointsOnBoundary.empty());
 };
@@ -720 +1164 @@
 bool Tesselation::IsEnd() const
 {
+        Info FunctionInfo(__func__);
   return (InternalPointer == PointsOnBoundary.end());
 };
@@ -729 +1174 @@
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
  */
-void
-Tesselation::GuessStartingTriangle()
-{
+void Tesselation::GuessStartingTriangle()
+{
+        Info FunctionInfo(__func__);
   // 4b. create a starting triangle
   // 4b1. create all distances
@@ -778 +1223 @@
           baseline->second.first->second->node->node,
           baseline->second.second->second->node->node);
-      Log() << Verbose(2) << "Plane vector of candidate triangle is ";
-      PlaneVector.Output();
-      Log() << Verbose(0) << endl;
+      Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl;
       // 4. loop over all points
       double sign = 0.;
@@ -796 +1239 @@
           if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
             continue;
-          Log() << Verbose(3) << "Projection of " << checker->second->node->Name
-              << " yields distance of " << distance << "." << endl;
+          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)
@@ -850 +1292 @@
       if (checker == PointsOnBoundary.end())
         {
-          Log() << Verbose(0) << "Looks like we have a candidate!" << endl;
+          Log() << Verbose(2) << "Looks like we have a candidate!" << endl;
           break;
         }
@@ -880 +1322 @@
   else
     {
-      Log() << Verbose(1) << "No starting triangle found." << endl;
-      exit(255);
+      eLog() << Verbose(0) << "No starting triangle found." << endl;
     }
 }
@@ -901 +1342 @@
 void Tesselation::TesselateOnBoundary(const PointCloud * const cloud)
 {
+        Info FunctionInfo(__func__);
   bool flag;
   PointMap::iterator winner;
@@ -919 +1361 @@
         // 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
-        Log() << Verbose(2) << "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];
-        Log() << Verbose(3) << " and has peak " << *peak << "." << endl;
+        Log() << Verbose(1) << " and has peak " << *peak << "." << endl;
 
         // prepare some auxiliary vectors
@@ -938 +1380 @@
           CenterVector.AddVector(BTS->endpoints[i]->node->node);
         CenterVector.Scale(1. / 3.);
-        Log() << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl;
+        Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl;
 
         // normal vector of triangle
@@ -945 +1387 @@
         BTS->GetNormalVector(NormalVector);
         NormalVector.CopyVector(&BTS->NormalVector);
-        Log() << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
+        Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl;
 
         // vector in propagation direction (out of triangle)
@@ -952 +1394 @@
         TempVector.CopyVector(&CenterVector);
         TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
-        //Log() << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
+        //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
         if (PropagationVector.ScalarProduct(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
           PropagationVector.Scale(-1.);
-        Log() << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl;
+        Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl;
         winner = PointsOnBoundary.end();
 
@@ -961 +1403 @@
         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
-            Log() << Verbose(3) << "Target point is " << *(target->second) << ":" << endl;
+            Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl;
 
             // first check direction, so that triangles don't intersect
@@ -968 +1410 @@
             VirtualNormalVector.ProjectOntoPlane(&NormalVector);
             TempAngle = VirtualNormalVector.Angle(&PropagationVector);
-            Log() << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << 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(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
+              Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
               continue;
             } else
-              Log() << Verbose(4) << "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)
@@ -979 +1421 @@
             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))) {
-              Log() << Verbose(4) << *(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))) {
-              Log() << Verbose(4) << *(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;
             }
@@ -1000 +1442 @@
             helper.ProjectOntoPlane(&TempVector);
             if (fabs(helper.NormSquared()) < MYEPSILON) {
-              Log() << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
+              Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl;
               continue;
             }
@@ -1017 +1459 @@
             // calculate angle
             TempAngle = NormalVector.Angle(&VirtualNormalVector);
-            Log() << Verbose(4) << "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;
-              Log() << Verbose(4) << "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...
@@ -1039 +1481 @@
                 SmallestAngle = TempAngle;
                 winner = target;
-                Log() << Verbose(4) << "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
-                Log() << Verbose(4) << "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
-              Log() << Verbose(4) << "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
@@ -1049 +1491 @@
         // 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()) {
-          Log() << Verbose(2) << "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;
@@ -1079 +1521 @@
           TrianglesOnBoundaryCount++;
         } else {
-          Log() << Verbose(1) << "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
-        Log() << Verbose(2) << "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);
 
@@ -1099 +1541 @@
 bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC)
 {
+        Info FunctionInfo(__func__);
   Vector Intersection, Normal;
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
-
-  Log() << Verbose(1) << "Begin of InsertStraddlingPoints" << endl;
 
   cloud->GoToFirst();
@@ -1117 +1558 @@
     }
     Walker = cloud->GetPoint();
-    Log() << Verbose(2) << "Current point is " << *Walker << "." << endl;
+    Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
-      Log() << Verbose(2) << "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 {
     }
-    Log() << Verbose(2) << "Closest triangle is " << *BTS << "." << endl;
+    Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl;
     // get the intersection point
     if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) {
-      Log() << Verbose(2) << "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!
-        Log() << Verbose(2) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;
+        Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;
       } else {
         // outside!
-        Log() << Verbose(2) << *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;
@@ -1147 +1588 @@
         Normal.CopyVector(&BTS->NormalVector);
         // add Walker to boundary points
-        Log() << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
+        Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl;
         AddFlag = true;
         if (AddBoundaryPoint(Walker,0))
@@ -1154 +1595 @@
           continue;
         // remove triangle
-        Log() << Verbose(2) << "Erasing triangle " << *BTS << "." << endl;
+        Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl;
         TrianglesOnBoundary.erase(BTS->Nr);
         delete(BTS);
@@ -1162 +1603 @@
           BPS[1] = OldPoints[i];
           NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-          Log() << Verbose(3) << "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++;
@@ -1173 +1614 @@
             if (NewLines[j]->IsConnectedTo(BLS[0])) {
               if (n>2) {
-                Log() << Verbose(1) << BLS[0] << " connects to all of the new lines?!" << endl;
+                eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl;
                 return false;
               } else
@@ -1184 +1625 @@
           BTS->GetNormalVector(Normal);
           Normal.Scale(-1.);
-          Log() << Verbose(2) << "Created new triangle " << *BTS << "." << endl;
+          Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl;
           TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
           TrianglesOnBoundaryCount++;
@@ -1198 +1639 @@
   // exit
   delete(Center);
-  Log() << Verbose(1) << "End of InsertStraddlingPoints" << endl;
   return true;
 };
@@ -1209 +1649 @@
 bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n)
 {
+        Info FunctionInfo(__func__);
   PointTestPair InsertUnique;
   BPS[n] = new class BoundaryPointSet(Walker);
@@ -1230 +1671 @@
 void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n)
 {
+        Info FunctionInfo(__func__);
   PointTestPair InsertUnique;
   TPS[n] = new class BoundaryPointSet(Candidate);
@@ -1237 +1679 @@
   } else {
     delete TPS[n];
-    Log() << Verbose(4) << "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;
   }
@@ -1250 +1692 @@
 void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const
 {
+        Info FunctionInfo(__func__);
   PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr);
   if (FindPoint != PointsOnBoundary.end())
@@ -1267 +1710 @@
   bool insertNewLine = true;
 
-  if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  if (FindLine != a->lines.end()) {
+    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
+
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-    Log() << Verbose(5) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
 
     for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
@@ -1277 +1721 @@
       if (FindLine->second->triangles.size() < 2) {
         insertNewLine = false;
-        Log() << Verbose(4) << "Using existing line " << *FindLine->second << endl;
+        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]);
+        if (CandidateLine != OpenLines.end()) {
+          Log() << Verbose(1) << " Removing line from OpenLines." << endl;
+          delete(CandidateLine->second);
+          OpenLines.erase(CandidateLine);
+        } else {
+          eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl;
+        }
 
         break;
@@ -1304 +1758 @@
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
 {
-  Log() << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
+        Info FunctionInfo(__func__);
+  Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
@@ -1312 +1767 @@
   // increase counter
   LinesOnBoundaryCount++;
+  // also add to open lines
+  CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
+  OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
 };
 
@@ -1319 +1777 @@
 void Tesselation::AddTesselationTriangle()
 {
+        Info FunctionInfo(__func__);
   Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
 
@@ -1337 +1796 @@
 void Tesselation::AddTesselationTriangle(const int nr)
 {
-  Log() << Verbose(1) << "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
@@ -1355 +1815 @@
 void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle)
 {
+        Info FunctionInfo(__func__);
   if (triangle == NULL)
     return;
   for (int i = 0; i < 3; i++) {
     if (triangle->lines[i] != NULL) {
-      Log() << Verbose(5) << "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()) {
-          Log() << Verbose(5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
+          Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
           RemoveTesselationLine(triangle->lines[i]);
       } else {
-        Log() << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle: ";
+        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++)
           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(5) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[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";
@@ -1382 +1844 @@
 
   if (TrianglesOnBoundary.erase(triangle->Nr))
-    Log() << Verbose(5) << "Removing triangle Nr. " << triangle->Nr << "." << endl;
+    Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl;
   delete(triangle);
 };
@@ -1392 +1854 @@
 void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line)
 {
+        Info FunctionInfo(__func__);
   int Numbers[2];
 
@@ -1412 +1875 @@
         for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
           if ((*Runner).second == line) {
-            Log() << Verbose(5) << "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;
@@ -1418 +1881 @@
       } else { // there's just a single line left
         if (line->endpoints[i]->lines.erase(line->Nr))
-          Log() << Verbose(5) << "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()) {
-        Log() << Verbose(5) << *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 {
-        Log() << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to: ";
+        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";
@@ -1437 +1900 @@
 
   if (LinesOnBoundary.erase(line->Nr))
-    Log() << Verbose(5) << "Removing line Nr. " << line->Nr << "." << endl;
+    Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl;
   delete(line);
 };
@@ -1448 +1911 @@
 void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point)
 {
+        Info FunctionInfo(__func__);
   if (point == NULL)
     return;
   if (PointsOnBoundary.erase(point->Nr))
-    Log() << Verbose(5) << "Removing point Nr. " << point->Nr << "." << endl;
+    Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl;
   delete(point);
 };
@@ -1466 +1930 @@
 int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const
 {
+        Info FunctionInfo(__func__);
   int adjacentTriangleCount = 0;
   class BoundaryPointSet *Points[3];
 
-  Log() << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
   // builds a triangle point set (Points) of the end points
   for (int i = 0; i < 3; i++) {
@@ -1488 +1952 @@
           for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
             TriangleMap *triangles = &FindLine->second->triangles;
-            Log() << Verbose(3) << "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)) {
@@ -1494 +1958 @@
               }
             }
-            Log() << Verbose(3) << "end." << endl;
+            Log() << Verbose(1) << "end." << endl;
           }
           // Only one of the triangle lines must be considered for the triangle count.
-          //Log() << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
           //return adjacentTriangleCount;
         }
@@ -1504 +1968 @@
   }
 
-  Log() << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-  Log() << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+  Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
   return adjacentTriangleCount;
 };
@@ -1519 +1982 @@
 class BoundaryTriangleSet * Tesselation::GetPresentTriangle(TesselPoint *Candidates[3])
 {
+        Info FunctionInfo(__func__);
   class BoundaryTriangleSet *triangle = NULL;
   class BoundaryPointSet *Points[3];
@@ -1548 +2012 @@
           }
           // Only one of the triangle lines must be considered for the triangle count.
-          //Log() << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
           //return adjacentTriangleCount;
         }
@@ -1569 +2033 @@
 void Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC)
 {
-  Log() << Verbose(1) << "Begin of FindStartingTriangle\n";
+        Info FunctionInfo(__func__);
   int i = 0;
-  TesselPoint* FirstPoint = NULL;
-  TesselPoint* SecondPoint = NULL;
   TesselPoint* MaxPoint[NDIM];
+  TesselPoint* Temporary;
   double maxCoordinate[NDIM];
-  Vector Oben;
+  BoundaryLineSet BaseLine;
   Vector helper;
   Vector Chord;
   Vector SearchDirection;
-
-  Oben.Zero();
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector SphereCenter;
+  Vector NormalVector;
+
+  NormalVector.Zero();
 
   for (i = 0; i < 3; i++) {
@@ -1593 +2060 @@
       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(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) {
             if ((*Runner)->node->x[i] > maxCoordinate[i]) {
-              Log() << Verbose(2) << "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);
@@ -1608 +2075 @@
   }
 
-  Log() << Verbose(2) << "Found maximum coordinates: ";
+  Log() << Verbose(1) << "Found maximum coordinates: ";
   for (int i=0;i<NDIM;i++)
     Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t";
@@ -1614 +2081 @@
 
   BTS = NULL;
-  CandidateList *OptCandidates = new CandidateList();
   for (int k=0;k<NDIM;k++) {
-    Oben.Zero();
-    Oben.x[k] = 1.;
-    FirstPoint = MaxPoint[k];
-    Log() << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;
+    NormalVector.Zero();
+    NormalVector.x[k] = 1.;
+    BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
+    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
 
     double ShortestAngle;
-    TesselPoint* OptCandidate = NULL;
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
 
-    FindSecondPointForTesselation(FirstPoint, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
-    SecondPoint = OptCandidate;
-    if (SecondPoint == NULL)  // have we found a second point?
+    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    if (Temporary == NULL)  // have we found a second point?
       continue;
-
-    helper.CopyVector(FirstPoint->node);
-    helper.SubtractVector(SecondPoint->node);
-    helper.Normalize();
-    Oben.ProjectOntoPlane(&helper);
-    Oben.Normalize();
-    helper.VectorProduct(&Oben);
+    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.
 
-    Chord.CopyVector(FirstPoint->node); // bring into calling function
-    Chord.SubtractVector(SecondPoint->node);
-    double radius = Chord.ScalarProduct(&Chord);
-    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
-    helper.CopyVector(&Oben);
-    helper.Scale(CircleRadius);
-    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+    SphereCenter.CopyVector(&NormalVector);
+    SphereCenter.Scale(CircleRadius);
+    SphereCenter.AddVector(&CircleCenter);
+    // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
-    Log() << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;
-    AddTesselationPoint(FirstPoint, 0);
-    Log() << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
-    AddTesselationPoint(SecondPoint, 1);
-    AddTesselationLine(TPS[0], TPS[1], 0);
-
-    //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    FindThirdPointForTesselation(Oben, SearchDirection, helper, BLS[0], NULL, *&OptCandidates, &ShortestAngle, RADIUS, LC);
-    Log() << Verbose(1) << "List of third Points is ";
-    for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
-        Log() << Verbose(0) << " " << *(*it)->point;
-    }
-    Log() << Verbose(0) << endl;
-
-    for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
-      // add third triangle point
-      AddTesselationPoint((*it)->point, 2);
-      // add the second and third line
-      AddTesselationLine(TPS[1], TPS[2], 1);
-      AddTesselationLine(TPS[0], TPS[2], 2);
-      // ... and triangles to the Maps of the Tesselation class
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTesselationTriangle();
-      // ... and calculate its normal vector (with correct orientation)
-      (*it)->OptCenter.Scale(-1.);
-      Log() << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
-      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
-      Log() << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
-      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
-
-      // if we do not reach the end with the next step of iteration, we need to setup a new first line
-      if (it != OptCandidates->end()--) {
-        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
-        SecondPoint = (*it)->point;
-        // adding point 1 and point 2 and the line between them
-        AddTesselationPoint(FirstPoint, 0);
-        AddTesselationPoint(SecondPoint, 1);
-        AddTesselationLine(TPS[0], TPS[1], 0);
-      }
-      Log() << Verbose(2) << "Projection is " << BTS->NormalVector.ScalarProduct(&Oben) << "." << endl;
-    }
+    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);
+    Log() << Verbose(0) << "List of third Points is:" << endl;
+    for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) {
+        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 {
       // 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);
-      }
-      OptCandidates->clear();
-    }
-  }
-
-  // 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);
-  Log() << Verbose(1) << "End of FindStartingTriangle\n";
+      OptCandidates.pointlist.clear();
+    }
+  }
 };
 
 /** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates.
  * This is supposed to prevent early closing of the tesselation.
- * \param *BaseRay baseline, i.e. not \a *OptCandidate
+ * \param CandidateLine CandidateForTesselation with baseline and shortestangle , i.e. not \a *OptCandidate
  * \param *ThirdNode third point in triangle, not in BoundaryLineSet::endpoints
- * \param ShortestAngle path length on this circle band for the current \a *ThirdNode
  * \param RADIUS radius of sphere
  * \param *LC LinkedCell structure
  * \return true - there is a better candidate (smaller angle than \a ShortestAngle), false - no better TesselPoint candidate found
  */
-bool Tesselation::HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const ThirdNode, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const
-{
-  bool result = false;
-  Vector CircleCenter;
-  Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
-  Vector SearchDirection;
-  Vector helper;
-  TesselPoint *OtherOptCandidate = NULL;
-  double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-  double radius, CircleRadius;
-  BoundaryLineSet *Line = NULL;
-  BoundaryTriangleSet *T = NULL;
-
-  Log() << Verbose(1) << "Begin of HasOtherBaselineBetterCandidate" << endl;
-
-  // check both other lines
-  PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr);
-  if (FindPoint != PointsOnBoundary.end()) {
-    for (int i=0;i<2;i++) {
-      LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr);
-      if (FindLine != (FindPoint->second)->lines.end()) {
-        Line = FindLine->second;
-        Log() << Verbose(1) << "Found line " << *Line << "." << endl;
-        if (Line->triangles.size() == 1) {
-          T = Line->triangles.begin()->second;
-          // construct center of circle
-          CircleCenter.CopyVector(Line->endpoints[0]->node->node);
-          CircleCenter.AddVector(Line->endpoints[1]->node->node);
-          CircleCenter.Scale(0.5);
-
-          // construct normal vector of circle
-          CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node);
-          CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node);
-
-          // calculate squared radius of circle
-          radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-          if (radius/4. < RADIUS*RADIUS) {
-            CircleRadius = RADIUS*RADIUS - radius/4.;
-            CirclePlaneNormal.Normalize();
-            //Log() << Verbose(2) << "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 = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
-            helper.Scale(sqrt(RADIUS*RADIUS - radius));
-            OldSphereCenter.AddVector(&helper);
-            OldSphereCenter.SubtractVector(&CircleCenter);
-            //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-
-            // construct SearchDirection
-            SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal);
-            helper.CopyVector(Line->endpoints[0]->node->node);
-            helper.SubtractVector(ThirdNode->node);
-            if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
-              SearchDirection.Scale(-1.);
-            SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-            SearchDirection.Normalize();
-            Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-            if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
-              // rotated the wrong way!
-              eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
-            }
-
-            // add third point
-            CandidateList *OptCandidates = new CandidateList();
-            FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, Line, ThirdNode, OptCandidates, &OtherShortestAngle, RADIUS, LC);
-            for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
-              if (((*it)->point == BaseRay->endpoints[0]->node) || ((*it)->point == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested
-                continue;
-              Log() << Verbose(1) << " Third point candidate is " << *(*it)->point
-              << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
-              Log() << Verbose(1) << " 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;
-              bool check=false;
-              int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates);
-              // If there is no triangle, add it regularly.
-              if (existentTrianglesCount == 0) {
-                SetTesselationPoint((*it)->point, 0);
-                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
-                SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
-
-                if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
-                  OtherOptCandidate = (*it)->point;
-                  check = true;
-                }
-              } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
-                SetTesselationPoint((*it)->point, 0);
-                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
-                SetTesselationPoint(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)) {
-                  OtherOptCandidate = (*it)->point;
-                  check = true;
-                }
-              }
-
-              if (check) {
-                if (ShortestAngle > OtherShortestAngle) {
-                  Log() << Verbose(1) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl;
-                  result = true;
-                  break;
-                }
-              }
-            }
-            delete(OptCandidates);
-            if (result)
-              break;
-          } else {
-            Log() << Verbose(1) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl;
-          }
-        } else {
-          eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl;
-        }
-      } else {
-        Log() << Verbose(2) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl;
-      }
-    }
-  } else {
-    eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl;
-  }
-
-  Log() << Verbose(1) << "End of HasOtherBaselineBetterCandidate" << endl;
-
-  return result;
-};
+//bool Tesselation::HasOtherBaselineBetterCandidate(CandidateForTesselation &CandidateLine, const TesselPoint * const ThirdNode, double RADIUS, const LinkedCell * const LC) const
+//{
+//      Info FunctionInfo(__func__);
+//  bool result = false;
+//  Vector CircleCenter;
+//  Vector CirclePlaneNormal;
+//  Vector OldSphereCenter;
+//  Vector SearchDirection;
+//  Vector helper;
+//  TesselPoint *OtherOptCandidate = NULL;
+//  double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+//  double radius, CircleRadius;
+//  BoundaryLineSet *Line = NULL;
+//  BoundaryTriangleSet *T = NULL;
+//
+//  // check both other lines
+//  PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr);
+//  if (FindPoint != PointsOnBoundary.end()) {
+//    for (int i=0;i<2;i++) {
+//      LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr);
+//      if (FindLine != (FindPoint->second)->lines.end()) {
+//        Line = FindLine->second;
+//        Log() << Verbose(0) << "Found line " << *Line << "." << endl;
+//        if (Line->triangles.size() == 1) {
+//          T = Line->triangles.begin()->second;
+//          // construct center of circle
+//          CircleCenter.CopyVector(Line->endpoints[0]->node->node);
+//          CircleCenter.AddVector(Line->endpoints[1]->node->node);
+//          CircleCenter.Scale(0.5);
+//
+//          // construct normal vector of circle
+//          CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node);
+//          CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node);
+//
+//          // calculate squared radius of circle
+//          radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+//          if (radius/4. < RADIUS*RADIUS) {
+//            CircleRadius = RADIUS*RADIUS - radius/4.;
+//            CirclePlaneNormal.Normalize();
+//            //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+//
+//            // 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 = Line->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(Line->endpoints[0]->node->node);
+//            helper.SubtractVector(ThirdNode->node);
+//            if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
+//              SearchDirection.Scale(-1.);
+//            SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+//            SearchDirection.Normalize();
+//            Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+//            if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+//              // rotated the wrong way!
+//              eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+//            }
+//
+//            // add third point
+//            FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, OptCandidates, ThirdNode, RADIUS, LC);
+//            for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); ++it) {
+//              if (((*it) == BaseRay->endpoints[0]->node) || ((*it) == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested
+//                continue;
+//              Log() << Verbose(0) << " Third point candidate is " << (*it)
+//              << " 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);
+//              PointCandidates[1] = BaseRay->endpoints[0]->node;
+//              PointCandidates[2] = BaseRay->endpoints[1]->node;
+//              bool check=false;
+//              int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates);
+//              // If there is no triangle, add it regularly.
+//              if (existentTrianglesCount == 0) {
+//                SetTesselationPoint((*it), 0);
+//                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
+//                SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
+//
+//                if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
+//                  OtherOptCandidate = (*it);
+//                  check = true;
+//                }
+//              } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
+//                SetTesselationPoint((*it), 0);
+//                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
+//                SetTesselationPoint(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)) {
+//                  OtherOptCandidate = (*it);
+//                  check = true;
+//                }
+//              }
+//
+//              if (check) {
+//                if (ShortestAngle > OtherShortestAngle) {
+//                  Log() << Verbose(0) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl;
+//                  result = true;
+//                  break;
+//                }
+//              }
+//            }
+//            delete(OptCandidates);
+//            if (result)
+//              break;
+//          } else {
+//            Log() << Verbose(0) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl;
+//          }
+//        } else {
+//          eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl;
+//        }
+//      } else {
+//        Log() << Verbose(1) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl;
+//      }
+//    }
+//  } else {
+//    eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl;
+//  }
+//
+//  return result;
+//};
 
 /** This function finds a triangle to a line, adjacent to an existing one.
  * @param out output stream for debugging
- * @param Line current baseline to search from
+ * @param CandidateLine current cadndiate baseline to search from
  * @param T current triangle which \a Line is edge of
  * @param RADIUS radius of the rolling ball
@@ -1866 +2292 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
-{
-  Log() << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
+bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
+{
+        Info FunctionInfo(__func__);
   bool result = true;
-  CandidateList *OptCandidates = new CandidateList();
 
   Vector CircleCenter;
   Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
+  Vector RelativeSphereCenter;
   Vector SearchDirection;
   Vector helper;
   TesselPoint *ThirdNode = NULL;
   LineMap::iterator testline;
-  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
   double radius, CircleRadius;
 
-  Log() << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
   for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
+    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) {
       ThirdNode = T.endpoints[i]->node;
+      break;
+    }
+  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdNode " << *ThirdNode << " of triangle " << T << "." << endl;
 
   // construct center of circle
-  CircleCenter.CopyVector(Line.endpoints[0]->node->node);
-  CircleCenter.AddVector(Line.endpoints[1]->node->node);
+  CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+  CircleCenter.AddVector(CandidateLine.BaseLine->endpoints[1]->node->node);
   CircleCenter.Scale(0.5);
 
   // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(Line.endpoints[0]->node->node);
-  CirclePlaneNormal.SubtractVector(Line.endpoints[1]->node->node);
+  CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+  CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // 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.;
     CirclePlaneNormal.Normalize();
-    //Log() << Verbose(2) << "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 = Line.endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
-    helper.Scale(sqrt(RADIUS*RADIUS - radius));
-    OldSphereCenter.AddVector(&helper);
-    OldSphereCenter.SubtractVector(&CircleCenter);
-    //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-
-    // construct SearchDirection
-    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(Line.endpoints[0]->node->node);
+    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+
+    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl;
+
+    // construct SearchDirection and an "outward pointer"
+    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
+    helper.CopyVector(&CircleCenter);
     helper.SubtractVector(ThirdNode->node);
     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
-    Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
@@ -1927 +2348 @@
 
     // add third point
-    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, OptCandidates, &ShortestAngle, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
 
   } else {
-    Log() << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
-  }
-
-  if (OptCandidates->begin() == OptCandidates->end()) {
+    Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl;
+  }
+
+  if (CandidateLine.pointlist.empty()) {
     eLog() << Verbose(2) << "Could not find a suitable candidate." << endl;
     return false;
   }
-  Log() << Verbose(1) << "Third Points are ";
-  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
-    Log() << Verbose(1) << " " << *(*it)->point << endl;
-  }
-
-  BoundaryLineSet *BaseRay = &Line;
-  for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {
-    Log() << Verbose(1) << " Third point candidate is " << *(*it)->point
-    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
-    Log() << Verbose(1) << " 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) ) {
-      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)) {
-          AddTesselationLine(TPS[0], TPS[1], 0);
-          AddTesselationLine(TPS[0], TPS[2], 1);
-          AddTesselationLine(TPS[1], TPS[2], 2);
-
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          AddTesselationTriangle();
-          (*it)->OptCenter.Scale(-1.);
-          BTS->GetNormalVector((*it)->OptCenter);
-          (*it)->OptCenter.Scale(-1.);
-
-          Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-            << " for this triangle ... " << endl;
-        //Log() << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
-        } else {
-          eLog() << Verbose(2) << "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 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)) {
-            AddTesselationLine(TPS[0], TPS[1], 0);
-            AddTesselationLine(TPS[0], TPS[2], 1);
-            AddTesselationLine(TPS[1], TPS[2], 2);
-
-            BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-            AddTesselationTriangle();  // add to global map
-
-            (*it)->OtherOptCenter.Scale(-1.);
-            BTS->GetNormalVector((*it)->OtherOptCenter);
-            (*it)->OtherOptCenter.Scale(-1.);
-
-            eLog() << Verbose(2) << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-            Log() << Verbose(1) << "We have "<< BaseRay->triangles.size() << " for line " << BaseRay << "." << endl;
-          } else {
-            eLog() << Verbose(2) << "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);
-  Log() << Verbose(0) << "End of FindNextSuitableTriangle\n";
+  Log() << Verbose(0) << "Third Points are: " << endl;
+  for (TesselPointList::iterator it = CandidateLine.pointlist.begin(); it != CandidateLine.pointlist.end(); ++it) {
+    Log() << Verbose(0) << " " << *(*it) << endl;
+  }
+
+  return true;
+
+//  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
+ * \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;
+
+  // fill the set of neighbours
+  TesselPointSet SetOfNeighbours;
+  SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
+  for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
+    SetOfNeighbours.insert(*Runner);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
+
+  // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
+  Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl;
+  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
+    Log() << Verbose(0) << **TesselRunner << endl;
+  TesselPointList::iterator Runner = connectedClosestPoints->begin();
+  TesselPointList::iterator Sprinter = Runner;
+  Sprinter++;
+  while(Sprinter != connectedClosestPoints->end()) {
+    // add the points
+    AddTesselationPoint(TurningPoint, 0);
+    AddTesselationPoint((*Runner), 1);
+    AddTesselationPoint((*Sprinter), 2);
+
+    // 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();
+    BTS->GetCenter(&Center);
+    Center.SubtractVector(&CandidateLine.OptCenter);
+    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+    BTS->GetNormalVector(Center);
+
+    Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << "." << endl;
+    Runner = Sprinter;
+    Sprinter++;
+    Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl;
+    if (Sprinter != connectedClosestPoints->end())
+      Log() << Verbose(0) << " There are still more triangles to add." << endl;
+  }
+  delete(connectedClosestPoints);
 };
 
@@ -2053 +2512 @@
 class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base)
 {
+        Info FunctionInfo(__func__);
   class BoundaryPointSet *Spot = NULL;
   class BoundaryLineSet *OtherBase;
@@ -2064 +2524 @@
   OtherBase = new class BoundaryLineSet(BPS,-1);
 
-  Log() << Verbose(3) << "INFO: Current base line is " << *Base << "." << endl;
-  Log() << Verbose(3) << "INFO: Other base line is " << *OtherBase << "." << endl;
+  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
@@ -2085 +2545 @@
   delete(ClosestPoint);
   if ((distance[0] * distance[1]) > 0)  { // have same sign?
-    Log() << Verbose(3) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1]  << ". " << *Base << "' rectangle is concave." << endl;
+    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];
@@ -2093 +2553 @@
     return Spot;
   } else {  // different sign, i.e. we are in between
-    Log() << Verbose(3) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;
+    Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;
     return NULL;
   }
@@ -2101 +2561 @@
 void Tesselation::PrintAllBoundaryPoints(ofstream *out) const
 {
+        Info FunctionInfo(__func__);
   // print all lines
-  Log() << Verbose(1) << "Printing all boundary points for debugging:" << endl;
+  Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl;
   for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++)
-    Log() << Verbose(2) << *(PointRunner->second) << endl;
+    Log() << Verbose(0) << *(PointRunner->second) << endl;
 };
 
 void Tesselation::PrintAllBoundaryLines(ofstream *out) const
 {
+        Info FunctionInfo(__func__);
   // print all lines
-  Log() << Verbose(1) << "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++)
-    Log() << Verbose(2) << *(LineRunner->second) << endl;
+    Log() << Verbose(0) << *(LineRunner->second) << endl;
 };
 
 void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const
 {
+        Info FunctionInfo(__func__);
   // print all triangles
-  Log() << Verbose(1) << "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++)
-    Log() << Verbose(2) << *(TriangleRunner->second) << endl;
+    Log() << Verbose(0) << *(TriangleRunner->second) << endl;
 };
 
@@ -2130 +2593 @@
 double Tesselation::PickFarthestofTwoBaselines(class BoundaryLineSet *Base)
 {
+        Info FunctionInfo(__func__);
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint[2];
@@ -2141 +2605 @@
   OtherBase = new class BoundaryLineSet(BPS,-1);
 
-  Log() << Verbose(3) << "INFO: Current base line is " << *Base << "." << endl;
-  Log() << Verbose(3) << "INFO: Other base line is " << *OtherBase << "." << endl;
+  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
@@ -2162 +2626 @@
 
   if (Distance.NormSquared() < MYEPSILON) { // check for intersection
-    Log() << Verbose(3) << "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
@@ -2172 +2636 @@
     }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-      Log() << Verbose(4) << "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));
     }
@@ -2178 +2642 @@
 
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
-      Log() << Verbose(2) << "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
-      Log() << Verbose(2) << "REJECT: Base line is higher: Nothing to do." << endl;
+      Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl;
       return 0.;
     }
@@ -2197 +2661 @@
 class BoundaryLineSet * Tesselation::FlipBaseline(class BoundaryLineSet *Base)
 {
+        Info FunctionInfo(__func__);
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
@@ -2203 +2668 @@
   int i,m;
 
-  Log() << Verbose(1) << "Begin of FlipBaseline" << endl;
-
   // calculate NormalVector for later use
   BaseLineNormal.Zero();
@@ -2212 +2675 @@
   }
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-    Log() << Verbose(4) << "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));
   }
@@ -2225 +2688 @@
   i=0;
   m=0;
-  Log() << Verbose(3) << "The four old lines are: ";
+  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
@@ -2233 +2696 @@
       }
   Log() << Verbose(0) << endl;
-  Log() << Verbose(3) << "The two old points are: ";
+  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
@@ -2259 +2722 @@
 
   // remove triangles and baseline removes itself
-  Log() << Verbose(3) << "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++) {
-    Log() << Verbose(3) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;
+    Log() << Verbose(0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;
     OldTriangleNrs[m++] = runner->second->Nr;
     RemoveTesselationTriangle(runner->second);
@@ -2273 +2736 @@
   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
-  Log() << Verbose(3) << "INFO: Created new baseline " << *NewLine << "." << endl;
+  Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl;
 
   // construct new triangles with flipped baseline
@@ -2288 +2751 @@
     BTS->GetNormalVector(BaseLineNormal);
     AddTesselationTriangle(OldTriangleNrs[0]);
-    Log() << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
+    Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl;
 
     BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
@@ -2296 +2759 @@
     BTS->GetNormalVector(BaseLineNormal);
     AddTesselationTriangle(OldTriangleNrs[1]);
-    Log() << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
+    Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl;
   } else {
-    Log() << Verbose(1) << "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;
   }
 
-  Log() << Verbose(1) << "End of FlipBaseline" << endl;
   return NewLine;
 };
@@ -2317 +2779 @@
 void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC)
 {
-  Log() << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl;
+        Info FunctionInfo(__func__);
   Vector AngleCheck;
   class TesselPoint* Candidate = NULL;
@@ -2338 +2800 @@
     Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
   }
-  Log() << Verbose(3) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :"
+  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;
 
@@ -2345 +2807 @@
       for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
         const LinkedNodes *List = LC->GetCurrentCell();
-        //Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
@@ -2376 +2838 @@
                 angle = AngleCheck.Angle(&Oben);
                 if (angle < Storage[0]) {
-                  //Log() << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                  Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
+                  //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+                  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;
-                  //Log() << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+                  //Log() << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
                 } else {
-                  //Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl;
+                  //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl;
                 }
               } else {
-                //Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
+                //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
               }
             } else {
-              //Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
+              //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
             }
           }
         } else {
-          Log() << Verbose(3) << "Linked cell list is empty." << endl;
+          Log() << Verbose(0) << "Linked cell list is empty." << endl;
         }
       }
-  Log() << Verbose(2) << "End of FindSecondPointForTesselation" << endl;
 };
 
@@ -2424 +2885 @@
  * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
  * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
- * @param BaseLine BoundaryLineSet with the current base line
+ * @param CandidateLine CandidateForTesselation with the current base line and list of candidates and ShortestAngle
  * @param ThirdNode third point to avoid in search
- * @param candidates list of equally good candidates to return
- * @param ShortestAngle the current path length on this circle band for the current OptCandidate
  * @param RADIUS radius of sphere
  * @param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, class BoundaryLineSet *BaseLine, const class TesselPoint  * const ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC) const
-{
+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
@@ -2440 +2900 @@
   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;
-  CandidateForTesselation *optCandidate = NULL;
-
-  Log() << Verbose(1) << "Begin of FindThirdPointForTesselation" << endl;
-
-  Log() << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+
+  Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
 
   // construct center of circle
-  CircleCenter.CopyVector(BaseLine->endpoints[0]->node->node);
-  CircleCenter.AddVector(BaseLine->endpoints[1]->node->node);
+  CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+  CircleCenter.AddVector(CandidateLine.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);
+  CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+  CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
+
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
 
   // calculate squared radius TesselPoint *ThirdNode,f circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
+  radius = CirclePlaneNormal.NormSquared()/4.;
+  if (radius < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius;
     CirclePlaneNormal.Normalize();
-    //Log() << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
     // test whether old center is on the band's plane
-    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      eLog() << Verbose(1) << "Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
-      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
-    }
-    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+      Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl;
 
       // check SearchDirection
-      //Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
         eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
       }
@@ -2486 +2949 @@
         for(int i=0;i<NDIM;i++) // store indices of this cell
         N[i] = LC->n[i];
-        //Log() << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+        //Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
       } else {
         eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
@@ -2492 +2955 @@
       }
       // then go through the current and all neighbouring cells and check the contained points for possible candidates
-      //Log() << Verbose(2) << "LC Intervals:";
+      //Log() << Verbose(1) << "LC Intervals:";
       for (int i=0;i<NDIM;i++) {
         Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
@@ -2503 +2966 @@
           for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
             const LinkedNodes *List = LC->GetCurrentCell();
-            //Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+            //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
             if (List != NULL) {
               for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
@@ -2509 +2972 @@
 
                 // check for three unique points
-                //Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->node << "." << endl;
-                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){
-
-                  // construct both new centers
-                  GetCenterofCircumcircle(&NewSphereCenter, *BaseLine->endpoints[0]->node->node, *BaseLine->endpoints[1]->node->node, *Candidate->node);
-                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
-
-                  if ((NewNormalVector.MakeNormalVector(BaseLine->endpoints[0]->node->node, BaseLine->endpoints[1]->node->node, Candidate->node))
-                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                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);
+                  Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;
+
+                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)
+                  && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
                   ) {
-                    helper.CopyVector(&NewNormalVector);
-                    //Log() << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
+                    Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
+                    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+                    Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+                    Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl;
                     if (radius < RADIUS*RADIUS) {
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
+                      if (fabs(radius - otherradius) > HULLEPSILON) {
+                        eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl;
+                      }
+                      // construct both new centers
+                      NewSphereCenter.CopyVector(&NewPlaneCenter);
+                      OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
+                      helper.CopyVector(&NewNormalVector);
                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
-                      //Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
+                      Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl;
                       NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
-                      //Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-
+                      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;
+                      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 (*ShortestAngle > (alpha - HULLEPSILON)) {
-                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
+                      if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
                         if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+                          CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
+                          CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
                         } else {
-                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
+                          CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
+                          CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
                         }
                         // if there is an equal candidate, add it to the list without clearing the list
-                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
-                          candidates->push_back(optCandidate);
-                          Log() << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
+                        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
-                          class CandidateForTesselation *remover = NULL;
-                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
-                            remover = *it;
-                            delete(remover);
-                          }
-                          candidates->clear();
-                          candidates->push_back(optCandidate);
-                          Log() << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
+                          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;
                         }
-                        *ShortestAngle = alpha;
-                        //Log() << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
+                        CandidateLine.ShortestAngle = alpha;
+                        Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl;
                       } else {
-                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
-                          //Log() << Verbose(2) << "REJECT: Old candidate " << *(optCandidate->point) << " with " << *ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << 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(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                          Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
                         }
                       }
-
                     } else {
-                      //Log() << Verbose(2) << "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 {
-                    //Log() << Verbose(2) << "REJECT: Three points from " << *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 (ThirdNode != NULL) {
-                    //Log() << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
+                    Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
                   } else {
-                    //Log() << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                    Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl;
                   }
                 }
@@ -2597 +3061 @@
   } else {
     if (ThirdNode != NULL)
-      Log() << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
+      Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
     else
-      Log() << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
-  }
-
-  //Log() << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
-  if (candidates->size() > 1) {
-    candidates->unique();
-    candidates->sort(SortCandidates);
-  }
-
-  Log() << Verbose(1) << "End of FindThirdPointForTesselation" << 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);
+  }
 };
 
@@ -2618 +3080 @@
 class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const
 {
+        Info FunctionInfo(__func__);
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  map<int, class BoundaryPointSet *> OrderMap;
-  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  PointMap OrderMap;
+  PointTestPair OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
@@ -2631 +3094 @@
           { // if insertion fails, we have common endpoint
             node = OrderTest.first->second;
-            Log() << Verbose(5) << "Common endpoint of lines " << *line1
+            Log() << Verbose(1) << "Common endpoint of lines " << *line1
                 << " and " << *line2 << " is: " << *node << "." << endl;
             j = 2;
@@ -2641 +3104 @@
 };
 
+/** Finds the boundary points that are closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
+ */
+DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+  PointMap::const_iterator FindPoint;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  if (LinesOnBoundary.empty()) {
+    eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl;
+    return NULL;
+  }
+
+  // gather all points close to the desired one
+  LC->SetIndexToVector(x); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  DistanceToPointMap * points = new DistanceToPointMap;
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //Log() << Verbose(1) << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        const LinkedNodes *List = LC->GetCurrentCell();
+        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            FindPoint = PointsOnBoundary.find((*Runner)->nr);
+            if (FindPoint != PointsOnBoundary.end()) {
+              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
+            }
+          }
+        } else {
+          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  // check whether we found some points
+  if (points->empty()) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    delete(points);
+    return NULL;
+  }
+  return points;
+};
+
+/** Finds the boundary line that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return closest BoundaryLineSet or NULL in degenerate case.
+ */
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+
+  // get closest points
+  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    return NULL;
+  }
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl;
+  BoundaryLineSet *ClosestLine = NULL;
+  double MinDistance = -1.;
+  Vector helper;
+  Vector Center;
+  Vector BaseLine;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+      // calculate closest point on line to desired point
+      helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
+      helper.Scale(0.5);
+      Center.CopyVector(x);
+      Center.SubtractVector(&helper);
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      Center.ProjectOntoPlane(&BaseLine);
+      const double distance = Center.NormSquared();
+      if ((ClosestLine == NULL) || (distance < MinDistance)) {
+        // additionally calculate intersection on line (whether it's on the line section or not)
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthA = helper.ScalarProduct(&BaseLine);
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthB = helper.ScalarProduct(&BaseLine);
+        if (lengthB*lengthA < 0) {  // if have different sign
+          ClosestLine = LineRunner->second;
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
+        } else {
+          Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
+        }
+      } else {
+        Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
+      }
+    }
+  }
+  delete(points);
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLine == NULL) {
+    Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
+    return NULL;
+  }
+  return ClosestLine;
+};
+
+
 /** Finds the triangle that is closest to a given Vector \a *x.
  * \param *out output stream for debugging
  * \param *x Vector to look from
- * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
- */
-list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const
-{
-  TesselPoint *trianglePoints[3];
-  TesselPoint *SecondPoint = NULL;
-  list<BoundaryTriangleSet*> *triangles = NULL;
-
-  if (LinesOnBoundary.empty()) {
-    Log() << Verbose(0) << "Error: There is no tesselation structure to compare the point with, please create one first.";
+ * \return BoundaryTriangleSet of nearest triangle or NULL.
+ */
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+{
+        Info FunctionInfo(__func__);
+
+        // get closest points
+        DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
     return NULL;
   }
-  Log() << Verbose(1) << "Finding closest Tesselpoint to " << *x << " ... " << endl;
-  trianglePoints[0] = FindClosestPoint(x, SecondPoint, LC);
-  
-  // check whether closest point is "too close" :), then it's inside
-  if (trianglePoints[0] == NULL) {
-    Log() << Verbose(2) << "Is the only point, no one else is closeby." << endl;
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
+  LineSet ClosestLines;
+  double MinDistance = 1e+16;
+  Vector BaseLineIntersection;
+  Vector Center;
+  Vector BaseLine;
+  Vector BaseLineCenter;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthBase = BaseLine.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+      const double lengthEndA = BaseLineIntersection.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthEndB = BaseLineIntersection.NormSquared();
+
+      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
+        const double lengthEnd = Min(lengthEndA, lengthEndB);
+        if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
+          ClosestLines.clear();
+          ClosestLines.insert(LineRunner->second);
+          MinDistance = lengthEnd;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl;
+        } else if  (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
+          ClosestLines.insert(LineRunner->second);
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl;
+        } else { // line is worse
+          Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
+        }
+      } else { // intersection is closer, calculate
+        // calculate closest point on line to desired point
+        BaseLineIntersection.CopyVector(x);
+        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        Center.CopyVector(&BaseLineIntersection);
+        Center.ProjectOntoPlane(&BaseLine);
+        BaseLineIntersection.SubtractVector(&Center);
+        const double distance = BaseLineIntersection.NormSquared();
+        if (Center.NormSquared() > BaseLine.NormSquared()) {
+          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
+        }
+        if ((ClosestLines.empty()) || (distance < MinDistance)) {
+          ClosestLines.insert(LineRunner->second);
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl;
+        } else {
+          Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
+        }
+      }
+    }
+  }
+  delete(points);
+
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLines.empty()) {
+    Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
-  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
-    Log() << Verbose(3) << "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 {
-    list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(trianglePoints[0], x);
-    if (connectedClosestPoints != NULL) {
-      trianglePoints[1] = connectedClosestPoints->front();
-      trianglePoints[2] = connectedClosestPoints->back();
-      for (int i=0;i<3;i++) {
-        if (trianglePoints[i] == NULL) {
-          eLog() << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
-        }
-        //Log() << Verbose(2) << "List of triangle points:" << endl;
-        //Log() << Verbose(3) << *trianglePoints[i] << endl;
-      }
-
-      triangles = FindTriangles(trianglePoints);
-      Log() << Verbose(2) << "List of possible triangles:" << endl;
-      for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-        Log() << Verbose(3) << **Runner << endl;
-
-      delete(connectedClosestPoints);
-    } else {
-      triangles = NULL;
-      Log() << Verbose(1) << "There is no circle of connected points!" << endl;
-    }
-  }
-  
-  if ((triangles == NULL) || (triangles->empty())) {
-    eLog() << Verbose(1) << "There is no nearest triangle. Please check the tesselation structure.";
-    delete(triangles);
-    return NULL;
-  } else
-    return triangles;
+  TriangleList * candidates = new TriangleList;
+  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
+    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
+    candidates->push_back(Runner->second);
+  }
+  return candidates;
 };
 
@@ -2727 +3317 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
-{
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
+{
+        Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
-  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
+  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
+  TriangleList candidates;
   Vector Center;
-
-  if (triangles == NULL)
+  Vector helper;
+
+  if ((triangles == NULL) || (triangles->empty()))
     return NULL;
 
-  if (triangles->size() == 1) { // there is no degenerate case
-    result = triangles->front();
-    Log() << Verbose(2) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
-  } else {
-    result = triangles->front();
-    result->GetCenter(&Center);
-    Center.SubtractVector(x);
-    Log() << Verbose(2) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
-    if (Center.ScalarProduct(&result->NormalVector) < 0) {
-      result = triangles->back();
-      Log() << Verbose(2) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
-      if (Center.ScalarProduct(&result->NormalVector) < 0) {
-        eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
-      }
+  // go through all and pick the one with the best alignment to x
+  double MinAlignment = 2.*M_PI;
+  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
+    (*Runner)->GetCenter(&Center);
+    helper.CopyVector(x);
+    helper.SubtractVector(&Center);
+    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
+    if (Alignment < MinAlignment) {
+      result = *Runner;
+      MinAlignment = Alignment;
+      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
+    } else {
+      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
     }
   }
   delete(triangles);
+
   return result;
 };
 
-/** Checks whether the provided Vector is within the tesselation structure.
+/** Checks whether the provided Vector is within the Tesselation structure.
+ * Basically calls Tesselation::GetDistanceToSurface() and checks the sign of the return value.
+ * @param point of which to check the position
+ * @param *LC LinkedCell structure
+ *
+ * @return true if the point is inside the Tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
+{
+  return (GetDistanceSquaredToSurface(Point, LC) < MYEPSILON);
+}
+
+/** Returns the distance to the surface given by the tesselation.
+ * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
+ * towards or away from the given \a &Point. Additionally, we check whether it's normal to the normal vector, i.e. on the
+ * closest triangle's plane. Then, we have to check whether \a Point is inside the triangle or not to determine whether it's
+ * an inside or outside point. This is done by calling BoundaryTriangleSet::GetIntersectionInsideTriangle().
+ * In the end we additionally find the point on the triangle who was smallest distance to \a Point:
+ *  -# Separate distance from point to center in vector in NormalDirection and on the triangle plane.
+ *  -# Check whether vector on triangle plane points inside the triangle or crosses triangle bounds.
+ *  -# If inside, take it to calculate closest distance
+ *  -# If not, take intersection with BoundaryLine as distance
+ *
+ * @note distance is squared despite it still contains a sign to determine in-/outside!
  *
  * @param point of which to check the position
  * @param *LC LinkedCell structure
  *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
-{
-  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
+ * @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__);
   Vector Center;
-
-  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(3) << "INFO: Central point of the triangle is " << Center << "." << endl;
-  Center.SubtractVector(&Point);
-  Log() << Verbose(3) << "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;
+  Vector helper;
+  Vector DistanceToCenter;
+  Vector Intersection;
+  double distance = 0.;
+
+  if (triangle == NULL) {// is boundary point or only point in point cloud?
+    Log() << Verbose(1) << "No triangle given!" << endl;
+    return -1.;
   } else {
-    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
-    return false;
-  }
-}
-
-/** Checks whether the provided TesselPoint is within the tesselation structure.
- *
- * @param *Point of which to check the position
- * @param *LC Linked Cell structure
- *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const
-{
-  return IsInnerPoint(*(Point->node), LC);
-}
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl;
+  }
+
+  triangle->GetCenter(&Center);
+  Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
+  DistanceToCenter.CopyVector(&Center);
+  DistanceToCenter.SubtractVector(&Point);
+  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl;
+
+  // check whether we are on boundary
+  if (fabs(DistanceToCenter.ScalarProduct(&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
+    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
+    if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
+      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
+      return 0.;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
+      return false;
+    }
+  } else {
+    // calculate smallest distance
+    distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection);
+    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
+
+    // then check direction to boundary
+    if (DistanceToCenter.ScalarProduct(&triangle->NormalVector) > MYEPSILON) {
+      Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
+      return -distance;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl;
+      return +distance;
+    }
+  }
+};
+
+/** Calculates distance 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::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
+{
+  BoundaryTriangleSet *triangle = FindClosestTriangleToVector(&Point, LC);
+  const double distance = GetDistanceSquaredToTriangle(Point, triangle);
+  return Min(distance, LC->RADIUS);
+};
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -2804 +3451 @@
  * @return set of the all points linked to the provided one
  */
-set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
-{
-  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
+TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+        TesselPointSet *connectedPoints = new TesselPointSet;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
   bool takePoint = false;
-
-  Log() << Verbose(3) << "Begin of GetAllConnectedPoints" << endl;
 
   // find the respective boundary point
@@ -2818 +3464 @@
     ReferencePoint = PointRunner->second;
   } else {
-    Log() << 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;
   }
@@ -2842 +3488 @@
 
    if (takePoint) {
-     Log() << Verbose(5) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
+     Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
      connectedPoints->insert(current);
    }
@@ -2849 +3495 @@
   }
 
-  if (connectedPoints->size() == 0) { // if have not found any points
+  if (connectedPoints->empty()) { // if have not found any points
     eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl;
     return NULL;
   }
 
-  Log() << Verbose(3) << "End of GetAllConnectedPoints" << endl;
   return connectedPoints;
 };
@@ -2866 +3511 @@
  *
  * @param *out output stream for debugging
+ * @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
  */
-list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(const TesselPoint* const Point, const Vector * const Reference) const
-{
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+{
+        Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(Point);
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  TesselPointList *connectedCircle = new TesselPointList;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
+  TriangleList *triangles = NULL;
+
+  if (SetOfNeighbours == NULL) {
+    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 {
+    eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl;
+    performCriticalExit();
+  }
+  PlaneNormal.Scale(1.0/triangles->size());
+  Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl;
+  PlaneNormal.Normalize();
+
+  // construct one orthogonal vector
+  if (Reference != NULL) {
+    AngleZero.CopyVector(Reference);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+  }
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
+    Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;
+    AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
+      performCriticalExit();
+    }
+  }
+  Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
+  if (AngleZero.NormSquared() > MYEPSILON)
+    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+  else
+    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+  Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
+
+  // go through all connected points and calculate angle
+  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);
+    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;
@@ -2881 +3609 @@
   Vector helper;
 
-  if (connectedPoints == NULL) {
-    Log() << Verbose(2) << "Could not find any connected points!" << endl;
+  if (SetOfNeighbours == NULL) {
+    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
     delete(connectedCircle);
     return NULL;
   }
-  Log() << Verbose(2) << "Begin of GetCircleOfConnectedPoints" << endl;
-
+
+  // check whether there's something to do
+  if (SetOfNeighbours->size() < 3) {
+    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+      connectedCircle->push_back(*TesselRunner);
+    return connectedCircle;
+  }
+
+  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
   // calculate central point
-  for (set<TesselPoint*>::const_iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
-    center.AddVector((*TesselRunner)->node);
+
+  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
+  TesselB++;
+  TesselC++;
+  TesselC++;
+  int counter = 0;
+  while (TesselC != SetOfNeighbours->end()) {
+    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
+    counter++;
+    TesselA++;
+    TesselB++;
+    TesselC++;
+    PlaneNormal.AddVector(&helper);
+  }
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
-  //  << "; scale factor " << 1.0/connectedPoints.size();
-  center.Scale(1.0/connectedPoints->size());
-  Log() << Verbose(4) << "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(4) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
+  //  << "; scale factor " << counter;
+  PlaneNormal.Scale(1.0/(double)counter);
+//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
+//
+//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
+//  PlaneNormal.CopyVector(Point->node);
+//  PlaneNormal.SubtractVector(&center);
+//  PlaneNormal.Normalize();
+  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
 
   // construct one orthogonal vector
@@ -2909 +3659 @@
   }
   if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
-    Log() << Verbose(4) << "Using alternatively " << *(*connectedPoints->begin())->node << " as angle 0 referencer." << endl;
-    AngleZero.CopyVector((*connectedPoints->begin())->node);
+    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);
@@ -2918 +3668 @@
     }
   }
-  Log() << Verbose(4) << "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);
-  Log() << Verbose(4) << "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 (set<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
+  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);
-    Log() << Verbose(3) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
-    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    if (angle > M_PI) // the correction is of no use here (and not desired)
+      angle = 2.*M_PI - angle;
+    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    if (!InserterTest.second) {
+      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
+      performCriticalExit();
+    }
   }
 
@@ -2938 +3695 @@
     connectedCircle->push_back(AngleRunner->second);
   }
-
-  delete(connectedPoints);
-
-  Log() << Verbose(2) << "End of GetCircleOfConnectedPoints" << endl;
 
   return connectedCircle;
@@ -2952 +3705 @@
  * @return list of the all points linked to the provided one
  */
-list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
+ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
+  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
@@ -2993 +3747 @@
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new list<TesselPoint*>;
+        connectedPath = new TesselPointList;
         triangle = NULL;
         CurrentLine = runner->second;
         StartLine = CurrentLine;
         CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
-        Log() << Verbose(3)<< "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
-          Log() << Verbose(3) << "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++) {
-            Log() << Verbose(3) << "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;
@@ -3013 +3767 @@
                 if (!TriangleRunner->second) {
                   TriangleRunner->second = true;
-                  Log() << Verbose(3) << "INFO: Connecting triangle is " << *triangle << "." << endl;
+                  Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl;
                   break;
                 } else {
-                  Log() << Verbose(3) << "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;
                 }
@@ -3031 +3785 @@
             if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
               CurrentLine = triangle->lines[i];
-              Log() << Verbose(3) << "INFO: Connecting line is " << *CurrentLine << "." << endl;
+              Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl;
               break;
             }
@@ -3045 +3799 @@
         } while (CurrentLine != StartLine);
         // last point is missing, as it's on start line
-        Log() << Verbose(3) << "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);
@@ -3051 +3805 @@
         ListOfPaths->push_back(connectedPath);
       } else {
-        Log() << Verbose(3) << "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;
       }
     }
@@ -3067 +3821 @@
  * @return list of the closed paths
  */
-list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
-  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point);
-  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
-  list<TesselPoint*> *newPath = NULL;
+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;
   int count = 0;
 
 
-  list<TesselPoint*>::iterator CircleRunner;
-  list<TesselPoint*>::iterator CircleStart;
-
-  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  TesselPointList::iterator CircleRunner;
+  TesselPointList::iterator CircleStart;
+
+  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
-    Log() << Verbose(2) << "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
@@ -3088 +3843 @@
 
     // go through list, look for reappearance of starting Point and create list
-    list<TesselPoint*>::iterator Marker = CircleStart;
+    TesselPointList::iterator Marker = CircleStart;
     for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
-        Log() << Verbose(3) << count+1 << ". closed path consists of: ";
-        newPath = new list<TesselPoint*>;
-        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        Log() << Verbose(1) << count+1 << ". closed path consists of: ";
+        newPath = new TesselPointList;
+        TesselPointList::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
@@ -3108 +3863 @@
     }
   }
-  Log() << Verbose(3) << "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
@@ -3128 +3883 @@
  * \return pointer to allocated list of triangles
  */
-set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
-{
-  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+        Info FunctionInfo(__func__);
+        TriangleSet *connectedTriangles = new TriangleSet;
 
   if (Point == NULL) {
@@ -3168 +3924 @@
     return 0.;
   } else
-    Log() << Verbose(2) << "Removing point " << *point << " from tesselated boundary ..." << endl;
+    Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl;
 
   // copy old location for the volume
@@ -3179 +3935 @@
   }
 
-  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
 
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  map<class BoundaryTriangleSet *, int> Candidates;
+  TriangleMap Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
     }
   }
@@ -3197 +3953 @@
   count=0;
   NormalVector.Zero();
-  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    Log() << Verbose(3) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
-    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
-    RemoveTesselationTriangle(Runner->first);
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl;
+    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->second);
     count++;
   }
   Log() << Verbose(1) << count << " triangles were removed." << endl;
 
-  list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
-  list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
-  map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
-  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+  TriangleMap::iterator NumberRunner = Candidates.begin();
+  TesselPointList::iterator StartNode, MiddleNode, EndNode;
   double angle;
   double smallestangle;
@@ -3222 +3978 @@
 
       // re-create all triangles by going through connected points list
-      list<class BoundaryLineSet *> NewLines;
+      LineList NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
@@ -3228 +3984 @@
         smallestangle = 0.;
         for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
-          Log() << Verbose(3) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+          Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
           // construct vectors to next and previous neighbour
           StartNode = MiddleNode;
@@ -3256 +4012 @@
         MiddleNode = EndNode;
         if (MiddleNode == connectedPath->end()) {
-          Log() << Verbose(1) << "CRITICAL: Could not find a smallest angle!" << endl;
-          exit(255);
+          eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl;
+          performCriticalExit();
         }
         StartNode = MiddleNode;
@@ -3266 +4022 @@
         if (EndNode == connectedPath->end())
           EndNode = connectedPath->begin();
-        Log() << Verbose(4) << "INFO: StartNode is " << **StartNode << "." << endl;
-        Log() << Verbose(4) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
-        Log() << Verbose(4) << "INFO: EndNode is " << **EndNode << "." << endl;
-        Log() << Verbose(3) << "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;
@@ -3275 +4031 @@
         triangle = GetPresentTriangle(TriangleCandidates);
         if (triangle != NULL) {
-          eLog() << Verbose(2) << "New triangle already present, skipping!" << endl;
+          eLog() << Verbose(0) << "New triangle already present, skipping!" << endl;
           StartNode++;
           MiddleNode++;
@@ -3287 +4043 @@
           continue;
         }
-        Log() << Verbose(5) << "Adding new triangle points."<< endl;
+        Log() << Verbose(3) << "Adding new triangle points."<< endl;
         AddTesselationPoint(*StartNode, 0);
         AddTesselationPoint(*MiddleNode, 1);
         AddTesselationPoint(*EndNode, 2);
-        Log() << Verbose(5) << "Adding new triangle lines."<< endl;
+        Log() << Verbose(3) << "Adding new triangle lines."<< endl;
         AddTesselationLine(TPS[0], TPS[1], 0);
         AddTesselationLine(TPS[0], TPS[2], 1);
@@ -3306 +4062 @@
         // prepare nodes for next triangle
         StartNode = EndNode;
-        Log() << Verbose(4) << "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
@@ -3313 +4069 @@
           break;
         } else if (connectedPath->size() < 2) { // something's gone wrong!
-          Log() << Verbose(1) << "CRITICAL: There are only two endpoints left!" << endl;
-          exit(255);
+          eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl;
+          performCriticalExit();
         } else {
           MiddleNode = StartNode;
@@ -3328 +4084 @@
       // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
       if (NewLines.size() > 1) {
-        list<class BoundaryLineSet *>::iterator Candidate;
+        LineList::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
-          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -3342 +4098 @@
           if (maxgain != 0) {
             volume += maxgain;
-            Log() << Verbose(3) << "Flipping baseline with highest volume" << **Candidate << "." << endl;
+            Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl;
             OtherBase = FlipBaseline(*Candidate);
             NewLines.erase(Candidate);
@@ -3353 +4109 @@
       delete(connectedPath);
     }
-    Log() << Verbose(1) << count << " triangles were created." << endl;
+    Log() << Verbose(0) << count << " triangles were created." << endl;
   } else {
     while (!ListOfClosedPaths->empty()) {
@@ -3361 +4117 @@
       delete(connectedPath);
     }
-    Log() << Verbose(1) << "No need to create any triangles." << endl;
+    Log() << Verbose(0) << "No need to create any triangles." << endl;
   }
   delete(ListOfClosedPaths);
 
-  Log() << Verbose(1) << "Removed volume is " << volume << "." << endl;
+  Log() << Verbose(0) << "Removed volume is " << volume << "." << endl;
 
   return volume;
@@ -3375 +4131 @@
  * Finds triangles belonging to the three provided points.
  *
- * @param *Points[3] list, is expected to contain three points
+ * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
  *
  * @return triangles which belong to the provided points, will be empty if there are none,
  *         will usually be one, in case of degeneration, there will be two
  */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
-{
-  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
+TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+        TriangleList *result = new TriangleList;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
   class BoundaryPointSet *TrianglePoints[3];
+  size_t NoOfWildcards = 0;
 
   for (int i = 0; i < 3; i++) {
-    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      TrianglePoints[i] = FindPoint->second;
+    if (Points[i] == NULL) {
+      NoOfWildcards++;
+      TrianglePoints[i] = NULL;
     } else {
-      TrianglePoints[i] = NULL;
-    }
-  }
-
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (TrianglePoints[i] != NULL) {
-      for (int j = i+1; j < 3; j++) {
-        if (TrianglePoints[j] != NULL) {
-          for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
-              (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
-              FindLine++) {
-            for (FindTriangle = FindLine->second->triangles.begin();
-                FindTriangle != FindLine->second->triangles.end();
-                FindTriangle++) {
-              if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
-                result->push_back(FindTriangle->second);
+      PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
+      if (FindPoint != PointsOnBoundary.end()) {
+        TrianglePoints[i] = FindPoint->second;
+      } else {
+        TrianglePoints[i] = NULL;
+      }
+    }
+  }
+
+  switch (NoOfWildcards) {
+    case 0: // checks lines between the points in the Points for their adjacent triangles
+      for (int i = 0; i < 3; i++) {
+        if (TrianglePoints[i] != NULL) {
+          for (int j = i+1; j < 3; j++) {
+            if (TrianglePoints[j] != NULL) {
+              for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
+                  (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
+                  FindLine++) {
+                for (FindTriangle = FindLine->second->triangles.begin();
+                    FindTriangle != FindLine->second->triangles.end();
+                    FindTriangle++) {
+                  if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+                    result->push_back(FindTriangle->second);
+                  }
+                }
               }
+              // Is it sufficient to consider one of the triangle lines for this.
+              return result;
             }
           }
-          // Is it sufficient to consider one of the triangle lines for this.
-          return result;
         }
       }
-    }
+      break;
+    case 1: // copy all triangles of the respective line
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] == NULL)
+          break;
+      for (FindLine = TrianglePoints[(i+1)%3]->lines.find(TrianglePoints[(i+2)%3]->node->nr); // is a multimap!
+          (FindLine != TrianglePoints[(i+1)%3]->lines.end()) && (FindLine->first == TrianglePoints[(i+2)%3]->node->nr);
+          FindLine++) {
+        for (FindTriangle = FindLine->second->triangles.begin();
+            FindTriangle != FindLine->second->triangles.end();
+            FindTriangle++) {
+          if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+            result->push_back(FindTriangle->second);
+          }
+        }
+      }
+      break;
+    }
+    case 2: // copy all triangles of the respective point
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] != NULL)
+          break;
+      for (LineMap::const_iterator line = TrianglePoints[i]->lines.begin(); line != TrianglePoints[i]->lines.end(); line++)
+        for (TriangleMap::const_iterator triangle = line->second->triangles.begin(); triangle != line->second->triangles.end(); triangle++)
+          result->push_back(triangle->second);
+      result->sort();
+      result->unique();
+      break;
+    }
+    case 3: // copy all triangles
+    {
+      for (TriangleMap::const_iterator triangle = TrianglesOnBoundary.begin(); triangle != TrianglesOnBoundary.end(); triangle++)
+        result->push_back(triangle->second);
+      break;
+    }
+    default:
+      eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl;
+      performCriticalExit();
+      break;
   }
 
   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>
 
 /**
@@ -3428 +4267 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedLines()
-{
-  map<int, class BoundaryLineSet *> AllLines;
-  map<int, int> * DegeneratedLines = new map<int, int>;
+IndexToIndex * Tesselation::FindAllDegeneratedLines()
+{
+        Info FunctionInfo(__func__);
+        UniqueLines AllLines;
+  IndexToIndex * DegeneratedLines = new IndexToIndex;
 
   // sanity check
   if (LinesOnBoundary.empty()) {
-    Log() << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
+    eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.";
     return DegeneratedLines;
   }
 
   LineMap::iterator LineRunner1;
-  pair<LineMap::iterator, bool> tester;
+  pair< UniqueLines::iterator, bool> tester;
   for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
-    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) );
+    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) );
     }
   }
@@ -3451 +4291 @@
   AllLines.clear();
 
-  Log() << Verbose(1) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
-  map<int,int>::iterator it;
-  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++)
-      Log() << Verbose(2) << (*it).first << " => " << (*it).second << endl;
+  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())
+      Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl;
+    else
+      eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl;
+  }
 
   return DegeneratedLines;
@@ -3465 +4311 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedTriangles()
-{
-  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
-  map<int, int> * DegeneratedTriangles = new map<int, int>;
+IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
+{
+        Info FunctionInfo(__func__);
+  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
+  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
 
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
@@ -3474 +4321 @@
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
 
-  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
@@ -3494 +4341 @@
   delete(DegeneratedLines);
 
-  Log() << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
+  IndexToIndex::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
-      Log() << Verbose(2) << (*it).first << " => " << (*it).second << endl;
+      Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
 
   return DegeneratedTriangles;
@@ -3508 +4355 @@
 void Tesselation::RemoveDegeneratedTriangles()
 {
-  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+        Info FunctionInfo(__func__);
+  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
 
-  Log() << Verbose(1) << "Begin of RemoveDegeneratedTriangles" << endl;
-
-  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
@@ -3572 +4418 @@
       // erase the pair
       count += (int) DegeneratedTriangles->erase(triangle->Nr);
-      Log() << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;
+      Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;
       RemoveTesselationTriangle(triangle);
       count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
-      Log() << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;
+      Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;
       RemoveTesselationTriangle(partnerTriangle);
     } else {
-      Log() << Verbose(1) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle
+      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;
@@ -3584 +4430 @@
   }
   delete(DegeneratedTriangles);
-
-  Log() << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl;
-  Log() << Verbose(1) << "End of RemoveDegeneratedTriangles" << endl;
+  if (count > 0)
+    LastTriangle = NULL;
+
+  Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl;
 }
 
@@ -3599 +4446 @@
 void Tesselation::AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC)
 {
-  Log() << Verbose(2) << "Begin of AddBoundaryPointByDegeneratedTriangle" << endl;
-
+        Info FunctionInfo(__func__);
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
@@ -3616 +4462 @@
     return;
   }
-  Log() << Verbose(2) << "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
@@ -3649 +4495 @@
 
   // create new triangle to connect point (connects automatically with the missing spot of the chosen line)
-  Log() << Verbose(5) << "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);
-  Log() << Verbose(5) << "Adding new triangle lines."<< endl;
+  Log() << Verbose(2) << "Adding new triangle lines."<< endl;
   AddTesselationLine(TPS[0], TPS[1], 0);
   AddTesselationLine(TPS[0], TPS[2], 1);
@@ -3660 +4506 @@
   BTS->GetNormalVector(TempTriangle->NormalVector);
   BTS->NormalVector.Scale(-1.);
-  Log() << Verbose(3) << "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
-  Log() << Verbose(5) << "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);
-  Log() << Verbose(5) << "Adding new triangle lines."<< endl;
+  Log() << Verbose(2) << "Adding new triangle lines."<< endl;
   AddTesselationLine(TPS[0], TPS[1], 0);
   AddTesselationLine(TPS[0], TPS[2], 1);
@@ -3674 +4520 @@
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   BTS->GetNormalVector(TempTriangle->NormalVector);
-  Log() << Verbose(3) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;
+  Log() << Verbose(1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;
   AddTesselationTriangle();
 
@@ -3681 +4527 @@
     if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
       if (BestLine == BTS->lines[i]){
-        Log() << Verbose(1) << "CRITICAL: BestLine is same as found line, something's wrong here!" << endl;
-        exit(255);
+        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) );
@@ -3689 +4535 @@
     }
   }
-
-  // exit
-  Log() << Verbose(2) << "End of AddBoundaryPointByDegeneratedTriangle" << endl;
 };
 
@@ -3701 +4544 @@
 void Tesselation::Output(const char *filename, const PointCloud * const cloud)
 {
+        Info FunctionInfo(__func__);
   ofstream *tempstream = NULL;
   string NameofTempFile;
@@ -3713 +4557 @@
       NameofTempFile.erase(npos, 1);
       NameofTempFile.append(TecplotSuffix);
-      Log() << Verbose(1) << "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);
@@ -3727 +4571 @@
       NameofTempFile.erase(npos, 1);
       NameofTempFile.append(Raster3DSuffix);
-      Log() << Verbose(1) << "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);
@@ -3739 +4583 @@
     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++) {
+    Log() << Verbose(0) << "Current point is " << *Runner->second << "." << endl;
+    map < int, Vector *> TriangleVectors;
+    // gather all NormalVectors
+    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)
+            Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        } else {
+          Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl;
+        }
+      }
+    // check whether there are two that are parallel
+    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
+          Log() << Verbose(1) << "Checking " << *VectorWalker->second<< " against " << *VectorRunner->second << ": " << SCP << endl;
+          if (fabs(SCP + 1.) < ParallelEpsilon) {
+            InsertionTester = EndpointCandidateList.insert((Runner->second));
+            if (InsertionTester.second)
+              Log() << Verbose(0) << " Adding " << *Runner->second << " to endpoint candidate list." << endl;
+            // and break out of both loops
+            VectorWalker = TriangleVectors.end();
+            VectorRunner = TriangleVectors.end();
+            break;
+          }
+        }
+  }
+
+  /// 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
+      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);
+        Log() << Verbose(1) << "Checking " << *OtherWalker << endl;
+        set < BoundaryPointSet *>::iterator Finder = EndpointCandidateList.find(OtherWalker);
+        if (Finder != EndpointCandidateList.end()) {  // found a connected partner
+          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 {
+          Log() << Verbose(1) << " is not connected to " << *Walker << endl;
+        }
+      }
+    }
+
+    Log() << Verbose(0) << "Final polygon is " << *Current << endl;
+    ListofDegeneratedPolygons.insert(Current);
+    Current = NULL;
+  }
+
+  const int counter = ListofDegeneratedPolygons.size();
+
+  Log() << Verbose(0) << "The following " << counter << " degenerated polygons have been found: " << endl;
+  for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++)
+    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) {
+      Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl;
+      delete(T);
+      continue;
+    }
+
+    // check whether number is even
+    // If this case occurs, we have to think about it!
+    // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
+    // connections to either polygon ...
+    if (T->size() % 2 != 0) {
+      eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl;
+      performCriticalExit();
+    }
+
+    TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator
+    /// 4a. Get NormalVector for one side (this is "front")
+    NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
+    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++;
+      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
+        Log() << Verbose(1) << " Removing ... " << endl;
+        TriangleNrs.push(triangle->Nr);
+        T->erase(TriangleWalker);
+        RemoveTesselationTriangle(triangle);
+      } else
+        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
+      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(TPS[0], TPS[1], 0);
+      AddTesselationLine(TPS[0], TPS[2], 1);
+      AddTesselationLine(TPS[1], TPS[2], 2);
+      if (TriangleNrs.empty())
+        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()) {
+      eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl;
+    }
+    delete(T);  // remove the triangleset
+  }
+
+  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  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++)
+      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;
+};