Changeset 6613ec

Timestamp:

Apr 21, 2010, 12:07:01 PM (15 years ago)

Author:

Frederik Heber <heber@…>

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, 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:

a67d19

Parents:

bc84ffc

git-author:

Frederik Heber <heber@…> (04/21/10 11:51:53)

git-committer:

Frederik Heber <heber@…> (04/21/10 12:07:01)

Message:

Rewrite of tesselation: treatment of degenerated triangles and polygons

several tricks have been necessary to make it work for heptane:

• AddDegeneratedTriangle() is not used, instead AddCandidateTriangle() gets told which OptCenter to use, i.e. which triangle side to add
• afterwards, in the case of a polygon (multiple candidates for a baseline) we have to set the candidate of internal edges in the polygon to the next point, such that AddTesselationLine() picks the right one
• The we fill other open lines, needed in case of a degenerated triangle
• If CheckDegeneracy() - which removes all candidates from the ListOfPoints, not just the first - results in true
• Then we call AddCandidateTriangle() for the other side of it
• and again in the case of a polygon, set the candidate for the internal edge accordingly.

• new
  • Tesselation::FindCandidatesforOpenLines() - checks for OpenLines with no candidates and calls FindNextSuitableTriangle() for each
  • Tesselation::FindDegeneratedCandidatesforOpenLines() - sets the candidate of an open line between two points to a given point (for above polygons)
  • enum centers is either Opt or OtherOpt for normal or back side of the triangle, used in AddCandidateTriangle()

• associated changes:
  • Tesselation::CheckDegeneracy() - needs current CandidateForTesselation for Baseline's endpoints and candidate list (not just the current trianglre in TPS), removes these instead of TPS
  • lots of identation changes made by eclipse, bust mostly sensible, due to coding on other laptop before
  • Tesselation::AddCandidateTriangle() - picks OptCenter according to given enum centers
  • Tesselation::AddTesselationLine() - in polygon case, a base line may have multiple candidates, though with same OptCenter, we check all of these against desired ThirdNode
  • FindNonConvexBorder() - removed call of CorrectAllDegeneratedPolygons(), as they should be created correctly now

Signed-off-by: Frederik Heber <heber@…>

Location:

src

Files:

• boundary.cpp (modified) (3 diffs)
• tesselation.cpp (modified) (289 diffs)
• tesselation.hpp (modified) (5 diffs)

src/boundary.cpp

@@ -960 +960 @@
   bool freeLC = false;
   bool status = false;
-  CandidateForTesselation *baseline;
-  LineMap::iterator testline;
+  CandidateForTesselation *baseline = NULL;
   bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
   bool TesselationFailFlag = false;
-  BoundaryTriangleSet *T = NULL;
 
   if (TesselStruct == NULL) {
@@ -997 +995 @@
 
     // 2b. find best candidate for each OpenLine
-    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
-      baseline = Runner->second;
-      if (baseline->pointlist.empty()) {
-        T = (((baseline->BaseLine->triangles.begin()))->second);
-        Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
-        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
-      }
-    }
+    TesselationFailFlag = TesselStruct->FindCandidatesforOpenLines(RADIUS, LCList);
 
     // 2c. print OpenLines with candidates again
@@ -1062 +1053 @@
   StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
 
-  // correct degenerated polygons
-  TesselStruct->CorrectAllDegeneratedPolygons();
-
-  // check envelope for consistency
-  status = CheckListOfBaselines(TesselStruct);
+//  // correct degenerated polygons
+//  TesselStruct->CorrectAllDegeneratedPolygons();
+//
+//  // check envelope for consistency
+//  status = CheckListOfBaselines(TesselStruct);
 
   // write final envelope

src/tesselation.cpp

@@ -25 +25 @@
  */
 BoundaryPointSet::BoundaryPointSet() :
-    LinesCount(0),
-    value(0.),
-    Nr(-1)
-{
-        Info FunctionInfo(__func__);
-        Log() << Verbose(1) << "Adding noname." << endl;
-};
+  LinesCount(0), value(0.), Nr(-1)
+{
+  Info FunctionInfo(__func__);
+  Log() << Verbose(1) << "Adding noname." << endl;
+}
+;
 
 /** Constructor of BoundaryPointSet with Tesselpoint.
@@ -37 +36 @@
  */
 BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
-  LinesCount(0),
-  node(Walker),
-  value(0.),
-  Nr(Walker->nr)
-{
-        Info FunctionInfo(__func__);
+  LinesCount(0), node(Walker), value(0.), Nr(Walker->nr)
+{
+  Info FunctionInfo(__func__);
   Log() << Verbose(1) << "Adding Node " << *Walker << endl;
-};
+}
+;
 
 /** Destructor of BoundaryPointSet.
@@ -52 +49 @@
 BoundaryPointSet::~BoundaryPointSet()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl;
   if (!lines.empty())
-    DoeLog(2) && (eLog()<< Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);
+    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);
   node = NULL;
-};
+}
+;
 
 /** Add a line to the LineMap of this point.
@@ -64 +62 @@
 void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
 {
-        Info FunctionInfo(__func__);
-  Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "."
-      << endl;
-  if (line->endpoints[0] == this)
-    {
-      lines.insert(LinePair(line->endpoints[1]->Nr, line));
-    }
-  else
-    {
-      lines.insert(LinePair(line->endpoints[0]->Nr, line));
-    }
+  Info FunctionInfo(__func__);
+  Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." << endl;
+  if (line->endpoints[0] == this) {
+    lines.insert(LinePair(line->endpoints[1]->Nr, line));
+  } else {
+    lines.insert(LinePair(line->endpoints[0]->Nr, line));
+  }
   LinesCount++;
-};
+}
+;
 
 /** output operator for BoundaryPointSet.
@@ -94 +89 @@
  */
 BoundaryLineSet::BoundaryLineSet() :
-    Nr(-1)
-{
-        Info FunctionInfo(__func__);
+  Nr(-1)
+{
+  Info FunctionInfo(__func__);
   for (int i = 0; i < 2; i++)
     endpoints[i] = NULL;
-};
+}
+;
 
 /** Constructor of BoundaryLineSet with two endpoints.
@@ -108 +104 @@
 BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // set number
   Nr = number;
@@ -120 +116 @@
   // clear triangles list
   Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
-};
+}
+;
 
 /** Constructor of BoundaryLineSet with two endpoints.
@@ -142 +139 @@
   // clear triangles list
   Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
-};
+}
+;
 
 /** Destructor for BoundaryLineSet.
@@ -150 +148 @@
 BoundaryLineSet::~BoundaryLineSet()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   int Numbers[2];
 
@@ -181 +179 @@
         //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         if (endpoints[i] != NULL) {
-          delete(endpoints[i]);
+          delete (endpoints[i]);
           endpoints[i] = NULL;
         }
@@ -188 +186 @@
   }
   if (!triangles.empty())
-    DoeLog(2) && (eLog()<< Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl);
-};
+    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl);
+}
+;
 
 /** Add triangle to TriangleMap of this boundary line.
@@ -196 +195 @@
 void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
   triangles.insert(TrianglePair(triangle->Nr, triangle));
-};
+}
+;
 
 /** Checks whether we have a common endpoint with given \a *line.
@@ -207 +207 @@
 bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
     return true;
   else
     return false;
-};
+}
+;
 
 /** Checks whether the adjacent triangles of a baseline are convex or not.
@@ -222 +223 @@
 bool BoundaryLineSet::CheckConvexityCriterion() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
   // get the two triangles
   if (triangles.size() != 2) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);
+    DoeLog(0) && (eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);
     return true;
   }
@@ -234 +235 @@
   BaseLineCenter.CopyVector(endpoints[0]->node->node);
   BaseLineCenter.AddVector(endpoints[1]->node->node);
-  BaseLineCenter.Scale(1./2.);
+  BaseLineCenter.Scale(1. / 2.);
   BaseLine.CopyVector(endpoints[0]->node->node);
   BaseLine.SubtractVector(endpoints[1]->node->node);
@@ -242 +243 @@
   NormalCheck.Zero();
   double sign = -1.;
-  int i=0;
+  int i = 0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for (TriangleMap::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);
@@ -250 +251 @@
     sign = -sign;
     if (runner->second->NormalVector.NormSquared() > MYEPSILON)
-      BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
+      BaseLineNormal.CopyVector(&runner->second->NormalVector); // yes, copy second on top of first
     else {
-      DoeLog(0) && (eLog()<< Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
+      DoeLog(0) && (eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
     }
     node = runner->second->GetThirdEndpoint(this);
@@ -259 +260 @@
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
-      helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
+      helper[i].MakeNormalVector(&BaseLine); // we want to compare the triangle's heights' angles!
       //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
     } else {
-      DoeLog(1) && (eLog()<< Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);
       return true;
     }
@@ -289 +290 @@
 bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
-        Info FunctionInfo(__func__);
-  for(int i=0;i<2;i++)
+  Info FunctionInfo(__func__);
+  for (int i = 0; i < 2; i++)
     if (point == endpoints[i])
       return true;
   return false;
-};
+}
+;
 
 /** Returns other endpoint of the line.
@@ -302 +304 @@
 class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   if (endpoints[0] == point)
     return endpoints[1];
@@ -309 +311 @@
   else
     return NULL;
-};
+}
+;
 
 /** output operator for BoundaryLineSet.
@@ -315 +318 @@
  * \param &a boundary line
  */
-ostream & operator <<(ostream &ost, const  BoundaryLineSet &a)
+ostream & operator <<(ostream &ost, const BoundaryLineSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]";
   return ost;
-};
+}
+;
 
 // ======================================== Triangles on Boundary =================================
@@ -328 +332 @@
   Nr(-1)
 {
-        Info FunctionInfo(__func__);
-  for (int i = 0; i < 3; i++)
-    {
-      endpoints[i] = NULL;
-      lines[i] = NULL;
-    }
-};
+  Info FunctionInfo(__func__);
+  for (int i = 0; i < 3; i++) {
+    endpoints[i] = NULL;
+    lines[i] = NULL;
+  }
+}
+;
 
 /** Constructor for BoundaryTriangleSet with three lines.
@@ -343 +347 @@
   Nr(number)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // set number
   // set lines
@@ -355 +359 @@
     // for all three lines
     for (int j = 0; j < 2; j++) { // for both endpoints
-      OrderMap.insert(pair<int, class BoundaryPointSet *> (
-          line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
+      OrderMap.insert(pair<int, class BoundaryPointSet *> (line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
       // and we don't care whether insertion fails
     }
@@ -368 +371 @@
   }
   if (Counter < 3) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl);
+    DoeLog(0) && (eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl);
     performCriticalExit();
   }
-};
+}
+;
 
 /** Destructor of BoundaryTriangleSet.
@@ -379 +383 @@
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   for (int i = 0; i < 3; i++) {
     if (lines[i] != NULL) {
@@ -386 +390 @@
       }
       if (lines[i]->triangles.empty()) {
-          //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-          delete (lines[i]);
-          lines[i] = NULL;
+        //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+        delete (lines[i]);
+        lines[i] = NULL;
       }
     }
   }
   //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl;
-};
+}
+;
 
 /** Calculates the normal vector for this triangle.
@@ -401 +406 @@
 void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // get normal vector
   NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);
@@ -409 +414 @@
     NormalVector.Scale(-1.);
   Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl;
-};
+}
+;
 
 /** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
@@ -430 +436 @@
 
   if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
     return false;
   }
@@ -441 +447 @@
     Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl;
     return true;
-  }   else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
+  } 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) {
+  } 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;
+  int i = 0;
   do {
-    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
-      helper.CopyVector(endpoints[(i+1)%3]->node->node);
-      helper.SubtractVector(endpoints[i%3]->node->node);
-      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-      const double s = CrossPoint.ScalarProduct(&helper)/helper.NormSquared();
+    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i % 3]->node->node, endpoints[(i + 1) % 3]->node->node, endpoints[(i + 2) % 3]->node->node, Intersection, &NormalVector)) {
+      helper.CopyVector(endpoints[(i + 1) % 3]->node->node);
+      helper.SubtractVector(endpoints[i % 3]->node->node);
+      CrossPoint.SubtractVector(endpoints[i % 3]->node->node); // cross point was returned as absolute vector
+      const double s = CrossPoint.ScalarProduct(&helper) / helper.NormSquared();
       Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl;
-      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
+      if ((s < -MYEPSILON) || ((s - 1.) > MYEPSILON)) {
         Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl;
-        i=4;
+        i = 4;
         break;
       }
       i++;
-    } else 
+    } else
       break;
-  } while (i<3);
-  if (i==3) {
+  } while (i < 3);
+  if (i == 3) {
     Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
     return true;
@@ -473 +479 @@
     return false;
   }
-};
+}
+;
 
 /** Finds the point on the triangle to the point \a *x.
@@ -502 +509 @@
   Vector InPlane;
   InPlane.CopyVector(x);
-  InPlane.SubtractVector(ClosestPoint);  // points from plane intersection to straight-down point
+  InPlane.SubtractVector(ClosestPoint); // points from plane intersection to straight-down point
   InPlane.ProjectOntoPlane(&NormalVector);
   InPlane.AddVector(ClosestPoint);
@@ -516 +523 @@
   Vector CrossPoint[3];
   Vector helper;
-  for (int i=0;i<3;i++) {
+  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);
+    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);
+    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();
+    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;
+    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)) {
@@ -539 +546 @@
   }
   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
+  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;
   }
@@ -548 +556 @@
     ClosestPoint->CopyVector(&InPlane);
     ShortestDistance = InPlane.DistanceSquared(x);
-  } else {  // also check endnodes
-    for (int i=0;i<3;i++) {
+  } 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)) {
@@ -559 +567 @@
   Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl;
   return ShortestDistance;
-};
+}
+;
 
 /** Checks whether lines is any of the three boundary lines this triangle contains.
@@ -567 +576 @@
 bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
 {
-        Info FunctionInfo(__func__);
-  for(int i=0;i<3;i++)
+  Info FunctionInfo(__func__);
+  for (int i = 0; i < 3; i++)
     if (line == lines[i])
       return true;
   return false;
-};
+}
+;
 
 /** Checks whether point is any of the three endpoints this triangle contains.
@@ -580 +590 @@
 bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
-        Info FunctionInfo(__func__);
-  for(int i=0;i<3;i++)
+  Info FunctionInfo(__func__);
+  for (int i = 0; i < 3; i++)
     if (point == endpoints[i])
       return true;
   return false;
-};
+}
+;
 
 /** Checks whether point is any of the three endpoints this triangle contains.
@@ -593 +604 @@
 bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
 {
-        Info FunctionInfo(__func__);
-  for(int i=0;i<3;i++)
+  Info FunctionInfo(__func__);
+  for (int i = 0; i < 3; i++)
     if (point == endpoints[i]->node)
       return true;
   return false;
-};
+}
+;
 
 /** Checks whether three given \a *Points coincide with triangle's endpoints.
@@ -606 +618 @@
 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])
-            || (endpoints[0] == Points[2])
-          ) && (
-            (endpoints[1] == Points[0])
-            || (endpoints[1] == Points[1])
-            || (endpoints[1] == Points[2])
-          ) && (
-            (endpoints[2] == Points[0])
-            || (endpoints[2] == Points[1])
-            || (endpoints[2] == Points[2])
-
-          ));
-};
+  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]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2])
+
+  ));
+}
+;
 
 /** Checks whether three given \a *Points coincide with triangle's endpoints.
@@ -629 +632 @@
 bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
 {
-        Info FunctionInfo(__func__);
-  return (((endpoints[0] == T->endpoints[0])
-            || (endpoints[0] == T->endpoints[1])
-            || (endpoints[0] == T->endpoints[2])
-          ) && (
-            (endpoints[1] == T->endpoints[0])
-            || (endpoints[1] == T->endpoints[1])
-            || (endpoints[1] == T->endpoints[2])
-          ) && (
-            (endpoints[2] == T->endpoints[0])
-            || (endpoints[2] == T->endpoints[1])
-            || (endpoints[2] == T->endpoints[2])
-
-          ));
-};
+  Info FunctionInfo(__func__);
+  return (((endpoints[0] == T->endpoints[0]) || (endpoints[0] == T->endpoints[1]) || (endpoints[0] == T->endpoints[2])) && ((endpoints[1] == T->endpoints[0]) || (endpoints[1] == T->endpoints[1]) || (endpoints[1] == T->endpoints[2])) && ((endpoints[2] == T->endpoints[0]) || (endpoints[2] == T->endpoints[1]) || (endpoints[2] == T->endpoints[2])
+
+  ));
+}
+;
 
 /** Returns the endpoint which is not contained in the given \a *line.
@@ -651 +645 @@
 class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // sanity check
   if (!ContainsBoundaryLine(line))
     return NULL;
-  for(int i=0;i<3;i++)
+  for (int i = 0; i < 3; i++)
     if (!line->ContainsBoundaryPoint(endpoints[i]))
       return endpoints[i];
   // actually, that' impossible :)
   return NULL;
-};
+}
+;
 
 /** Calculates the center point of the triangle.
@@ -668 +663 @@
 void BoundaryTriangleSet::GetCenter(Vector * const center) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   center->Zero();
-  for(int i=0;i<3;i++)
+  for (int i = 0; i < 3; i++)
     center->AddVector(endpoints[i]->node->node);
-  center->Scale(1./3.);
+  center->Scale(1. / 3.);
   Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
@@ -683 +678 @@
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
-//  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
-//      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
+  //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
+  //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
   return ost;
-};
+}
+;
 
 // ======================================== Polygons on Boundary =================================
@@ -696 +692 @@
 {
   Info FunctionInfo(__func__);
-};
+}
+;
 
 /** Destructor of BoundaryPolygonSet.
@@ -707 +704 @@
   endpoints.clear();
   Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl;
-};
+}
+;
 
 /** Calculates the normal vector for this triangle.
@@ -721 +719 @@
   Vector *TotalNormal = new Vector;
   PointSet::const_iterator Runner[3];
-  for (int i=0;i<3; i++) {
+  for (int i = 0; i < 3; i++) {
     Runner[i] = endpoints.begin();
-    for (int j = 0; j<i; j++) { // go as much further
+    for (int j = 0; j < i; j++) { // go as much further
       Runner[i]++;
       if (Runner[i] == endpoints.end()) {
-        DoeLog(0) && (eLog()<< Verbose(0) << "There are less than three endpoints in the polygon!" << endl);
+        DoeLog(0) && (eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl);
         performCriticalExit();
       }
@@ -732 +730 @@
   }
   TotalNormal->Zero();
-  int counter=0;
-  for (; Runner[2] != endpoints.end(); ) {
+  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
+    for (int i = 0; i < 3; i++) // increase each of them
       Runner[i]++;
     TotalNormal->AddVector(&TemporaryNormal);
   }
-  TotalNormal->Scale(1./(double)counter);
+  TotalNormal->Scale(1. / (double) counter);
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
@@ -747 +745 @@
 
   return TotalNormal;
-};
+}
+;
 
 /** Calculates the center point of the triangle.
@@ -758 +757 @@
   center->Zero();
   int counter = 0;
-  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
     center->AddVector((*Runner)->node->node);
     counter++;
   }
-  center->Scale(1./(double)counter);
+  center->Scale(1. / (double) counter);
   Log() << Verbose(1) << "Center is at " << *center << "." << endl;
 }
@@ -774 +773 @@
   Info FunctionInfo(__func__);
   return ContainsPresentTupel(triangle->endpoints, 3);
-};
+}
+;
 
 /** Checks whether the polygons contains both endpoints of the line.
@@ -784 +784 @@
   Info FunctionInfo(__func__);
   return ContainsPresentTupel(line->endpoints, 2);
-};
+}
+;
 
 /** Checks whether point is any of the three endpoints this triangle contains.
@@ -793 +794 @@
 {
   Info FunctionInfo(__func__);
-  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
     Log() << Verbose(0) << "Checking against " << **Runner << endl;
     if (point == (*Runner)) {
@@ -802 +803 @@
   Log() << Verbose(0) << " Not contained." << endl;
   return false;
-};
+}
+;
 
 /** Checks whether point is any of the three endpoints this triangle contains.
@@ -811 +813 @@
 {
   Info FunctionInfo(__func__);
-  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
+  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
     if (point == (*Runner)->node) {
       Log() << Verbose(0) << " Contained." << endl;
@@ -818 +820 @@
   Log() << Verbose(0) << " Not contained." << endl;
   return false;
-};
+}
+;
 
 /** Checks whether given array of \a *Points coincide with polygons's endpoints.
@@ -830 +833 @@
   int counter = 0;
   Log() << Verbose(1) << "Polygon is " << *this << endl;
-  for(int i=0;i<dim;i++) {
+  for (int i = 0; i < dim; i++) {
     Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl;
     if (ContainsBoundaryPoint(Points[i])) {
@@ -841 +844 @@
   else
     return false;
-};
+}
+;
 
 /** Checks whether given PointList coincide with polygons's endpoints.
@@ -852 +856 @@
   size_t counter = 0;
   Log() << Verbose(1) << "Polygon is " << *this << endl;
-  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
+  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
     Log() << Verbose(1) << " Testing endpoint " << **Runner << endl;
     if (ContainsBoundaryPoint(*Runner))
@@ -862 +866 @@
   else
     return false;
-};
+}
+;
 
 /** Checks whether given set of \a *Points coincide with polygons's endpoints.
@@ -870 +875 @@
 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const
 {
-  return ContainsPresentTupel((const PointSet)P->endpoints);
-};
+  return ContainsPresentTupel((const PointSet) P->endpoints);
+}
+;
 
 /** Gathers all the endpoints' triangles in a unique set.
@@ -879 +885 @@
 {
   Info FunctionInfo(__func__);
-  pair <TriangleSet::iterator, bool> Tester;
+  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++) {
+  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)) {
@@ -895 +901 @@
   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.
@@ -904 +911 @@
 {
   Info FunctionInfo(__func__);
-  pair <PointSet::iterator, bool> Tester;
+  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++) {
+  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)
@@ -917 +924 @@
 
   return true;
-};
+}
+;
 
 /** output operator for BoundaryPolygonSet.
@@ -926 +934 @@
 {
   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<< "]";
+  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 ===========================================
@@ -945 +954 @@
   node = NULL;
   nr = -1;
-  Name =  NULL;
-};
+  Name = NULL;
+}
+;
 
 /** Destructor for class TesselPoint.
@@ -953 +963 @@
 {
   //Info FunctionInfo(__func__);
-};
+}
+;
 
 /** Prints LCNode to screen.
  */
-ostream & operator << (ostream &ost, const TesselPoint &a)
+ostream & operator <<(ostream &ost, const TesselPoint &a)
 {
   ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
   return ost;
-};
+}
+;
 
 /** Prints LCNode to screen.
  */
-ostream & TesselPoint::operator << (ostream &ost)
-{
-        Info FunctionInfo(__func__);
+ostream & TesselPoint::operator <<(ostream &ost)
+{
+  Info FunctionInfo(__func__);
   ost << "[" << (nr) << "|" << this << "]";
   return ost;
-};
-
+}
+;
 
 // =========================================================== class POINTCLOUD ============================================
@@ -979 +991 @@
 PointCloud::PointCloud()
 {
-        //Info FunctionInfo(__func__);
-};
+  //Info FunctionInfo(__func__);
+}
+;
 
 /** Destructor for class PointCloud.
@@ -986 +999 @@
 PointCloud::~PointCloud()
 {
-        //Info FunctionInfo(__func__);
-};
+  //Info FunctionInfo(__func__);
+}
+;
 
 // ============================ CandidateForTesselation =============================
@@ -993 +1007 @@
 /** Constructor of class CandidateForTesselation.
  */
-CandidateForTesselation::CandidateForTesselation (BoundaryLineSet* line) :
-  BaseLine(line),
-  ThirdPoint(NULL),
-  T(NULL),
-  ShortestAngle(2.*M_PI),
-  OtherShortestAngle(2.*M_PI)
-{
-        Info FunctionInfo(__func__);
-};
-
+CandidateForTesselation::CandidateForTesselation(BoundaryLineSet* line) :
+  BaseLine(line), ThirdPoint(NULL), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
+{
+  Info FunctionInfo(__func__);
+}
+;
 
 /** Constructor of class CandidateForTesselation.
  */
-CandidateForTesselation::CandidateForTesselation (TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
-    BaseLine(line),
-    ThirdPoint(point),
-    T(NULL),
-    ShortestAngle(2.*M_PI),
-    OtherShortestAngle(2.*M_PI)
-{
-        Info FunctionInfo(__func__);
+CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
+  BaseLine(line), ThirdPoint(point), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
+{
+  Info FunctionInfo(__func__);
   OptCenter.CopyVector(&OptCandidateCenter);
   OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
-};
+}
+;
 
 /** Destructor for class CandidateForTesselation.
  */
-CandidateForTesselation::~CandidateForTesselation() {
-};
+CandidateForTesselation::~CandidateForTesselation()
+{
+}
+;
 
 /** Checks validity of a given sphere of a candidate line.
@@ -1033 +1042 @@
   Info FunctionInfo(__func__);
 
-  const double radiusSquared = RADIUS*RADIUS;
+  const double radiusSquared = RADIUS * RADIUS;
   list<const Vector *> VectorList;
   VectorList.push_back(&OptCenter);
@@ -1039 +1048 @@
 
   if (!pointlist.empty())
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list " << *(*pointlist.begin()) << " and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
   else
     DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
   // check baseline for OptCenter and OtherOptCenter being on sphere's surface
   for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-    for (int i=0;i<2;i++) {
+    for (int i = 0; i < 2; i++) {
       const double distance = fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->node) - radiusSquared);
       if (distance > HULLEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << setprecision(13) << distance << "." << endl);
+        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
         return false;
       }
@@ -1054 +1063 @@
 
   // check Candidates for OptCenter and OtherOptCenter being on sphere's surface
-  for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.begin(); ++Runner) {
+  for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
     const TesselPoint *Walker = *Runner;
     for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
       const double distance = fabs((*VRunner)->DistanceSquared(Walker->node) - radiusSquared);
       if (distance > HULLEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << Walker << " is out of sphere at " << *(*VRunner) << " by " << setprecision(13) << distance << "." << endl);
+        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
         return false;
+      } else {
+        Log() << Verbose(1) << "Candidate " << *Walker << " is inside by " << distance << "." << endl;
       }
     }
@@ -1070 +1081 @@
     // get all points inside the sphere
     TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));
+
+    Log() << Verbose(1) << "The following atoms are inside sphere at " << OtherOptCenter << ":" << endl;
+    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
+      Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&OtherOptCenter) << "." << endl;
+
     // remove baseline's endpoints and candidates
-    for (int i=0;i<2;i++)
+    for (int i = 0; i < 2; i++) {
+      Log() << Verbose(1) << "INFO: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "." << endl;
       ListofPoints->remove(BaseLine->endpoints[i]->node);
-    for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner)
+    }
+    for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
+      Log() << Verbose(1) << "INFO: removing candidate tesselpoint " << *(*Runner) << "." << endl;
       ListofPoints->remove(*Runner);
+    }
     if (!ListofPoints->empty()) {
-      cout << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl;
+      DoeLog(1) && (eLog() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
       flag = false;
       DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl);
@@ -1082 +1102 @@
         DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << endl);
     }
-    delete(ListofPoints);
+    delete (ListofPoints);
 
     // check with animate_sphere.tcl VMD script
     if (ThirdPoint != NULL) {
-      cout << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr+1 << " " << BaseLine->endpoints[1]->Nr+1 << " " << ThirdPoint->Nr+1 << " " << RADIUS << " ";
-      cout << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " ";
-      cout << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl;
+      Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl;
     } else {
-      cout << Verbose(1) << "Check by: ... missing third point ..." << endl;
-      cout << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr+1 << " " << BaseLine->endpoints[1]->Nr+1 << " ??? " << RADIUS << " ";
-      cout << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " ";
-      cout << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl;
+      Log() << Verbose(1) << "Check by: ... missing third point ..." << endl;
+      Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl;
     }
   }
   return flag;
-};
+}
+;
 
 /** output operator for CandidateForTesselation.
@@ -1103 +1120 @@
  * \param &a boundary line
  */
-ostream & operator <<(ostream &ost, const  CandidateForTesselation &a)
+ostream & operator <<(ostream &ost, const CandidateForTesselation &a)
 {
   ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->Name << "," << a.BaseLine->endpoints[1]->node->Name << "] with ";
@@ -1116 +1133 @@
     for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++)
       ost << *(*Runner) << " ";
-    ost << " at angle " << (a.ShortestAngle)<< ".";
+    ost << " at angle " << (a.ShortestAngle) << ".";
   }
 
   return ost;
-};
-
+}
+;
 
 // =========================================================== class TESSELATION ===========================================
@@ -1128 +1145 @@
  */
 Tesselation::Tesselation() :
-  PointsOnBoundaryCount(0),
-  LinesOnBoundaryCount(0),
-  TrianglesOnBoundaryCount(0),
-  LastTriangle(NULL),
-  TriangleFilesWritten(0),
-  InternalPointer(PointsOnBoundary.begin())
-{
-        Info FunctionInfo(__func__);
+  PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin())
+{
+  Info FunctionInfo(__func__);
 }
 ;
@@ -1144 +1156 @@
 Tesselation::~Tesselation()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl;
   for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
@@ -1151 +1163 @@
       runner->second = NULL;
     } else
-      DoeLog(1) && (eLog()<< Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl);
   }
   Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
@@ -1159 +1171 @@
 /** PointCloud implementation of GetCenter
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 Vector * Tesselation::GetCenter(ofstream *out) const
 {
-        Info FunctionInfo(__func__);
-  Vector *Center = new Vector(0.,0.,0.);
-  int num=0;
+  Info FunctionInfo(__func__);
+  Vector *Center = new Vector(0., 0., 0.);
+  int num = 0;
   for (GoToFirst(); (!IsEnd()); GoToNext()) {
     Center->AddVector(GetPoint()->node);
     num++;
   }
-  Center->Scale(1./num);
+  Center->Scale(1. / num);
   return Center;
-};
+}
+;
 
 /** PointCloud implementation of GoPoint
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 TesselPoint * Tesselation::GetPoint() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   return (InternalPointer->second->node);
-};
+}
+;
 
 /** PointCloud implementation of GetTerminalPoint.
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 TesselPoint * Tesselation::GetTerminalPoint() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   PointMap::const_iterator Runner = PointsOnBoundary.end();
   Runner--;
   return (Runner->second->node);
-};
+}
+;
 
 /** PointCloud implementation of GoToNext.
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 void Tesselation::GoToNext() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   if (InternalPointer != PointsOnBoundary.end())
     InternalPointer++;
-};
+}
+;
 
 /** PointCloud implementation of GoToPrevious.
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 void Tesselation::GoToPrevious() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   if (InternalPointer != PointsOnBoundary.begin())
     InternalPointer--;
-};
+}
+;
 
 /** PointCloud implementation of GoToFirst.
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 void Tesselation::GoToFirst() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   InternalPointer = PointsOnBoundary.begin();
-};
+}
+;
 
 /** PointCloud implementation of GoToLast.
@@ -1227 +1245 @@
 void Tesselation::GoToLast() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   InternalPointer = PointsOnBoundary.end();
   InternalPointer--;
-};
+}
+;
 
 /** PointCloud implementation of IsEmpty.
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 bool Tesselation::IsEmpty() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   return (PointsOnBoundary.empty());
-};
+}
+;
 
 /** PointCloud implementation of IsLast.
  * Uses PointsOnBoundary and STL stuff.
- */   
+ */
 bool Tesselation::IsEnd() const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   return (InternalPointer == PointsOnBoundary.end());
-};
-
+}
+;
 
 /** Gueses first starting triangle of the convex envelope.
@@ -1258 +1278 @@
 void Tesselation::GuessStartingTriangle()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // 4b. create a starting triangle
   // 4b1. create all distances
@@ -1268 +1288 @@
 
   // with A chosen, take each pair B,C and sort
-  if (A != PointsOnBoundary.end())
-    {
-      B = A;
-      B++;
-      for (; B != PointsOnBoundary.end(); B++)
-        {
-          C = B;
-          C++;
-          for (; C != PointsOnBoundary.end(); C++)
-            {
-              tmp = A->second->node->node->DistanceSquared(B->second->node->node);
-              distance = tmp * tmp;
-              tmp = A->second->node->node->DistanceSquared(C->second->node->node);
-              distance += tmp * tmp;
-              tmp = B->second->node->node->DistanceSquared(C->second->node->node);
-              distance += tmp * tmp;
-              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
-            }
-        }
-    }
+  if (A != PointsOnBoundary.end()) {
+    B = A;
+    B++;
+    for (; B != PointsOnBoundary.end(); B++) {
+      C = B;
+      C++;
+      for (; C != PointsOnBoundary.end(); C++) {
+        tmp = A->second->node->node->DistanceSquared(B->second->node->node);
+        distance = tmp * tmp;
+        tmp = A->second->node->node->DistanceSquared(C->second->node->node);
+        distance += tmp * tmp;
+        tmp = B->second->node->node->DistanceSquared(C->second->node->node);
+        distance += tmp * tmp;
+        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
+      }
+    }
+  }
   //    // listing distances
   //    Log() << Verbose(1) << "Listing DistanceMMap:";
@@ -1297 +1314 @@
   // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
   DistanceMultiMap::iterator baseline = DistanceMMap.begin();
-  for (; baseline != DistanceMMap.end(); baseline++)
-    {
-      // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
-      // 2. next, we have to check whether all points reside on only one side of the triangle
-      // 3. construct plane vector
-      PlaneVector.MakeNormalVector(A->second->node->node,
-          baseline->second.first->second->node->node,
-          baseline->second.second->second->node->node);
-      Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl;
-      // 4. loop over all points
-      double sign = 0.;
-      PointMap::iterator checker = PointsOnBoundary.begin();
-      for (; checker != PointsOnBoundary.end(); checker++)
-        {
-          // (neglecting A,B,C)
-          if ((checker == A) || (checker == baseline->second.first) || (checker
-              == baseline->second.second))
-            continue;
-          // 4a. project onto plane vector
-          TrialVector.CopyVector(checker->second->node->node);
-          TrialVector.SubtractVector(A->second->node->node);
-          distance = TrialVector.ScalarProduct(&PlaneVector);
-          if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
-            continue;
-          Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl;
-          tmp = distance / fabs(distance);
-          // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
-          if ((sign != 0) && (tmp != sign))
-            {
-              // 4c. If so, break 4. loop and continue with next candidate in 1. loop
-              Log() << Verbose(2) << "Current candidates: "
-                  << A->second->node->Name << ","
-                  << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leaves "
-                  << checker->second->node->Name << " outside the convex hull."
-                  << endl;
-              break;
-            }
-          else
-            { // note the sign for later
-              Log() << Verbose(2) << "Current candidates: "
-                  << A->second->node->Name << ","
-                  << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leave "
-                  << checker->second->node->Name << " inside the convex hull."
-                  << endl;
-              sign = tmp;
-            }
-          // 4d. Check whether the point is inside the triangle (check distance to each node
-          tmp = checker->second->node->node->DistanceSquared(A->second->node->node);
-          int innerpoint = 0;
-          if ((tmp < A->second->node->node->DistanceSquared(
-              baseline->second.first->second->node->node)) && (tmp
-              < A->second->node->node->DistanceSquared(
-                  baseline->second.second->second->node->node)))
-            innerpoint++;
-          tmp = checker->second->node->node->DistanceSquared(
-              baseline->second.first->second->node->node);
-          if ((tmp < baseline->second.first->second->node->node->DistanceSquared(
-              A->second->node->node)) && (tmp
-              < baseline->second.first->second->node->node->DistanceSquared(
-                  baseline->second.second->second->node->node)))
-            innerpoint++;
-          tmp = checker->second->node->node->DistanceSquared(
-              baseline->second.second->second->node->node);
-          if ((tmp < baseline->second.second->second->node->node->DistanceSquared(
-              baseline->second.first->second->node->node)) && (tmp
-              < baseline->second.second->second->node->node->DistanceSquared(
-                  A->second->node->node)))
-            innerpoint++;
-          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
-          if (innerpoint == 3)
-            break;
-        }
-      // 5. come this far, all on same side? Then break 1. loop and construct triangle
-      if (checker == PointsOnBoundary.end())
-        {
-          Log() << Verbose(2) << "Looks like we have a candidate!" << endl;
-          break;
-        }
-    }
-  if (baseline != DistanceMMap.end())
-    {
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = A->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = A->second;
-      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-
-      // 4b3. insert created triangle
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-      TrianglesOnBoundaryCount++;
-      for (int i = 0; i < NDIM; i++)
-        {
-          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
-          LinesOnBoundaryCount++;
-        }
-
-      Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
-    }
-  else
-    {
-      DoeLog(0) && (eLog()<< Verbose(0) << "No starting triangle found." << endl);
-    }
+  for (; baseline != DistanceMMap.end(); baseline++) {
+    // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+    // 2. next, we have to check whether all points reside on only one side of the triangle
+    // 3. construct plane vector
+    PlaneVector.MakeNormalVector(A->second->node->node, baseline->second.first->second->node->node, baseline->second.second->second->node->node);
+    Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl;
+    // 4. loop over all points
+    double sign = 0.;
+    PointMap::iterator checker = PointsOnBoundary.begin();
+    for (; checker != PointsOnBoundary.end(); checker++) {
+      // (neglecting A,B,C)
+      if ((checker == A) || (checker == baseline->second.first) || (checker == baseline->second.second))
+        continue;
+      // 4a. project onto plane vector
+      TrialVector.CopyVector(checker->second->node->node);
+      TrialVector.SubtractVector(A->second->node->node);
+      distance = TrialVector.ScalarProduct(&PlaneVector);
+      if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
+        continue;
+      Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl;
+      tmp = distance / fabs(distance);
+      // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
+      if ((sign != 0) && (tmp != sign)) {
+        // 4c. If so, break 4. loop and continue with next candidate in 1. loop
+        Log() << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name << "," << baseline->second.second->second->node->Name << " leaves " << checker->second->node->Name << " outside the convex hull." << endl;
+        break;
+      } else { // note the sign for later
+        Log() << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name << "," << baseline->second.second->second->node->Name << " leave " << checker->second->node->Name << " inside the convex hull." << endl;
+        sign = tmp;
+      }
+      // 4d. Check whether the point is inside the triangle (check distance to each node
+      tmp = checker->second->node->node->DistanceSquared(A->second->node->node);
+      int innerpoint = 0;
+      if ((tmp < A->second->node->node->DistanceSquared(baseline->second.first->second->node->node)) && (tmp < A->second->node->node->DistanceSquared(baseline->second.second->second->node->node)))
+        innerpoint++;
+      tmp = checker->second->node->node->DistanceSquared(baseline->second.first->second->node->node);
+      if ((tmp < baseline->second.first->second->node->node->DistanceSquared(A->second->node->node)) && (tmp < baseline->second.first->second->node->node->DistanceSquared(baseline->second.second->second->node->node)))
+        innerpoint++;
+      tmp = checker->second->node->node->DistanceSquared(baseline->second.second->second->node->node);
+      if ((tmp < baseline->second.second->second->node->node->DistanceSquared(baseline->second.first->second->node->node)) && (tmp < baseline->second.second->second->node->node->DistanceSquared(A->second->node->node)))
+        innerpoint++;
+      // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+      if (innerpoint == 3)
+        break;
+    }
+    // 5. come this far, all on same side? Then break 1. loop and construct triangle
+    if (checker == PointsOnBoundary.end()) {
+      Log() << Verbose(2) << "Looks like we have a candidate!" << endl;
+      break;
+    }
+  }
+  if (baseline != DistanceMMap.end()) {
+    BPS[0] = baseline->second.first->second;
+    BPS[1] = baseline->second.second->second;
+    BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+    BPS[0] = A->second;
+    BPS[1] = baseline->second.second->second;
+    BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+    BPS[0] = baseline->second.first->second;
+    BPS[1] = A->second;
+    BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+
+    // 4b3. insert created triangle
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+    TrianglesOnBoundaryCount++;
+    for (int i = 0; i < NDIM; i++) {
+      LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
+      LinesOnBoundaryCount++;
+    }
+
+    Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+  } else {
+    DoeLog(0) && (eLog() << Verbose(0) << "No starting triangle found." << endl);
+  }
 }
 ;
@@ -1424 +1407 @@
 void Tesselation::TesselateOnBoundary(const PointCloud * const cloud)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   bool flag;
   PointMap::iterator winner;
@@ -1493 +1476 @@
             TempAngle = VirtualNormalVector.Angle(&PropagationVector);
             Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
-            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
+            if (TempAngle > (M_PI / 2.)) { // no bends bigger than Pi/2 (90 degrees)
               Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
               continue;
@@ -1534 +1517 @@
             TempVector.AddVector(baseline->second->endpoints[1]->node->node);
             TempVector.AddVector(target->second->node->node);
-            TempVector.Scale(1./3.);
+            TempVector.Scale(1. / 3.);
             TempVector.SubtractVector(Center);
             // make it always point outward
@@ -1603 +1586 @@
           TrianglesOnBoundaryCount++;
         } else {
-          DoeLog(2) && (eLog()<< Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl);
+          DoeLog(2) && (eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl);
         }
 
@@ -1612 +1595 @@
 
   // exit
-  delete(Center);
-};
+  delete (Center);
+}
+;
 
 /** Inserts all points outside of the tesselated surface into it by adding new triangles.
@@ -1623 +1607 @@
 bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Vector Intersection, Normal;
   TesselPoint *Walker = NULL;
@@ -1633 +1617 @@
   cloud->GoToFirst();
   BoundaryPoints = new LinkedCell(this, 5.);
-  while (!cloud->IsEnd()) {  // we only have to go once through all points, as boundary can become only bigger
+  while (!cloud->IsEnd()) { // we only have to go once through all points, as boundary can become only bigger
     if (AddFlag) {
-      delete(BoundaryPoints);
+      delete (BoundaryPoints);
       BoundaryPoints = new LinkedCell(this, 5.);
       AddFlag = false;
@@ -1664 +1648 @@
         class BoundaryPointSet *OldPoints[3], *NewPoint;
         // store the three old lines and old points
-        for (int i=0;i<3;i++) {
+        for (int i = 0; i < 3; i++) {
           OldLines[i] = BTS->lines[i];
           OldPoints[i] = BTS->endpoints[i];
@@ -1672 +1656 @@
         Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl;
         AddFlag = true;
-        if (AddBoundaryPoint(Walker,0))
+        if (AddBoundaryPoint(Walker, 0))
           NewPoint = BPS[0];
         else
@@ -1679 +1663 @@
         Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl;
         TrianglesOnBoundary.erase(BTS->Nr);
-        delete(BTS);
+        delete (BTS);
         // create three new boundary lines
-        for (int i=0;i<3;i++) {
+        for (int i = 0; i < 3; i++) {
           BPS[0] = NewPoint;
           BPS[1] = OldPoints[i];
@@ -1690 +1674 @@
         }
         // create three new triangle with new point
-        for (int i=0;i<3;i++) { // find all baselines
+        for (int i = 0; i < 3; i++) { // find all baselines
           BLS[0] = OldLines[i];
           int n = 1;
-          for (int j=0;j<3;j++) {
+          for (int j = 0; j < 3; j++) {
             if (NewLines[j]->IsConnectedTo(BLS[0])) {
-              if (n>2) {
-                DoeLog(2) && (eLog()<< Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl);
+              if (n > 2) {
+                DoeLog(2) && (eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl);
                 return false;
               } else
@@ -1713 +1697 @@
       }
     } else { // something is wrong with FindClosestTriangleToPoint!
-      DoeLog(1) && (eLog()<< Verbose(1) << "The closest triangle did not produce an intersection!" << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl);
       return false;
     }
@@ -1720 +1704 @@
 
   // exit
-  delete(Center);
+  delete (Center);
   return true;
-};
+}
+;
 
 /** Adds a point to the tesselation::PointsOnBoundary list.
@@ -1731 +1716 @@
 bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   PointTestPair InsertUnique;
   BPS[n] = new class BoundaryPointSet(Walker);
@@ -1739 +1724 @@
     return true;
   } else {
-    delete(BPS[n]);
+    delete (BPS[n]);
     BPS[n] = InsertUnique.first->second;
     return false;
@@ -1753 +1738 @@
 void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   PointTestPair InsertUnique;
   TPS[n] = new class BoundaryPointSet(Candidate);
@@ -1774 +1759 @@
 void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr);
   if (FindPoint != PointsOnBoundary.end())
@@ -1780 +1765 @@
   else
     TPS[n] = NULL;
-};
+}
+;
 
 /** Function tries to add line from current Points in BPS to BoundaryLineSet.
@@ -1791 +1777 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) {
+void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
+{
   bool insertNewLine = true;
-
   LineMap::iterator FindLine = a->lines.find(b->node->nr);
   BoundaryLineSet *WinningLine = NULL;
@@ -1799 +1785 @@
     Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
 
-    pair<LineMap::iterator,LineMap::iterator> FindPair;
+    pair<LineMap::iterator, LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
 
-    for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
+    for (FindLine = FindPair.first; (FindLine != FindPair.second) && (insertNewLine); FindLine++) {
       Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl;
       // If there is a line with less than two attached triangles, we don't need a new line.
@@ -1812 +1798 @@
           Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate." << endl;
         // get open line
-        if ((!Finder->second->pointlist.empty()) && (*(Finder->second->pointlist.begin()) == candidate->node)
-            && (OptCenter == NULL || *OptCenter == Finder->second->OptCenter)) { // stop searching if candidate matches
-          insertNewLine = false;
-          WinningLine = FindLine->second;
-          break;
+        for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
+          if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || *OptCenter == Finder->second->OptCenter)) { // stop searching if candidate matches
+            insertNewLine = false;
+            WinningLine = FindLine->second;
+            break;
+          }
         }
       }
@@ -1840 +1827 @@
 void Tesselation::AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
-  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
+  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); // this also adds the line to the local maps
   // add line to global map
   LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
@@ -1851 +1838 @@
   // also add to open lines
   CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
-  OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
-};
+  OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
+}
+;
 
 /** Uses an existing line for a new triangle.
@@ -1868 +1856 @@
   BPS[1] = Line->endpoints[1];
   BLS[n] = Line;
-
   // 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);
+    delete (CandidateLine->second);
     OpenLines.erase(CandidateLine);
   } else {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
-  }
-};
+    DoeLog(1) && (eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
+  }
+}
+;
 
 /** Function adds triangle to global list.
@@ -1885 +1873 @@
 void Tesselation::AddTesselationTriangle()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
 
@@ -1896 +1884 @@
 
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
-};
+}
+;
 
 /** Function adds triangle to global list.
@@ -1904 +1893 @@
 void Tesselation::AddTesselationTriangle(const int nr)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl;
 
@@ -1914 +1903 @@
 
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
-};
+}
+;
 
 /** Removes a triangle from the tesselation.
@@ -1923 +1913 @@
 void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   if (triangle == NULL)
     return;
@@ -1931 +1921 @@
       triangle->lines[i]->triangles.erase(triangle->Nr);
       if (triangle->lines[i]->triangles.empty()) {
-          Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
-          RemoveTesselationLine(triangle->lines[i]);
+        Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
+        RemoveTesselationLine(triangle->lines[i]);
       } else {
         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++)
+        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(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
-//          for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
-//            Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
-//          Log() << Verbose(0) << endl;
-//        }
+        //        for (int j=0;j<2;j++) {
+        //          Log() << Verbose(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
+        //          for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
+        //            Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
+        //          Log() << Verbose(0) << endl;
+        //        }
       }
-      triangle->lines[i] = NULL;  // free'd or not: disconnect
+      triangle->lines[i] = NULL; // free'd or not: disconnect
     } else
-      DoeLog(1) && (eLog()<< Verbose(1) << "This line " << i << " has already been free'd." << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl);
   }
 
   if (TrianglesOnBoundary.erase(triangle->Nr))
     Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl;
-  delete(triangle);
-};
+  delete (triangle);
+}
+;
 
 /** Removes a line from the tesselation.
@@ -1962 +1953 @@
 void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   int Numbers[2];
 
@@ -1996 +1987 @@
       } else {
         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++)
+        for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
           Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
         Log() << Verbose(0) << endl;
       }
-      line->endpoints[i] = NULL;  // free'd or not: disconnect
+      line->endpoints[i] = NULL; // free'd or not: disconnect
     } else
-      DoeLog(1) && (eLog()<< Verbose(1) << "Endpoint " << i << " has already been free'd." << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl);
   }
   if (!line->triangles.empty())
-    DoeLog(2) && (eLog()<< Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl);
+    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl);
 
   if (LinesOnBoundary.erase(line->Nr))
     Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl;
-  delete(line);
-};
+  delete (line);
+}
+;
 
 /** Removes a point from the tesselation.
@@ -2019 +2011 @@
 void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   if (point == NULL)
     return;
   if (PointsOnBoundary.erase(point->Nr))
     Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl;
-  delete(point);
-};
-
+  delete (point);
+}
+;
 
 /** Checks validity of a given sphere of a candidate line.
  * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
- * Note that endpoints are stored in Tesselation::TPS.
- * \param *OtherOptCenter center of the other triangle
+ * We check CandidateForTesselation::OtherOptCenter
+ * \param &CandidateLine contains other degenerated candidates which we have to subtract as well
  * \param RADIUS radius of sphere
  * \param *LC LinkedCell structure with other atoms
  * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
  */
-bool Tesselation::CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const
+bool Tesselation::CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) const
 {
   Info FunctionInfo(__func__);
@@ -2043 +2035 @@
   bool flag = true;
 
-  cout << Verbose(1) << "Check by: draw sphere {" << OtherOptCenter->x[0] << " " << OtherOptCenter->x[1] << " " << OtherOptCenter->x[2] << "} radius " << RADIUS << " resolution 30" << endl;
-
+  Log() << Verbose(1) << "Check by: draw sphere {" << CandidateLine.OtherOptCenter.x[0] << " " << CandidateLine.OtherOptCenter.x[1] << " " << CandidateLine.OtherOptCenter.x[2] << "} radius " << RADIUS << " resolution 30" << endl;
   // get all points inside the sphere
-  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, OtherOptCenter);
-
-  Log() << Verbose(1) << "The following atoms are inside sphere at " << *OtherOptCenter << ":" << endl;
+  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);
+
+  Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl;
   for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-    Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
+    Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl;
 
   // remove triangles's endpoints
-  for (int i=0;i<3;i++)
-    ListofPoints->remove(TPS[i]->node);
+  for (int i = 0; i < 2; i++)
+    ListofPoints->remove(CandidateLine.BaseLine->endpoints[i]->node);
+
+  // remove other candidates
+  for (TesselPointList::const_iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); ++Runner)
+    ListofPoints->remove(*Runner);
 
   // check for other points
@@ -2060 +2055 @@
     Log() << Verbose(1) << "CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere." << endl;
     flag = false;
-    Log() << Verbose(1) << "External atoms inside of sphere at " << *OtherOptCenter << ":" << endl;
+    Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl;
     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
-  }
-  delete(ListofPoints);
+      Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl;
+  }
+  delete (ListofPoints);
 
   return flag;
-};
-
+}
+;
 
 /** Checks whether the triangle consisting of the three points is already present.
@@ -2081 +2076 @@
 int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   int adjacentTriangleCount = 0;
   class BoundaryPointSet *Points[3];
@@ -2121 +2116 @@
   Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
   return adjacentTriangleCount;
-};
+}
+;
 
 /** Checks whether the triangle consisting of the three points is already present.
@@ -2133 +2129 @@
 class BoundaryTriangleSet * Tesselation::GetPresentTriangle(TesselPoint *Candidates[3])
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   class BoundaryTriangleSet *triangle = NULL;
   class BoundaryPointSet *Points[3];
@@ -2171 +2167 @@
 
   return triangle;
-};
-
+}
+;
 
 /** Finds the starting triangle for FindNonConvexBorder().
@@ -2184 +2180 @@
 void Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   int i = 0;
   TesselPoint* MaxPoint[NDIM];
@@ -2193 +2189 @@
   Vector Chord;
   Vector SearchDirection;
-  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  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;
@@ -2206 +2202 @@
 
   // 1. searching topmost point with respect to each axis
-  for (int i=0;i<NDIM;i++) { // each axis
-    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
-    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
-      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
+  for (int i = 0; i < NDIM; i++) { // each axis
+    LC->n[i] = LC->N[i] - 1; // current axis is topmost cell
+    for (LC->n[(i + 1) % NDIM] = 0; LC->n[(i + 1) % NDIM] < LC->N[(i + 1) % NDIM]; LC->n[(i + 1) % NDIM]++)
+      for (LC->n[(i + 2) % NDIM] = 0; LC->n[(i + 2) % NDIM] < LC->N[(i + 2) % NDIM]; LC->n[(i + 2) % NDIM]++) {
         const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
         //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
         if (List != NULL) {
-          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) {
+          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             if ((*Runner)->node->x[i] > maxCoordinate[i]) {
               Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
@@ -2221 +2217 @@
           }
         } else {
-          DoeLog(1) && (eLog()<< Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
+          DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
         }
       }
@@ -2227 +2223 @@
 
   Log() << Verbose(1) << "Found maximum coordinates: ";
-  for (int i=0;i<NDIM;i++)
+  for (int i = 0; i < NDIM; i++)
     Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t";
   Log() << Verbose(0) << endl;
 
   BTS = NULL;
-  for (int k=0;k<NDIM;k++) {
+  for (int k = 0; k < NDIM; k++) {
     NormalVector.Zero();
     NormalVector.x[k] = 1.;
@@ -2260 +2256 @@
 
     double radius = CirclePlaneNormal.NormSquared();
-    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+    double CircleRadius = sqrt(RADIUS * RADIUS - radius / 4.);
 
     NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
     NormalVector.Normalize();
-    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+    ShortestAngle = 2. * M_PI; // This will indicate the quadrant.
 
     SphereCenter.CopyVector(&NormalVector);
@@ -2272 +2268 @@
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection.MakeNormalVector(&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
@@ -2283 +2279 @@
     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;
+      Log() << Verbose(0) << " " << *(*it) << endl;
     }
     if (!OptCandidates.pointlist.empty()) {
@@ -2302 +2298 @@
     delete BaseLine;
   }
-};
+}
+;
 
 /** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates.
@@ -2453 +2450 @@
 bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
 {
-        Info FunctionInfo(__func__);
-  bool result = true;
-
+  Info FunctionInfo(__func__);
   Vector CircleCenter;
   Vector CirclePlaneNormal;
@@ -2465 +2460 @@
   double radius, CircleRadius;
 
-  for (int i=0;i<3;i++)
+  for (int i = 0; i < 3; i++)
     if ((T.endpoints[i] != CandidateLine.BaseLine->endpoints[0]) && (T.endpoints[i] != CandidateLine.BaseLine->endpoints[1])) {
       ThirdPoint = T.endpoints[i];
@@ -2485 +2480 @@
   // calculate squared radius of circle
   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
+  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.;
+    CircleRadius = RADIUS * RADIUS - radius / 4.;
     CirclePlaneNormal.Normalize();
     Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
@@ -2505 +2500 @@
     if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
-      DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
     }
 
@@ -2516 +2511 @@
 
   if (CandidateLine.pointlist.empty()) {
-    DoeLog(2) && (eLog()<< Verbose(2) << "Could not find a suitable candidate." << endl);
+    DoeLog(2) && (eLog() << Verbose(2) << "Could not find a suitable candidate." << endl);
     return false;
   }
@@ -2525 +2520 @@
 
   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 {
-////          DoeLog(1) && (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 {
-////            DoeLog(1) && (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)) {
-//      DoeLog(1) && (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;
-};
+}
+;
+
+/** Walks through Tesselation::OpenLines() and finds candidates for newly created ones.
+ * \param *&LCList atoms in LinkedCell list
+ * \param RADIUS radius of the virtual sphere
+ * \return true - for all open lines without candidates so far, a candidate has been found,
+ *         false - at least one open line without candidate still
+ */
+bool Tesselation::FindCandidatesforOpenLines(const double RADIUS, const LinkedCell *&LCList)
+{
+  bool TesselationFailFlag = true;
+  CandidateForTesselation *baseline = NULL;
+  BoundaryTriangleSet *T = NULL;
+
+  for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) {
+    baseline = Runner->second;
+    if (baseline->pointlist.empty()) {
+      T = (((baseline->BaseLine->triangles.begin()))->second);
+      Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
+      TesselationFailFlag = TesselationFailFlag && FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
+    }
+  }
+  return TesselationFailFlag;
+}
+;
 
 /** Adds the present line and candidate point from \a &CandidateLine to the Tesselation.
@@ -2641 +2577 @@
   Sprinter = Runner;
   Sprinter++;
-  while(Sprinter != connectedClosestPoints->end()) {
+  while (Sprinter != connectedClosestPoints->end()) {
+    Log() << Verbose(0) << "Current Runner is " << *(*Runner) << " and sprinter is " << *(*Sprinter) << "." << endl;
+
     AddTesselationPoint(TurningPoint, 0);
     AddTesselationPoint(*Runner, 1);
     AddTesselationPoint(*Sprinter, 2);
 
+    AddCandidateTriangle(CandidateLine, Opt);
+
+    Runner = Sprinter;
+    Sprinter++;
+    if (Sprinter != connectedClosestPoints->end()) {
+      // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
+      FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OptCenter);
+      Log() << Verbose(0) << " There are still more triangles to add." << endl;
+    }
+    // pick candidates for other open lines as well
+    FindCandidatesforOpenLines(RADIUS, LC);
+
     // check whether we add a degenerate or a normal triangle
-    if (CheckDegeneracy(&CandidateLine.OtherOptCenter, RADIUS, LC)) {
+    if (CheckDegeneracy(CandidateLine, RADIUS, LC)) {
       // add normal and degenerate triangles
       Log() << Verbose(1) << "Triangle of endpoints " << *TPS[0] << "," << *TPS[1] << " and " << *TPS[2] << " is degenerated, adding both sides." << endl;
-      AddDegeneratedTriangle(CandidateLine, RADIUS, LC);
-    } else {
-      // add this triangle
-      AddCandidateTriangle(CandidateLine);
-    }
-
-    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);
+      AddCandidateTriangle(CandidateLine, OtherOpt);
+
+      if (Sprinter != connectedClosestPoints->end()) {
+        // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
+        FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OtherOptCenter);
+      }
+      // pick candidates for other open lines as well
+      FindCandidatesforOpenLines(RADIUS, LC);
+    }
+  }
+  delete (connectedClosestPoints);
+};
+
+/** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate.
+ * \param *Sprinter next candidate to which internal open lines are set
+ * \param *OptCenter OptCenter for this candidate
+ */
+void Tesselation::FindDegeneratedCandidatesforOpenLines(TesselPoint * const Sprinter, const Vector * const OptCenter)
+{
+  Info FunctionInfo(__func__);
+
+  pair<LineMap::iterator, LineMap::iterator> FindPair = TPS[0]->lines.equal_range(TPS[2]->node->nr);
+  for (LineMap::const_iterator FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
+    Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl;
+    // If there is a line with less than two attached triangles, we don't need a new line.
+    if (FindLine->second->triangles.size() == 1) {
+      CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
+      if (!Finder->second->pointlist.empty())
+        Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl;
+      else {
+        Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter) << endl;
+        Finder->second->pointlist.push_back(Sprinter);
+        Finder->second->ShortestAngle = 0.;
+        Finder->second->OptCenter.CopyVector(OptCenter);
+      }
+    }
+  }
 };
 
@@ -2680 +2654 @@
 
   /// 1. Create or pick the lines for the first triangle
-  // for each amount of open lines we have a different case:
-  // case 0: no triangles at this line and not closed
-  // case 1: no triangles at new line is closed
-  // case 2: one line with one triangle attached, one will used this line
-  // case 3: one line with two triangles attached, new line will be used by both degenerate triangles
-  // case 4: two lines with one triangle each, each new triangle will pick one
   Log() << Verbose(0) << "INFO: Creating/Picking lines for first triangle ..." << endl;
-  for (int i=0;i<3;i++) {
+  for (int i = 0; i < 3; i++) {
     BLS[i] = NULL;
-    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
-    AddTesselationLine(&CandidateLine.OptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl;
+    AddTesselationLine(&CandidateLine.OptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
   }
 
@@ -2712 +2680 @@
   /// 3. Gather candidates for each new line
   Log() << Verbose(0) << "INFO: Adding candidates to new lines ..." << endl;
-  for (int i=0;i<3;i++) {
-    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+  for (int i = 0; i < 3; i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl;
     CandidateCheck = OpenLines.find(BLS[i]);
     if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
@@ -2724 +2692 @@
   /// 4. Create or pick the lines for the second triangle
   Log() << Verbose(0) << "INFO: Creating/Picking lines for second triangle ..." << endl;
-  for (int i=0;i<3;i++) {
-    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
-    AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+  for (int i = 0; i < 3; i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl;
+    AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
   }
 
@@ -2746 +2714 @@
   /// 6. Adding triangle to new lines
   Log() << Verbose(0) << "INFO: Adding second triangles to new lines ..." << endl;
-  for (int i=0;i<3;i++) {
-    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+  for (int i = 0; i < 3; i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl;
     CandidateCheck = OpenLines.find(BLS[i]);
     if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
@@ -2754 +2722 @@
     }
   }
-};
+}
+;
 
 /** Adds a triangle to the Tesselation structure from three given TesselPoint's.
  * Note that endpoints are in Tesselation::TPS.
  * \param CandidateLine CandidateForTesselation structure contains other information
- */
-void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine)
+ * \param type which opt center to add (i.e. which side) and thus which NormalVector to take
+ */
+void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine, enum centers type)
 {
   Info FunctionInfo(__func__);
   Vector Center;
+  Vector *OptCenter = (type == Opt) ? &CandidateLine.OptCenter : &CandidateLine.OtherOptCenter;
 
   // add the lines
-  AddTesselationLine(&CandidateLine.OptCenter, TPS[2], TPS[0], TPS[1], 0);
-  AddTesselationLine(&CandidateLine.OptCenter, TPS[1], TPS[0], TPS[2], 1);
-  AddTesselationLine(&CandidateLine.OptCenter, TPS[0], TPS[1], TPS[2], 2);
+  AddTesselationLine(OptCenter, TPS[2], TPS[0], TPS[1], 0);
+  AddTesselationLine(OptCenter, TPS[1], TPS[0], TPS[2], 1);
+  AddTesselationLine(OptCenter, TPS[0], TPS[1], TPS[2], 2);
 
   // add the triangles
@@ -2776 +2747 @@
   // create normal vector
   BTS->GetCenter(&Center);
-  Center.SubtractVector(&CandidateLine.OptCenter);
-  BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+  Center.SubtractVector(OptCenter);
+  BTS->SphereCenter.CopyVector(OptCenter);
   BTS->GetNormalVector(Center);
 
   // give some verbose output about the whole procedure
   if (CandidateLine.T != NULL)
-    Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
+    Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
   else
-    Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
-};
+    Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
+}
+;
 
 /** Checks whether the quadragon of the two triangles connect to \a *Base is convex.
@@ -2796 +2768 @@
 class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   class BoundaryPointSet *Spot = NULL;
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint;
 
-  int m=0;
-  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j=0;j<3;j++) // all of their endpoints and baselines
+  int m = 0;
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
       if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
         BPS[m++] = runner->second->endpoints[j];
-  OtherBase = new class BoundaryLineSet(BPS,-1);
+  OtherBase = new class BoundaryLineSet(BPS, -1);
 
   Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl;
@@ -2815 +2787 @@
 
   // delete the temporary other base line
-  delete(OtherBase);
+  delete (OtherBase);
 
   // get the distance vector from Base line to OtherBase line
@@ -2822 +2794 @@
   BaseLine.CopyVector(Base->endpoints[1]->node->node);
   BaseLine.SubtractVector(Base->endpoints[0]->node->node);
-  for (int i=0;i<2;i++) {
+  for (int i = 0; i < 2; i++) {
     DistanceToIntersection[i].CopyVector(ClosestPoint);
     DistanceToIntersection[i].SubtractVector(Base->endpoints[i]->node->node);
     distance[i] = BaseLine.ScalarProduct(&DistanceToIntersection[i]);
   }
-  delete(ClosestPoint);
-  if ((distance[0] * distance[1]) > 0)  { // have same sign?
-    Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1]  << ". " << *Base << "' rectangle is concave." << endl;
+  delete (ClosestPoint);
+  if ((distance[0] * distance[1]) > 0) { // have same sign?
+    Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl;
     if (distance[0] < distance[1]) {
       Spot = Base->endpoints[0];
@@ -2836 +2808 @@
     }
     return Spot;
-  } else {  // different sign, i.e. we are in between
+  } else { // different sign, i.e. we are in between
     Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;
     return NULL;
   }
 
-};
+}
+;
 
 void Tesselation::PrintAllBoundaryPoints(ofstream *out) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // print all lines
   Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl;
-  for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++)
+  for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin(); PointRunner != PointsOnBoundary.end(); PointRunner++)
     Log() << Verbose(0) << *(PointRunner->second) << endl;
-};
+}
+;
 
 void Tesselation::PrintAllBoundaryLines(ofstream *out) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // print all lines
   Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl;
   for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
     Log() << Verbose(0) << *(LineRunner->second) << endl;
-};
+}
+;
 
 void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // print all triangles
   Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl;
   for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
     Log() << Verbose(0) << *(TriangleRunner->second) << endl;
-};
+}
+;
 
 /** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher".
@@ -2877 +2853 @@
 double Tesselation::PickFarthestofTwoBaselines(class BoundaryLineSet *Base)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint[2];
   double volume;
 
-  int m=0;
-  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j=0;j<3;j++) // all of their endpoints and baselines
+  int m = 0;
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
       if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
         BPS[m++] = runner->second->endpoints[j];
-  OtherBase = new class BoundaryLineSet(BPS,-1);
+  OtherBase = new class BoundaryLineSet(BPS, -1);
 
   Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl;
@@ -2905 +2881 @@
 
   // delete the temporary other base line and the closest points
-  delete(ClosestPoint[0]);
-  delete(ClosestPoint[1]);
-  delete(OtherBase);
+  delete (ClosestPoint[0]);
+  delete (ClosestPoint[1]);
+  delete (OtherBase);
 
   if (Distance.NormSquared() < MYEPSILON) { // check for intersection
@@ -2916 +2892 @@
     BaseLineNormal.Zero();
     if (Base->triangles.size() < 2) {
-      DoeLog(1) && (eLog()<< Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
       return 0.;
     }
@@ -2923 +2899 @@
       BaseLineNormal.AddVector(&(runner->second->NormalVector));
     }
-    BaseLineNormal.Scale(1./2.);
+    BaseLineNormal.Scale(1. / 2.);
 
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
@@ -2929 +2905 @@
       // calculate volume summand as a general tetraeder
       return volume;
-    } else {  // Base higher than OtherBase -> do nothing
+    } else { // Base higher than OtherBase -> do nothing
       Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl;
       return 0.;
     }
   }
-};
+}
+;
 
 /** For a given baseline and its two connected triangles, flips the baseline.
@@ -2945 +2922 @@
 class BoundaryLineSet * Tesselation::FlipBaseline(class BoundaryLineSet *Base)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
   Vector BaseLineNormal;
   int OldTriangleNrs[2], OldBaseLineNr;
-  int i,m;
+  int i, m;
 
   // calculate NormalVector for later use
   BaseLineNormal.Zero();
   if (Base->triangles.size() < 2) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
     return NULL;
   }
@@ -2962 +2939 @@
     BaseLineNormal.AddVector(&(runner->second->NormalVector));
   }
-  BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
+  BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
 
   // get the two triangles
   // gather four endpoints and four lines
-  for (int j=0;j<4;j++)
+  for (int j = 0; j < 4; j++)
     OldLines[j] = NULL;
-  for (int j=0;j<2;j++)
+  for (int j = 0; j < 2; j++)
     OldPoints[j] = NULL;
-  i=0;
-  m=0;
+  i = 0;
+  m = 0;
   Log() << Verbose(0) << "The four old lines are: ";
-  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j=0;j<3;j++) // all of their endpoints and baselines
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
       if (runner->second->lines[j] != Base) { // pick not the central baseline
         OldLines[i++] = runner->second->lines[j];
@@ -2981 +2958 @@
   Log() << Verbose(0) << endl;
   Log() << Verbose(0) << "The two old points are: ";
-  for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-    for (int j=0;j<3;j++) // all of their endpoints and baselines
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
       if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints
         OldPoints[m++] = runner->second->endpoints[j];
@@ -2990 +2967 @@
 
   // check whether everything is in place to create new lines and triangles
-  if (i<4) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "We have not gathered enough baselines!" << endl);
+  if (i < 4) {
+    DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
     return NULL;
   }
-  for (int j=0;j<4;j++)
+  for (int j = 0; j < 4; j++)
     if (OldLines[j] == NULL) {
-      DoeLog(1) && (eLog()<< Verbose(1) << "We have not gathered enough baselines!" << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
       return NULL;
     }
-  for (int j=0;j<2;j++)
+  for (int j = 0; j < 2; j++)
     if (OldPoints[j] == NULL) {
-      DoeLog(1) && (eLog()<< Verbose(1) << "We have not gathered enough endpoints!" << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl);
       return NULL;
     }
@@ -3008 +2985 @@
   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++) {
+  m = 0;
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
     Log() << Verbose(0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;
     OldTriangleNrs[m++] = runner->second->Nr;
@@ -3023 +3000 @@
 
   // construct new triangles with flipped baseline
-  i=-1;
+  i = -1;
   if (OldLines[0]->IsConnectedTo(OldLines[2]))
-    i=2;
+    i = 2;
   if (OldLines[0]->IsConnectedTo(OldLines[3]))
-    i=3;
-  if (i!=-1) {
+    i = 3;
+  if (i != -1) {
     BLS[0] = OldLines[0];
     BLS[1] = OldLines[i];
@@ -3037 +3014 @@
     Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl;
 
-    BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
+    BLS[0] = (i == 2 ? OldLines[3] : OldLines[2]);
     BLS[1] = OldLines[1];
     BLS[2] = NewLine;
@@ -3045 +3022 @@
     Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl;
   } else {
-    DoeLog(0) && (eLog()<< Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl);
+    DoeLog(0) && (eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl);
     return NULL;
   }
 
   return NewLine;
-};
-
+}
+;
 
 /** Finds the second point of starting triangle.
@@ -3063 +3040 @@
 void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   Vector AngleCheck;
   class TesselPoint* Candidate = NULL;
@@ -3072 +3049 @@
   int Nupper[NDIM];
 
-  if (LC->SetIndexToNode(a)) {  // get cell for the starting point
-    for(int i=0;i<NDIM;i++) // store indices of this cell
+  if (LC->SetIndexToNode(a)) { // get cell for the starting point
+    for (int i = 0; i < NDIM; i++) // store indices of this cell
       N[i] = LC->n[i];
   } else {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl);
     return;
   }
   // then go through the current and all neighbouring cells and check the contained points for possible candidates
-  for (int i=0;i<NDIM;i++) {
-    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-  }
-  Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :"
-    << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl;
+  for (int i = 0; i < NDIM; i++) {
+    Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
+    Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
+  }
+  Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl;
 
   for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
@@ -3119 +3095 @@
               norm = aCandidate.Norm();
               // second point shall have smallest angle with respect to Oben vector
-              if (norm < RADIUS*2.) {
+              if (norm < RADIUS * 2.) {
                 angle = AngleCheck.Angle(&Oben);
                 if (angle < Storage[0]) {
@@ -3141 +3117 @@
         }
       }
-};
-
+}
+;
 
 /** This recursive function finds a third point, to form a triangle with two given ones.
@@ -3174 +3150 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
-{
-        Info FunctionInfo(__func__);
-  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+void Tesselation::FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
+{
+  Info FunctionInfo(__func__);
+  Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
   Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
   Vector SphereCenter;
-  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
-  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
-  Vector NewNormalVector;   // normal vector of the Candidate's triangle
+  Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+  Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+  Vector NewNormalVector; // normal vector of the Candidate's triangle
   Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
   Vector RelativeOldSphereCenter;
@@ -3213 +3189 @@
 
   // calculate squared radius TesselPoint *ThirdPoint,f circle
-  radius = CirclePlaneNormal.NormSquared()/4.;
-  if (radius < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius;
+  radius = CirclePlaneNormal.NormSquared() / 4.;
+  if (radius < RADIUS * RADIUS) {
+    CircleRadius = RADIUS * RADIUS - radius;
     CirclePlaneNormal.Normalize();
     Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
@@ -3221 +3197 @@
     // test whether old center is on the band's plane
     if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      DoeLog(1) && (eLog()<< Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl);
       RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
     }
@@ -3230 +3206 @@
       // check SearchDirection
       Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
-        DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
+        DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);
       }
 
       // get cell for the starting point
       if (LC->SetIndexToVector(&CircleCenter)) {
-        for(int i=0;i<NDIM;i++) // store indices of this cell
-        N[i] = LC->n[i];
+        for (int i = 0; i < NDIM; i++) // store indices of this cell
+          N[i] = LC->n[i];
         //Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
       } else {
-        DoeLog(1) && (eLog()<< Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl);
+        DoeLog(1) && (eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl);
         return;
       }
       // then go through the current and all neighbouring cells and check the contained points for possible candidates
       //Log() << Verbose(1) << "LC Intervals:";
-      for (int i=0;i<NDIM;i++) {
-        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+      for (int i = 0; i < NDIM; i++) {
+        Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
+        Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
         //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
       }
@@ -3262 +3238 @@
                 // check for three unique points
                 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) ){
+                if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node)) {
 
                   // find center on the plane
@@ -3268 +3244 @@
                   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)
-                  ) {
+                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node) && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)) {
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
                     radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
@@ -3276 +3250 @@
                     Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
                     Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl;
-                    if (radius < RADIUS*RADIUS) {
+                    if (radius < RADIUS * RADIUS) {
                       otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
                       if (fabs(radius - otherradius) > HULLEPSILON) {
-                        DoeLog(1) && (eLog()<< Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl);
+                        DoeLog(1) && (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
@@ -3285 +3259 @@
                       OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
                       helper.CopyVector(&NewNormalVector);
-                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                      helper.Scale(sqrt(RADIUS * RADIUS - radius));
                       Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl;
                       NewSphereCenter.AddVector(&helper);
@@ -3311 +3285 @@
                         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;
+                          Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
                         } else {
                           // remove all candidates from the list and then the list itself
                           CandidateLine.pointlist.clear();
                           CandidateLine.pointlist.push_back(Candidate);
-                          Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with "
-                            << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
+                          Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;
                         }
                         CandidateLine.ShortestAngle = alpha;
@@ -3346 +3318 @@
           }
     } else {
-      DoeLog(1) && (eLog()<< Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
     }
   } else {
@@ -3355 +3327 @@
   }
 
-  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
-    DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl);
-    performCriticalExit();
-  }
-
   DoLog(1) && (Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl);
   if (CandidateLine.pointlist.size() > 1) {
@@ -3365 +3332 @@
     CandidateLine.pointlist.sort(); //SortCandidates);
   }
-};
+
+  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
+    DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl);
+    performCriticalExit();
+  }
+}
+;
 
 /** Finds the endpoint two lines are sharing.
@@ -3374 +3347 @@
 class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
@@ -3381 +3354 @@
   for (int i = 0; i < 2; i++)
     // for both lines
-    for (int j = 0; j < 2; j++)
-      { // for both endpoints
-        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
-            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
-        if (!OrderTest.second)
-          { // if insertion fails, we have common endpoint
-            node = OrderTest.first->second;
-            Log() << Verbose(1) << "Common endpoint of lines " << *line1
-                << " and " << *line2 << " is: " << *node << "." << endl;
-            j = 2;
-            i = 2;
-            break;
-          }
+    for (int j = 0; j < 2; j++) { // for both endpoints
+      OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
+      if (!OrderTest.second) { // if insertion fails, we have common endpoint
+        node = OrderTest.first->second;
+        Log() << Verbose(1) << "Common endpoint of lines " << *line1 << " and " << *line2 << " is: " << *node << "." << endl;
+        j = 2;
+        i = 2;
+        break;
       }
+    }
   return node;
-};
+}
+;
 
 /** Finds the boundary points that are closest to a given Vector \a *x.
@@ -3410 +3380 @@
 
   if (LinesOnBoundary.empty()) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl);
     return NULL;
   }
@@ -3416 +3386 @@
   // 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
+  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);
@@ -3432 +3401 @@
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end()) {
-              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+              points->insert(DistanceToPointPair(FindPoint->second->node->node->DistanceSquared(x), FindPoint->second));
               Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
             }
           }
         } else {
-          DoeLog(1) && (eLog()<< Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
+          DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
         }
       }
@@ -3443 +3412 @@
   // check whether we found some points
   if (points->empty()) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
-    delete(points);
+    DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
+    delete (points);
     return NULL;
   }
   return points;
-};
+}
+;
 
 /** Finds the boundary line that is closest to a given Vector \a *x.
@@ -3458 +3428 @@
 {
   Info FunctionInfo(__func__);
-
   // get closest points
-  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x, LC);
   if (points == NULL) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
     return NULL;
   }
@@ -3495 +3464 @@
         helper.SubtractVector(&Center);
         const double lengthB = helper.ScalarProduct(&BaseLine);
-        if (lengthB*lengthA < 0) {  // if have different sign
+        if (lengthB * lengthA < 0) { // if have different sign
           ClosestLine = LineRunner->second;
           MinDistance = distance;
@@ -3507 +3476 @@
     }
   }
-  delete(points);
+  delete (points);
   // check whether closest line is "too close" :), then it's inside
   if (ClosestLine == NULL) {
@@ -3514 +3483 @@
   }
   return ClosestLine;
-};
+}
+;
 
 /** Finds the triangle that is closest to a given Vector \a *x.
@@ -3523 +3493 @@
 TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
 {
-        Info FunctionInfo(__func__);
-
-        // get closest points
-        DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  Info FunctionInfo(__func__);
+  // get closest points
+  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x, LC);
   if (points == NULL) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
     return NULL;
   }
@@ -3555 +3524 @@
       const double lengthEndB = BaseLineIntersection.NormSquared();
 
-      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
+      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
@@ -3562 +3531 @@
           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
+        } 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;
@@ -3577 +3546 @@
         const double distance = BaseLineIntersection.NormSquared();
         if (Center.NormSquared() > BaseLine.NormSquared()) {
-          DoeLog(0) && (eLog()<< Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl);
+          DoeLog(0) && (eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl);
         }
         if ((ClosestLines.empty()) || (distance < MinDistance)) {
@@ -3589 +3558 @@
     }
   }
-  delete(points);
+  delete (points);
 
   // check whether closest line is "too close" :), then it's inside
@@ -3599 +3568 @@
   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);
-  }
+      candidates->push_back(Runner->second);
+    }
   return candidates;
-};
+}
+;
 
 /** Finds closest triangle to a point.
@@ -3613 +3583 @@
 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
   TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
@@ -3624 +3594 @@
 
   // go through all and pick the one with the best alignment to x
-  double MinAlignment = 2.*M_PI;
+  double MinAlignment = 2. * M_PI;
   for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
     (*Runner)->GetCenter(&Center);
@@ -3638 +3608 @@
     }
   }
-  delete(triangles);
+  delete (triangles);
 
   return result;
-};
+}
+;
 
 /** Checks whether the provided Vector is within the Tesselation structure.
@@ -3653 +3624 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point,this,LC);
+  TriangleIntersectionList Intersections(&Point, this, LC);
 
   return Intersections.IsInside();
-};
+}
+;
 
 /** Returns the distance to the surface given by the tesselation.
@@ -3727 +3699 @@
     }
   }
-};
+}
+;
 
 /** Calculates minimum distance from \a&Point to a tesselated surface.
@@ -3738 +3711 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point,this,LC);
+  TriangleIntersectionList Intersections(&Point, this, LC);
 
   return Intersections.GetSmallestDistance();
-};
+}
+;
 
 /** Calculates minimum distance from \a&Point to a tesselated surface.
@@ -3752 +3726 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point,this,LC);
+  TriangleIntersectionList Intersections(&Point, this, LC);
 
   return Intersections.GetClosestTriangle();
-};
+}
+;
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -3765 +3740 @@
 TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
 {
-        Info FunctionInfo(__func__);
-        TesselPointSet *connectedPoints = new TesselPointSet;
+  Info FunctionInfo(__func__);
+  TesselPointSet *connectedPoints = new TesselPointSet;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
   bool takePoint = false;
-
   // find the respective boundary point
   PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
@@ -3776 +3750 @@
     ReferencePoint = PointRunner->second;
   } else {
-    DoeLog(2) && (eLog()<< Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
+    DoeLog(2) && (eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
     ReferencePoint = NULL;
   }
@@ -3782 +3756 @@
   // little trick so that we look just through lines connect to the BoundaryPoint
   // OR fall-back to look through all lines if there is no such BoundaryPoint
-  const LineMap *Lines;;
+  const LineMap *Lines;
+  ;
   if (ReferencePoint != NULL)
     Lines = &(ReferencePoint->lines);
@@ -3789 +3764 @@
   LineMap::const_iterator findLines = Lines->begin();
   while (findLines != Lines->end()) {
-   takePoint = false;
-
-   if (findLines->second->endpoints[0]->Nr == Point->nr) {
-     takePoint = true;
-     current = findLines->second->endpoints[1]->node;
-   } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
-     takePoint = true;
-     current = findLines->second->endpoints[0]->node;
-   }
-
-   if (takePoint) {
-     Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
-     connectedPoints->insert(current);
-   }
-
-   findLines++;
+    takePoint = false;
+
+    if (findLines->second->endpoints[0]->Nr == Point->nr) {
+      takePoint = true;
+      current = findLines->second->endpoints[1]->node;
+    } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
+      takePoint = true;
+      current = findLines->second->endpoints[0]->node;
+    }
+
+    if (takePoint) {
+      Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
+      connectedPoints->insert(current);
+    }
+
+    findLines++;
   }
 
   if (connectedPoints->empty()) { // if have not found any points
-    DoeLog(1) && (eLog()<< Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "We have not found any connected points to " << *Point << "." << endl);
     return NULL;
   }
 
   return connectedPoints;
-};
-
+}
+;
 
 /** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
@@ -3830 +3805 @@
 TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
   TesselPointList *connectedCircle = new TesselPointList;
@@ -3837 +3812 @@
   Vector OrthogonalVector;
   Vector helper;
-  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
+  const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL };
   TriangleList *triangles = NULL;
 
   if (SetOfNeighbours == NULL) {
-    DoeLog(2) && (eLog()<< Verbose(2) << "Could not find any connected points!" << endl);
-    delete(connectedCircle);
+    DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);
+    delete (connectedCircle);
     return NULL;
   }
@@ -3852 +3827 @@
       PlaneNormal.AddVector(&(*Runner)->NormalVector);
   } else {
-    DoeLog(0) && (eLog()<< Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);
+    DoeLog(0) && (eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);
     performCriticalExit();
   }
-  PlaneNormal.Scale(1.0/triangles->size());
+  PlaneNormal.Scale(1.0 / triangles->size());
   Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl;
   PlaneNormal.Normalize();
@@ -3865 +3840 @@
     AngleZero.ProjectOntoPlane(&PlaneNormal);
   }
-  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
     Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;
     AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
@@ -3871 +3846 @@
     AngleZero.ProjectOntoPlane(&PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
-      DoeLog(0) && (eLog()<< Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
+      DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
       performCriticalExit();
     }
@@ -3889 +3864 @@
     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++) {
+    anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
+  }
+
+  for (map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
     connectedCircle->push_back(AngleRunner->second);
   }
@@ -3922 +3897 @@
 
   if (SetOfNeighbours == NULL) {
-    DoeLog(2) && (eLog()<< Verbose(2) << "Could not find any connected points!" << endl);
-    delete(connectedCircle);
+    DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);
+    delete (connectedCircle);
     return NULL;
   }
@@ -3936 +3911 @@
   Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
   // calculate central point
-
   TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
   TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
@@ -3955 +3929 @@
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
   //  << "; scale factor " << counter;
-  PlaneNormal.Scale(1.0/(double)counter);
-//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
-//
-//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
-//  PlaneNormal.CopyVector(Point->node);
-//  PlaneNormal.SubtractVector(&center);
-//  PlaneNormal.Normalize();
+  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;
 
@@ -3970 +3944 @@
     AngleZero.ProjectOntoPlane(&PlaneNormal);
   }
-  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
     Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;
     AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
@@ -3976 +3950 @@
     AngleZero.ProjectOntoPlane(&PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
-      DoeLog(0) && (eLog()<< Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
+      DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
       performCriticalExit();
     }
@@ -3988 +3962 @@
 
   // go through all connected points and calculate angle
-  pair <map<double, TesselPoint*>::iterator, bool > InserterTest;
+  pair<map<double, TesselPoint*>::iterator, bool> InserterTest;
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
@@ -3995 +3969 @@
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     if (angle > M_PI) // the correction is of no use here (and not desired)
-      angle = 2.*M_PI - angle;
+      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)));
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
     if (!InserterTest.second) {
-      DoeLog(0) && (eLog()<< Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl);
+      DoeLog(0) && (eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl);
       performCriticalExit();
     }
   }
 
-  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
+  for (map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
     connectedCircle->push_back(AngleRunner->second);
   }
@@ -4019 +3993 @@
 ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
+  list<TesselPointList *> *ListOfPaths = new list<TesselPointList *> ;
   TesselPointList *connectedPath = NULL;
   Vector center;
@@ -4033 +4007 @@
   class BoundaryLineSet *CurrentLine = NULL;
   class BoundaryLineSet *StartLine = NULL;
-
   // find the respective boundary point
   PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
@@ -4039 +4012 @@
     ReferencePoint = PointRunner->second;
   } else {
-    DoeLog(1) && (eLog()<< Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
     return NULL;
   }
 
-  map <class BoundaryLineSet *, bool> TouchedLine;
-  map <class BoundaryTriangleSet *, bool> TouchedTriangle;
-  map <class BoundaryLineSet *, bool>::iterator LineRunner;
-  map <class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
+  map<class BoundaryLineSet *, bool> TouchedLine;
+  map<class BoundaryTriangleSet *, bool> TouchedTriangle;
+  map<class BoundaryLineSet *, bool>::iterator LineRunner;
+  map<class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
   for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
-    TouchedLine.insert( pair <class BoundaryLineSet *, bool>(Runner->second, false) );
+    TouchedLine.insert(pair<class BoundaryLineSet *, bool> (Runner->second, false));
     for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++)
-      TouchedTriangle.insert( pair <class BoundaryTriangleSet *, bool>(Sprinter->second, false) );
+      TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool> (Sprinter->second, false));
   }
   if (!ReferencePoint->lines.empty()) {
@@ -4056 +4029 @@
       LineRunner = TouchedLine.find(runner->second);
       if (LineRunner == TouchedLine.end()) {
-        DoeLog(1) && (eLog()<< Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl);
+        DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl);
       } else if (!LineRunner->second) {
         LineRunner->second = true;
@@ -4064 +4037 @@
         StartLine = CurrentLine;
         CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
-        Log() << Verbose(1)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;
+        Log() << Verbose(1) << "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;
         do {
           // push current one
@@ -4086 +4059 @@
                 }
               } else {
-                DoeLog(1) && (eLog()<< Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl);
+                DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl);
                 triangle = NULL;
               }
@@ -4094 +4067 @@
             break;
           // find next line
-          for (int i=0;i<3;i++) {
+          for (int i = 0; i < 3; i++) {
             if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
               CurrentLine = triangle->lines[i];
@@ -4103 +4076 @@
           LineRunner = TouchedLine.find(CurrentLine);
           if (LineRunner == TouchedLine.end())
-            DoeLog(1) && (eLog()<< Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl);
+            DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl);
           else
             LineRunner->second = true;
@@ -4121 +4094 @@
     }
   } else {
-    DoeLog(1) && (eLog()<< Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl);
   }
 
@@ -4135 +4108 @@
 ListOfTesselPointList * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
-  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
+  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *> ;
   TesselPointList *connectedPath = NULL;
   TesselPointList *newPath = NULL;
   int count = 0;
-
-
   TesselPointList::iterator CircleRunner;
   TesselPointList::iterator CircleStart;
 
-  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  for (list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
@@ -4153 +4124 @@
     // go through list, look for reappearance of starting Point and count
     CircleStart = connectedPath->begin();
-
     // go through list, look for reappearance of starting Point and create list
     TesselPointList::iterator Marker = CircleStart;
@@ -4159 +4129 @@
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
-        Log() << Verbose(1) << count+1 << ". closed path consists of: ";
+        Log() << Verbose(1) << count + 1 << ". closed path consists of: ";
         newPath = new TesselPointList;
         TesselPointList::iterator CircleSprinter = Marker;
@@ -4181 +4151 @@
     connectedPath = *(ListofPaths->begin());
     ListofPaths->remove(connectedPath);
-    delete(connectedPath);
-  }
-  delete(ListofPaths);
+    delete (connectedPath);
+  }
+  delete (ListofPaths);
 
   // exit
   return ListofClosedPaths;
-};
-
+}
+;
 
 /** Gets all belonging triangles for a given BoundaryPointSet.
@@ -4197 +4167 @@
 TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
 {
-        Info FunctionInfo(__func__);
-        TriangleSet *connectedTriangles = new TriangleSet;
+  Info FunctionInfo(__func__);
+  TriangleSet *connectedTriangles = new TriangleSet;
 
   if (Point == NULL) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Point given is NULL." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "Point given is NULL." << endl);
   } else {
     // go through its lines and insert all triangles
     for (LineMap::const_iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
       for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
-      connectedTriangles->insert(TriangleRunner->second);
-    }
+        connectedTriangles->insert(TriangleRunner->second);
+      }
   }
 
   return connectedTriangles;
-};
-
+}
+;
 
 /** Removes a boundary point from the envelope while keeping it closed.
@@ -4225 +4195 @@
  * \return volume added to the volume inside the tesselated surface by the removal
  */
-double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point) {
+double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point)
+{
   class BoundaryLineSet *line = NULL;
   class BoundaryTriangleSet *triangle = NULL;
@@ -4233 +4204 @@
 
   if (point == NULL) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl);
     return 0.;
   } else
@@ -4243 +4214 @@
   // get list of connected points
   if (point->lines.empty()) {
-    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl);
     return 0.;
   }
@@ -4252 +4223 @@
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
-    count+=LineRunner->second->triangles.size();
+    count += LineRunner->second->triangles.size();
   TriangleMap Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
@@ -4258 +4229 @@
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
+      Candidates.insert(TrianglePair(triangle->Nr, triangle));
     }
   }
 
   // remove all triangles
-  count=0;
+  count = 0;
   NormalVector.Zero();
   for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
@@ -4280 +4251 @@
   double smallestangle;
   Vector Point, Reference, OrthogonalVector;
-  if (count > 2) {  // less than three triangles, then nothing will be created
+  if (count > 2) { // less than three triangles, then nothing will be created
     class TesselPoint *TriangleCandidates[3];
     count = 0;
-    for ( ; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) {  // go through all closed paths
+    for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
       if (ListAdvance != ListOfClosedPaths->end())
         ListAdvance++;
 
       connectedPath = *ListRunner;
-
       // re-create all triangles by going through connected points list
       LineList NewLines;
-      for (;!connectedPath->empty();) {
+      for (; !connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
         EndNode = connectedPath->end();
@@ -4317 +4287 @@
           angle = GetAngle(Point, Reference, OrthogonalVector);
           //if (angle < M_PI)  // no wrong-sided triangles, please?
-            if(fabs(angle - M_PI) < fabs(smallestangle - M_PI)) {  // get straightest angle (i.e. construct those triangles with smallest area first)
-              smallestangle = angle;
-              EndNode = MiddleNode;
-            }
+          if (fabs(angle - M_PI) < fabs(smallestangle - M_PI)) { // get straightest angle (i.e. construct those triangles with smallest area first)
+            smallestangle = angle;
+            EndNode = MiddleNode;
+          }
         }
         MiddleNode = EndNode;
         if (MiddleNode == connectedPath->end()) {
-          DoeLog(0) && (eLog()<< Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl);
+          DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl);
           performCriticalExit();
         }
@@ -4343 +4313 @@
         triangle = GetPresentTriangle(TriangleCandidates);
         if (triangle != NULL) {
-          DoeLog(0) && (eLog()<< Verbose(0) << "New triangle already present, skipping!" << endl);
+          DoeLog(0) && (eLog() << Verbose(0) << "New triangle already present, skipping!" << endl);
           StartNode++;
           MiddleNode++;
@@ -4355 +4325 @@
           continue;
         }
-        Log() << Verbose(3) << "Adding new triangle points."<< endl;
+        Log() << Verbose(3) << "Adding new triangle points." << endl;
         AddTesselationPoint(*StartNode, 0);
         AddTesselationPoint(*MiddleNode, 1);
         AddTesselationPoint(*EndNode, 2);
-        Log() << Verbose(3) << "Adding new triangle lines."<< endl;
+        Log() << Verbose(3) << "Adding new triangle lines." << endl;
         AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
         AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
@@ -4381 +4351 @@
           break;
         } else if (connectedPath->size() < 2) { // something's gone wrong!
-          DoeLog(0) && (eLog()<< Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl);
+          DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl);
           performCriticalExit();
         } else {
@@ -4401 +4371 @@
         do {
           maxgain = 0;
-          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for (LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -4419 +4389 @@
 
       ListOfClosedPaths->remove(connectedPath);
-      delete(connectedPath);
+      delete (connectedPath);
     }
     Log() << Verbose(0) << count << " triangles were created." << endl;
@@ -4427 +4397 @@
       connectedPath = *ListRunner;
       ListOfClosedPaths->remove(connectedPath);
-      delete(connectedPath);
+      delete (connectedPath);
     }
     Log() << Verbose(0) << "No need to create any triangles." << endl;
   }
-  delete(ListOfClosedPaths);
+  delete (ListOfClosedPaths);
 
   Log() << Verbose(0) << "Removed volume is " << volume << "." << endl;
 
   return volume;
-};
-
-
+}
+;
 
 /**
@@ -4450 +4419 @@
 TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
 {
-        Info FunctionInfo(__func__);
-        TriangleList *result = new TriangleList;
+  Info FunctionInfo(__func__);
+  TriangleList *result = new TriangleList;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
@@ -4475 +4444 @@
       for (int i = 0; i < 3; i++) {
         if (TrianglePoints[i] != NULL) {
-          for (int j = i+1; j < 3; j++) {
+          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++) {
+              (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);
@@ -4497 +4463 @@
     case 1: // copy all triangles of the respective line
     {
-      int i=0;
+      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++) {
+      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);
@@ -4516 +4479 @@
     case 2: // copy all triangles of the respective point
     {
-      int i=0;
+      int i = 0;
       for (; i < 3; i++)
         if (TrianglePoints[i] != NULL)
@@ -4534 +4497 @@
     }
     default:
-      DoeLog(0) && (eLog()<< Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl);
+      DoeLog(0) && (eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl);
       performCriticalExit();
       break;
@@ -4542 +4505 @@
 }
 
-struct BoundaryLineSetCompare {
-  bool operator() (const BoundaryLineSet * const a, const BoundaryLineSet * const b) {
+struct BoundaryLineSetCompare
+{
+  bool operator()(const BoundaryLineSet * const a, const BoundaryLineSet * const b)
+  {
     int lowerNra = -1;
     int lowerNrb = -1;
@@ -4561 +4526 @@
     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;
+    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;
-  };
+  }
+  ;
 };
 
@@ -4581 +4547 @@
 IndexToIndex * Tesselation::FindAllDegeneratedLines()
 {
-        Info FunctionInfo(__func__);
-        UniqueLines AllLines;
+  Info FunctionInfo(__func__);
+  UniqueLines AllLines;
   IndexToIndex * DegeneratedLines = new IndexToIndex;
 
   // sanity check
   if (LinesOnBoundary.empty()) {
-    DoeLog(2) && (eLog()<< Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.");
+    DoeLog(2) && (eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.");
     return DegeneratedLines;
   }
-
   LineMap::iterator LineRunner1;
-  pair< UniqueLines::iterator, bool> tester;
+  pair<UniqueLines::iterator, bool> tester;
   for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
-    tester = AllLines.insert( LineRunner1->second );
+    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) );
+      DegeneratedLines->insert(pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr));
+      DegeneratedLines->insert(pair<int, int> ((*tester.first)->Nr, LineRunner1->second->Nr));
     }
   }
@@ -4611 +4576 @@
       Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl;
     else
-      DoeLog(1) && (eLog()<< Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl);
+      DoeLog(1) && (eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl);
   }
 
@@ -4625 +4590 @@
 IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
   IndexToIndex * DegeneratedTriangles = new IndexToIndex;
-
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
   LineMap::iterator Liner;
@@ -4643 +4607 @@
     for (TriangleRunner1 = line1->triangles.begin(); TriangleRunner1 != line1->triangles.end(); ++TriangleRunner1) {
       for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) {
-        if ((TriangleRunner1->second != TriangleRunner2->second)
-          && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
-          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr) );
-          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr) );
+        if ((TriangleRunner1->second != TriangleRunner2->second) && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
+          DegeneratedTriangles->insert(pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr));
+          DegeneratedTriangles->insert(pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr));
         }
       }
     }
   }
-  delete(DegeneratedLines);
+  delete (DegeneratedLines);
 
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
   IndexToIndex::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
-      Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
+    Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
 
   return DegeneratedTriangles;
@@ -4667 +4630 @@
 void Tesselation::RemoveDegeneratedTriangles()
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
-  int count  = 0;
-
-  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
-    TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
-  ) {
+  int count = 0;
+
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner) {
     finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
     if (finder != TrianglesOnBoundary.end())
@@ -4692 +4653 @@
         trianglesShareLine = trianglesShareLine || triangle->lines[i] == partnerTriangle->lines[j];
 
-    if (trianglesShareLine
-      && (triangle->endpoints[1]->LinesCount > 2)
-      && (triangle->endpoints[2]->LinesCount > 2)
-      && (triangle->endpoints[0]->LinesCount > 2)
-    ) {
+    if (trianglesShareLine && (triangle->endpoints[1]->LinesCount > 2) && (triangle->endpoints[2]->LinesCount > 2) && (triangle->endpoints[0]->LinesCount > 2)) {
       // check whether we have to fix lines
       BoundaryTriangleSet *Othertriangle = NULL;
@@ -4716 +4673 @@
             // the line of triangle receives the degenerated ones
             triangle->lines[i]->triangles.erase(Othertriangle->Nr);
-            triangle->lines[i]->triangles.insert( TrianglePair( partnerTriangle->Nr, partnerTriangle) );
-            for (int k=0;k<3;k++)
+            triangle->lines[i]->triangles.insert(TrianglePair(partnerTriangle->Nr, partnerTriangle));
+            for (int k = 0; k < 3; k++)
               if (triangle->lines[i] == Othertriangle->lines[k]) {
                 Othertriangle->lines[k] = partnerTriangle->lines[j];
@@ -4723 +4680 @@
               }
             // the line of partnerTriangle receives the non-degenerated ones
-            partnerTriangle->lines[j]->triangles.erase( partnerTriangle->Nr);
-            partnerTriangle->lines[j]->triangles.insert( TrianglePair( Othertriangle->Nr, Othertriangle) );
+            partnerTriangle->lines[j]->triangles.erase(partnerTriangle->Nr);
+            partnerTriangle->lines[j]->triangles.insert(TrianglePair(Othertriangle->Nr, Othertriangle));
             partnerTriangle->lines[j] = triangle->lines[i];
           }
@@ -4736 +4693 @@
       RemoveTesselationTriangle(partnerTriangle);
     } else {
-      Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle
-        << " and its partner " << *partnerTriangle << " because it is essential for at"
-        << " least one of the endpoints to be kept in the tesselation structure." << endl;
-    }
-  }
-  delete(DegeneratedTriangles);
+      Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle << " and its partner " << *partnerTriangle << " because it is essential for at" << " least one of the endpoints to be kept in the tesselation structure." << endl;
+    }
+  }
+  delete (DegeneratedTriangles);
   if (count > 0)
     LastTriangle = NULL;
@@ -4758 +4713 @@
 void Tesselation::AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
@@ -4771 +4726 @@
     NearestBoundaryPoint = PointRunner->second;
   } else {
-    DoeLog(1) && (eLog()<< Verbose(1) << "I cannot find the boundary point." << endl);
+    DoeLog(1) && (eLog() << Verbose(1) << "I cannot find the boundary point." << endl);
     return;
   }
@@ -4789 +4744 @@
     CenterToPoint.SubtractVector(point->node);
     angle = CenterToPoint.Angle(&BaseLine);
-    if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
+    if (fabs(angle - M_PI / 2.) < fabs(BestAngle - M_PI / 2.)) {
       BestAngle = angle;
       BestLine = Runner->second;
@@ -4799 +4754 @@
   BestLine->triangles.erase(TempTriangle->Nr);
   int nr = -1;
-  for (int i=0;i<3; i++) {
+  for (int i = 0; i < 3; i++) {
     if (TempTriangle->lines[i] == BestLine) {
       nr = i;
@@ -4807 +4762 @@
 
   // create new triangle to connect point (connects automatically with the missing spot of the chosen line)
-  Log() << Verbose(2) << "Adding new triangle points."<< endl;
+  Log() << Verbose(2) << "Adding new triangle points." << endl;
   AddTesselationPoint((BestLine->endpoints[0]->node), 0);
   AddTesselationPoint((BestLine->endpoints[1]->node), 1);
   AddTesselationPoint(point, 2);
-  Log() << Verbose(2) << "Adding new triangle lines."<< endl;
+  Log() << Verbose(2) << "Adding new triangle lines." << endl;
   AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
   AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
@@ -4822 +4777 @@
 
   // create other side of this triangle and close both new sides of the first created triangle
-  Log() << Verbose(2) << "Adding new triangle points."<< endl;
+  Log() << Verbose(2) << "Adding new triangle points." << endl;
   AddTesselationPoint((BestLine->endpoints[0]->node), 0);
   AddTesselationPoint((BestLine->endpoints[1]->node), 1);
   AddTesselationPoint(point, 2);
-  Log() << Verbose(2) << "Adding new triangle lines."<< endl;
+  Log() << Verbose(2) << "Adding new triangle lines." << endl;
   AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
   AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
@@ -4836 +4791 @@
 
   // add removed triangle to the last open line of the second triangle
-  for (int i=0;i<3;i++) { // look for the same line as BestLine (only it's its degenerated companion)
+  for (int i = 0; i < 3; i++) { // look for the same line as BestLine (only it's its degenerated companion)
     if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
-      if (BestLine == BTS->lines[i]){
-        DoeLog(0) && (eLog()<< Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl);
+      if (BestLine == BTS->lines[i]) {
+        DoeLog(0) && (eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl);
         performCriticalExit();
       }
-      BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) );
+      BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle));
       TempTriangle->lines[nr] = BTS->lines[i];
       break;
     }
   }
-};
+}
+;
 
 /** Writes the envelope to file.
@@ -4856 +4812 @@
 void Tesselation::Output(const char *filename, const PointCloud * const cloud)
 {
-        Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   ofstream *tempstream = NULL;
   string NameofTempFile;
@@ -4862 +4818 @@
 
   if (LastTriangle != NULL) {
-    sprintf(NumberName, "-%04d-%s_%s_%s", (int)TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name);
+    sprintf(NumberName, "-%04d-%s_%s_%s", (int) TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name);
     if (DoTecplotOutput) {
       string NameofTempFile(filename);
       NameofTempFile.append(NumberName);
-      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-      NameofTempFile.erase(npos, 1);
+      for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
       NameofTempFile.append(TecplotSuffix);
       Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
@@ -4874 +4830 @@
       tempstream->close();
       tempstream->flush();
-      delete(tempstream);
+      delete (tempstream);
     }
 
@@ -4880 +4836 @@
       string NameofTempFile(filename);
       NameofTempFile.append(NumberName);
-      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-      NameofTempFile.erase(npos, 1);
+      for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
       NameofTempFile.append(Raster3DSuffix);
       Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
@@ -4889 +4845 @@
       tempstream->close();
       tempstream->flush();
-      delete(tempstream);
+      delete (tempstream);
     }
   }
   if (DoTecplotOutput || DoRaster3DOutput)
     TriangleFilesWritten++;
-};
-
-struct BoundaryPolygonSetCompare {
-  bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const {
+}
+;
+
+struct BoundaryPolygonSetCompare
+{
+  bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const
+  {
     if (s1->endpoints.size() < s2->endpoints.size())
       return true;
@@ -4926 +4885 @@
 {
   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;
+  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;
+    map<int, Vector *> TriangleVectors;
     // gather all NormalVectors
     Log() << Verbose(1) << "Gathering triangles ..." << endl;
@@ -4940 +4898 @@
       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)) );
+          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;
@@ -4949 +4907 @@
     // 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++)
+    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;
+          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));
@@ -4965 +4923 @@
         }
   }
-  delete(DegeneratedTriangles);
-
+  delete (DegeneratedTriangles);
   /// 3. Find connected endpoint candidates and put them into a polygon
   UniquePolygonSet ListofDegeneratedPolygons;
@@ -4972 +4929 @@
   BoundaryPointSet *OtherWalker = NULL;
   BoundaryPolygonSet *Current = NULL;
-  stack <BoundaryPointSet*> ToCheckConnecteds;
+  stack<BoundaryPointSet*> ToCheckConnecteds;
   while (!EndpointCandidateList.empty()) {
     Walker = *(EndpointCandidateList.begin());
-    if (Current == NULL) {  // create a new polygon with current candidate
+    if (Current == NULL) { // create a new polygon with current candidate
       Log() << Verbose(0) << "Starting new polygon set at point " << *Walker << endl;
       Current = new BoundaryPolygonSet;
@@ -4990 +4947 @@
         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
+        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
+          EndpointCandidateList.erase(Finder); // remove from candidates
+          ToCheckConnecteds.push(OtherWalker); // but check its partners too
         } else {
           Log() << Verbose(1) << " is not connected to " << *Walker << endl;
@@ -5015 +4972 @@
   /// 4. Go through all these degenerated polygons
   for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++) {
-    stack <int> TriangleNrs;
+    stack<int> TriangleNrs;
     Vector NormalVector;
     /// 4a. Gather all triangles of this polygon
@@ -5023 +4980 @@
     if (T->size() == 2) {
       Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl;
-      delete(T);
+      delete (T);
       continue;
     }
@@ -5032 +4989 @@
     // connections to either polygon ...
     if (T->size() % 2 != 0) {
-      DoeLog(0) && (eLog()<< Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl);
+      DoeLog(0) && (eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl);
       performCriticalExit();
     }
-
-    TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator
+    TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator
     /// 4a. Get NormalVector for one side (this is "front")
     NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
@@ -5059 +5015 @@
     /// 4c. Copy all "front" triangles but with inverse NormalVector
     TriangleWalker = T->begin();
-    while (TriangleWalker != T->end()) {  // go through all front triangles
+    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++)
@@ -5067 +5023 @@
       AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
       if (TriangleNrs.empty())
-        DoeLog(0) && (eLog()<< Verbose(0) << "No more free triangle numbers!" << endl);
+        DoeLog(0) && (eLog() << Verbose(0) << "No more free triangle numbers!" << endl);
       BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ...
       AddTesselationTriangle(); // ... and add
@@ -5076 +5032 @@
     }
     if (!TriangleNrs.empty()) {
-      DoeLog(0) && (eLog()<< Verbose(0) << "There have been less triangles created than removed!" << endl);
-    }
-    delete(T);  // remove the triangleset
-  }
-
+      DoeLog(0) && (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);
-
+    Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
+  delete (SimplyDegeneratedTriangles);
   /// 5. exit
   UniquePolygonSet::iterator PolygonRunner;
   while (!ListofDegeneratedPolygons.empty()) {
     PolygonRunner = ListofDegeneratedPolygons.begin();
-    delete(*PolygonRunner);
+    delete (*PolygonRunner);
     ListofDegeneratedPolygons.erase(PolygonRunner);
   }
 
   return counter;
-};
+}
+;

src/tesselation.hpp

@@ -89 +89 @@
 
 #define ListOfTesselPointList list<list <TesselPoint *> *>
+
+enum centers {Opt, OtherOpt};
 
 /********************************************** declarations *******************************/
@@ -284 +286 @@
     void AddTesselationTriangle();
     void AddTesselationTriangle(const int nr);
-    void AddCandidateTriangle(CandidateForTesselation &CandidateLine);
+    void AddCandidateTriangle(CandidateForTesselation &CandidateLine, enum centers type);
     void AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC);
     void AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC);
@@ -290 +292 @@
     void RemoveTesselationLine(class BoundaryLineSet *line);
     void RemoveTesselationPoint(class BoundaryPointSet *point);
-    bool CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const;
+    bool CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) const;
 
 
@@ -298 +300 @@
     void FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const;
     bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
+    bool FindCandidatesforOpenLines(const double RADIUS, const LinkedCell *&LCList);
     int CheckPresenceOfTriangle(class TesselPoint *Candidates[3]) const;
     class BoundaryTriangleSet * GetPresentTriangle(TesselPoint *Candidates[3]);
@@ -372 +375 @@
 
     //bool HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const OptCandidate, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const;
+    void FindDegeneratedCandidatesforOpenLines(TesselPoint * const Sprinter, const Vector * const OptCenter);
 };