Changes in src/Tesselation/tesselation.cpp [e2396e:8592c9]

File:

• src/Tesselation/tesselation.cpp (modified) (107 diffs)

src/Tesselation/tesselation.cpp

@@ -35 +35 @@
 #include "CodePatterns/MemDebug.hpp"
 
-#include <algorithm>
-#include <boost/foreach.hpp>
 #include <fstream>
 #include <iomanip>
-#include <iterator>
 #include <sstream>
 
@@ -69 +66 @@
 class molecule;
 
-const char *TecplotSuffix = ".dat";
-const char *Raster3DSuffix = ".r3d";
-const char *VRMLSUffix = ".wrl";
-
-const double ParallelEpsilon = 1e-3;
+const char *TecplotSuffix=".dat";
+const char *Raster3DSuffix=".r3d";
+const char *VRMLSUffix=".wrl";
+
+const double ParallelEpsilon=1e-3;
 const double Tesselation::HULLEPSILON = 1e-9;
 
@@ -79 +76 @@
  */
 Tesselation::Tesselation() :
-    PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin())
+  PointsOnBoundaryCount(0),
+  LinesOnBoundaryCount(0),
+  TrianglesOnBoundaryCount(0),
+  LastTriangle(NULL),
+  TriangleFilesWritten(0),
+  InternalPointer(PointsOnBoundary.begin())
 {
   //Info FunctionInfo(__func__);
@@ -111 +113 @@
 {
   // create linkedcell
-  LinkedCell_deprecated *LinkedList = new LinkedCell_deprecated(cloud, 2. * SPHERERADIUS);
-
-  // check for at least three points
-  {
-    bool ThreePointsFound = true;
-    cloud.GoToFirst();
-    for (size_t i = 0; i < 3; ++i, cloud.GoToNext())
+        LinkedCell_deprecated *LinkedList = new LinkedCell_deprecated(cloud, 2.*SPHERERADIUS);
+
+        // check for at least three points
+        {
+          bool ThreePointsFound = true;
+          cloud.GoToFirst();
+    for (size_t i=0;i<3;++i, cloud.GoToNext())
       ThreePointsFound &= (!cloud.IsEnd());
     cloud.GoToFirst();
@@ -124 +126 @@
       return;
     }
-  }
-
-  // find a starting triangle
+        }
+
+        // find a starting triangle
   FindStartingTriangle(SPHERERADIUS, LinkedList);
 
@@ -140 +142 @@
         //the line is there, so there is a triangle, but only one.
         const bool TesselationFailFlag = FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList);
-        ASSERT(TesselationFailFlag, "Tesselation::operator() - no suitable candidate triangle found.");
+        ASSERT( TesselationFailFlag,
+            "Tesselation::operator() - no suitable candidate triangle found.");
       }
     }
 
     // 2b. search for smallest ShortestAngle among all candidates
-    double ShortestAngle = 4. * M_PI;
+    double ShortestAngle = 4.*M_PI;
     for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) {
       if (Runner->second->ShortestAngle < ShortestAngle) {
@@ -152 +155 @@
       }
     }
-    if ((ShortestAngle == 4. * M_PI) || (baseline->pointlist.empty()))
+    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
       OneLoopWithoutSuccessFlag = false;
     else {
@@ -158 +161 @@
     }
   }
-
-  delete LinkedList;
 }
 
@@ -169 +170 @@
 double Tesselation::getVolumeOfConvexEnvelope(const bool IsAngstroem) const
 {
-  // calculate center of gravity
-  Vector center;
-  if (!PointsOnBoundary.empty()) {
-    for (PointMap::const_iterator iter = PointsOnBoundary.begin();
-        iter != PointsOnBoundary.end(); ++iter)
-      center += iter->second->node->getPosition();
-    center *= 1./(double)PointsOnBoundary.size();
-  }
+  double volume = 0.;
+  Vector x;
+  Vector y;
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
-  double volume = 0.;
-  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
-    const double TetrahedronVolume = CalculateVolumeofGeneralTetraeder(
-        runner->second->endpoints[0]->getPosition(),
-        runner->second->endpoints[1]->getPosition(),
-        runner->second->endpoints[2]->getPosition(),
-        center);
-    LOG(1, "INFO: volume of tetrahedron is " << setprecision(10) << TetrahedronVolume
-        << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
-    volume += TetrahedronVolume;
-  }
-  LOG(0, "RESULT: The summed volume is " << setprecision(6) << volume
-      << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++)
+    { // go through every triangle, calculate volume of its pyramid with CoG as peak
+      x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1);
+      const double G = runner->second->getArea();
+      x = runner->second->getPlane().getNormal();
+      x.Scale(runner->second->getEndpoint(1).ScalarProduct(x));
+      const double h = x.Norm(); // distance of CoG to triangle
+      const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
+      LOG(1, "INFO: Area of triangle is " << setprecision(10) << G << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
+          << h << " and the volume is " << PyramidVolume << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+      volume += PyramidVolume;
+    }
+  LOG(0, "RESULT: The summed volume is " << setprecision(6)
+      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
 
   return volume;
@@ -208 +207 @@
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
-  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
-    const double area = runner->second->getArea();
-    LOG(1, "INFO: Area of triangle is " << setprecision(10) << area << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2.");
-    surfacearea += area;
-  }
-  LOG(0, "RESULT: The summed surface area is " << setprecision(6) << surfacearea << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+  for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++)
+    { // go through every triangle, calculate volume of its pyramid with CoG as peak
+                const double area = runner->second->getArea();
+                LOG(1, "INFO: Area of triangle is " << setprecision(10) << area << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2.");
+      surfacearea += area;
+    }
+  LOG(0, "RESULT: The summed surface area is " << setprecision(6)
+      << surfacearea << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
 
   return surfacearea;
 }
+
 
 /** Gueses first starting triangle of the convex envelope.
@@ -248 +251 @@
         tmp = B->second->node->DistanceSquared(C->second->node->getPosition());
         distance += tmp * tmp;
-        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(B, C)));
+        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
       }
     }
@@ -268 +271 @@
     // 2. next, we have to check whether all points reside on only one side of the triangle
     // 3. construct plane vector
-    PlaneVector = Plane(A->second->node->getPosition(), baseline->second.first->second->node->getPosition(), baseline->second.second->second->node->getPosition()).getNormal();
+    PlaneVector = Plane(A->second->node->getPosition(),
+                        baseline->second.first->second->node->getPosition(),
+                        baseline->second.second->second->node->getPosition()).getNormal();
     LOG(2, "Plane vector of candidate triangle is " << PlaneVector);
     // 4. loop over all points
@@ -383 +388 @@
         // prepare some auxiliary vectors
         Vector BaseLineCenter, BaseLine;
-        BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) + (baseline->second->endpoints[1]->node->getPosition()));
+        BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                (baseline->second->endpoints[1]->node->getPosition()));
         BaseLine = (baseline->second->endpoints[0]->node->getPosition()) - (baseline->second->endpoints[1]->node->getPosition());
 
@@ -401 +407 @@
         // vector in propagation direction (out of triangle)
         // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
-        PropagationVector = Plane(BaseLine, NormalVector, 0).getNormal();
+        PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
         TempVector = CenterVector - (baseline->second->endpoints[0]->node->getPosition()); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
         //LOG(0, "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << ".");
@@ -454 +460 @@
             // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
             flag = true;
-            VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()), (baseline->second->endpoints[1]->node->getPosition()), (target->second->node->getPosition())).getNormal();
-            TempVector = (1. / 3.) * ((baseline->second->endpoints[0]->node->getPosition()) + (baseline->second->endpoints[1]->node->getPosition()) + (target->second->node->getPosition()));
+            VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()),
+                                        (baseline->second->endpoints[1]->node->getPosition()),
+                                        (target->second->node->getPosition())).getNormal();
+            TempVector = (1./3.) * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                    (baseline->second->endpoints[1]->node->getPosition()) +
+                                    (target->second->node->getPosition()));
             TempVector -= (*Center);
             // make it always point outward
@@ -467 +477 @@
               winner = target;
               LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle between normal vectors.");
+            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
+              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
+              helper = (target->second->node->getPosition()) - BaseLineCenter;
+              helper.ProjectOntoPlane(BaseLine);
+              // ...the one with the smaller angle is the better candidate
+              TempVector = (target->second->node->getPosition()) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              TempAngle = TempVector.Angle(helper);
+              TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              if (TempAngle < TempVector.Angle(helper)) {
+                TempAngle = NormalVector.Angle(VirtualNormalVector);
+                SmallestAngle = TempAngle;
+                winner = target;
+                LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction.");
+              } else
+                LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction.");
             } else
-              if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
-                // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
-                helper = (target->second->node->getPosition()) - BaseLineCenter;
-                helper.ProjectOntoPlane(BaseLine);
-                // ...the one with the smaller angle is the better candidate
-                TempVector = (target->second->node->getPosition()) - BaseLineCenter;
-                TempVector.ProjectOntoPlane(VirtualNormalVector);
-                TempAngle = TempVector.Angle(helper);
-                TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
-                TempVector.ProjectOntoPlane(VirtualNormalVector);
-                if (TempAngle < TempVector.Angle(helper)) {
-                  TempAngle = NormalVector.Angle(VirtualNormalVector);
-                  SmallestAngle = TempAngle;
-                  winner = target;
-                  LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction.");
-                } else
-                  LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction.");
-              } else
-                LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors.");
+              LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors.");
           }
         } // end of loop over all boundary points
@@ -553 +562 @@
 
   cloud.GoToFirst();
-  PointCloudAdaptor<Tesselation, MapValueIterator<Tesselation::iterator> > newcloud(this, cloud.GetName());
+  PointCloudAdaptor< Tesselation, MapValueIterator<Tesselation::iterator> > newcloud(this, cloud.GetName());
   BoundaryPoints = new LinkedCell_deprecated(newcloud, 5.);
   while (!cloud.IsEnd()) { // we only have to go once through all points, as boundary can become only bigger
@@ -720 +729 @@
  * @param *b second endpoint
  * @param n index of Tesselation::BLS giving the line with both endpoints
- * @return true - inserted a new line, false - present line used
- */
-bool 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;
@@ -738 +746 @@
       if (FindLine->second->triangles.size() == 1) {
         CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
-        ASSERT(Finder != OpenLines.end(), "Tesselation::AddTesselationLine() - " + toString(*FindLine->second) + " is not a new line and not in OpenLines.");
-        if (Finder->second != NULL) {
-          if (!Finder->second->pointlist.empty())
-            LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << ".");
-          else {
-            LOG(4, "ACCEPT: line " << *(FindLine->second) << " is open with no candidate.");
+        if (!Finder->second->pointlist.empty())
+          LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << ".");
+        else
+          LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with no candidate.");
+        // get open line
+        for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
+          if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches
+            LOG(4, "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << ".");
             insertNewLine = false;
             WinningLine = FindLine->second;
+            break;
+          } else {
+            LOG(5, "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << ".");
           }
-          // get open line
-          for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
-            if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON)) { // stop searching if candidate matches
-              LOG(4, "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << ".");
-              insertNewLine = false;
-              WinningLine = FindLine->second;
-              break;
-            } else {
-              LOG(5, "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << ".");
-            }
-          }
-        } else {
-          LOG(4, "ACCEPT: There are no candidates for present open line, use what we have.");
-          insertNewLine = false;
-          WinningLine = FindLine->second;
         }
       }
@@ -772 +770 @@
     AddExistingTesselationTriangleLine(WinningLine, n);
   }
-
-  return insertNewLine;
 }
 ;
@@ -798 +794 @@
   // 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));
 }
 ;
@@ -887 +883 @@
       } else {
         output << "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++)
           output << "\t[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t";
-        OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *>(triangle->lines[i], NULL));
       }
       LOG(3, "DEBUG: " << output.str());
@@ -923 +919 @@
   else
     Numbers[1] = -1;
-
-  // erase from OpenLines if present
-  {
-    CandidateMap::iterator openlineiter = OpenLines.find(line);
-    if (openlineiter != OpenLines.end())
-      OpenLines.erase(openlineiter);
-  }
 
   for (int i = 0; i < 2; i++) {
@@ -948 +937 @@
         LOG(4, "DEBUG: " << *line->endpoints[i] << " has no more lines it's attached to, erasing.");
         RemoveTesselationPoint(line->endpoints[i]);
-      } else
-        if (DoLog(0)) {
-          std::stringstream output;
-          output << "DEBUG: " << *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++)
-            output << "[" << *(LineRunner->second) << "] \t";
-          LOG(4, output.str());
-        }
+      } else if (DoLog(0)) {
+        std::stringstream output;
+        output << "DEBUG: " << *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++)
+          output << "[" << *(LineRunner->second) << "] \t";
+        LOG(4, output.str());
+      }
       line->endpoints[i] = NULL; // free'd or not: disconnect
     } else
@@ -997 +985 @@
   //Info FunctionInfo(__func__);
 
-  LOG(3, "DEBUG: Checking degeneracy by whether sphere contains no others points ...");
+  LOG(3, "DEBUG: Checking whether sphere contains no others points ...");
   bool flag = true;
 
@@ -1004 +992 @@
   TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);
 
-  LOG(3, "DEBUG: CheckDegeneracy: There are " << ListofPoints->size() << " points inside the sphere.");
-  LOG(4, "DEBUG: The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":");
+  LOG(3, "DEBUG: The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":");
   for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-    LOG(4, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
+    LOG(3, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
 
   // remove triangles's endpoints
@@ -1021 +1008 @@
     LOG(3, "DEBUG: CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere.");
     flag = false;
-    LOG(4, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":");
+    LOG(3, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":");
     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      LOG(4, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
-  }
-  delete ListofPoints;
+      LOG(3, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
+  }
+  delete (ListofPoints);
 
   return flag;
@@ -1168 +1155 @@
 
   // 1. searching topmost point with respect to each axis
-  LOG(2, "INFO: Searching for topmost pointer from each axis.");
   for (int i = 0; i < NDIM; i++) { // each axis
     LC->n[i] = LC->N[i] - 1; // current axis is topmost cell
@@ -1190 +1176 @@
   }
 
-  if (DoLog(3)) {
+  if (DoLog(1)) {
     std::stringstream output;
     output << "Found maximum coordinates: ";
@@ -1204 +1190 @@
     BaseLine = new BoundaryLineSet();
     BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
-    LOG(2, "INFO: Coordinates of start node at " << *BaseLine->endpoints[0]->node << ".");
+    LOG(2, "DEBUG: Coordinates of start node at " << *BaseLine->endpoints[0]->node << ".");
 
     double ShortestAngle;
@@ -1217 +1203 @@
     }
     BaseLine->endpoints[1] = new BoundaryPointSet(Temporary);
-    LOG(2, "INFO: Second node is at " << *Temporary << ".");
+    LOG(1, "INFO: Second node is at " << *Temporary << ".");
 
     // construct center of circle
     CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
-    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ".");
+    LOG(1, "INFO: CircleCenter is at " << CircleCenter << ".");
 
     // construct normal vector of circle
     CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
-    LOG(4, "DEBUG: CirclePlaneNormal is at " << CirclePlaneNormal << ".");
+    LOG(1, "INFO: CirclePlaneNormal is at " << CirclePlaneNormal << ".");
 
     double radius = CirclePlaneNormal.NormSquared();
@@ -1232 +1218 @@
     NormalVector.ProjectOntoPlane(CirclePlaneNormal);
     NormalVector.Normalize();
-    LOG(4, "DEBUG: NormalVector is at " << NormalVector << ".");
+    LOG(1, "INFO: NormalVector is at " << NormalVector << ".");
     ShortestAngle = 2. * M_PI; // This will indicate the quadrant.
 
-    SphereCenter = (CircleRadius * NormalVector) + CircleCenter;
+    SphereCenter = (CircleRadius * NormalVector) +  CircleCenter;
     // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
 
     // look in one direction of baseline for initial candidate
     try {
-      SearchDirection = Plane(CirclePlaneNormal, NormalVector, 0).getNormal(); // whether we look "left" first or "right" first is not important ...
-    } catch (LinearAlgebraException &e) {
-      ELOG(1, "Vectors are linear dependent: " << CirclePlaneNormal << ", " << NormalVector << ".");
+      SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal();  // whether we look "left" first or "right" first is not important ...
+    } catch(LinearAlgebraException) {
+      ELOG(1, "Vectors are linear dependent: "
+          << CirclePlaneNormal << ", " << NormalVector << ".");
       delete BaseLine;
       continue;
@@ -1248 +1235 @@
 
     // adding point 1 and point 2 and add the line between them
-    LOG(3, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << ".");
+    LOG(2, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << ".");
 
     //LOG(1, "INFO: OldSphereCenter is at " << helper << ".");
@@ -1257 +1244 @@
       for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++)
         output << *(*it);
-      LOG(3, "DEBUG: List of third Points is: " << output.str());
+      LOG(2, "DEBUG: List of third Points is: " << output.str());
     }
     if (!OptCandidates.pointlist.empty()) {
@@ -1419 +1406 @@
 //  return result;
 //};
+
 /** This function finds a triangle to a line, adjacent to an existing one.
  * @param out output stream for debugging
@@ -1450 +1438 @@
 
   // construct center of circle
-  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
 
   // construct normal vector of circle
-  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
 
   // calculate squared radius of circle
@@ -1463 +1453 @@
     CircleRadius = RADIUS * RADIUS - radius / 4.;
     CirclePlaneNormal.Normalize();
-    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
-
-    LOG(4, "DEBUG: OldSphereCenter is at " << T.SphereCenter << ".");
+    LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+
+    LOG(3, "DEBUG: OldSphereCenter is at " << T.SphereCenter << ".");
 
     // construct SearchDirection and an "outward pointer"
-    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal, 0).getNormal();
+    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
     helper = CircleCenter - (ThirdPoint->node->getPosition());
-    if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards!
+    if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
-    LOG(4, "DEBUG: SearchDirection is " << SearchDirection << ".");
+    LOG(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
     if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
@@ -1482 +1472 @@
 
   } else {
-    LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!");
+    LOG(3, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!");
   }
 
@@ -1515 +1505 @@
     baseline = Runner->second;
     if (baseline->pointlist.empty()) {
-      ASSERT((baseline->BaseLine->triangles.size() == 1), "Open line without exactly one attached triangle");
+      ASSERT((baseline->BaseLine->triangles.size() == 1),"Open line without exactly one attached triangle");
       T = (((baseline->BaseLine->triangles.begin()))->second);
       LOG(4, "DEBUG: Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T);
@@ -1548 +1538 @@
   TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->getPosition());
 
-  if (DoLog(3)) {
+  {
     std::stringstream output;
     for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
       output << **TesselRunner;
-    LOG(3, "DEBUG: List of Candidates for Turning Point " << *TurningPoint << ":" << output.str());
+    LOG(3, "DEBUG: List of Candidates for Turning Point " << *TurningPoint << ":");
   }
 
@@ -1593 +1583 @@
   }
   delete (connectedClosestPoints);
-}
-;
+};
 
 /** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate.
@@ -1614 +1603 @@
       else {
         LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter));
-        Finder->second->T = BTS; // is last triangle
+        Finder->second->T = BTS;  // is last triangle
         Finder->second->pointlist.push_back(Sprinter);
         Finder->second->ShortestAngle = 0.;
@@ -1621 +1610 @@
     }
   }
-}
-;
+};
 
 /** If a given \a *triangle is degenerated, this adds both sides.
@@ -1658 +1646 @@
   // give some verbose output about the whole procedure
   if (CandidateLine.T != NULL)
-    LOG(2, "INFO: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << ".");
+    LOG(2, "DEBUG: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << ".");
   else
-    LOG(2, "INFO: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle.");
+    LOG(2, "DEBUG: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle.");
   triangle = BTS;
 
@@ -1744 +1732 @@
 ;
 
-bool Tesselation::isConvex() const
-{
-  bool status = true;
-  // go through all lines on boundary
-  for (LineMap::const_iterator lineiter = LinesOnBoundary.begin();
-      lineiter != LinesOnBoundary.end(); ++lineiter) {
-    if (!lineiter->second->CheckConvexityCriterion()) {
-      LOG(2, "DEBUG: Line " << *lineiter->second << " is not convex.");
-      status = false;
-    }
-  }
-  return status;
-}
-
 /** Checks whether the quadragon of the two triangles connect to \a *Base is convex.
  * We look whether the closest point on \a *Base with respect to the other baseline is outside
@@ -1894 +1868 @@
     }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-      LOG(4, "DEBUG: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
+    LOG(4, "DEBUG: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
       BaseLineNormal += (runner->second->NormalVector);
     }
@@ -1923 +1897 @@
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
-  Vector BaseLineNormal[2];
-  Vector OtherEndpoint[2]; // fourth point to either triangle
+  Vector BaseLineNormal;
   int OldTriangleNrs[2], OldBaseLineNr;
   int i, m;
 
   // calculate NormalVector for later use
+  BaseLineNormal.Zero();
   if (Base->triangles.size() < 2) {
     ELOG(1, "Less than two triangles are attached to this baseline!");
     return NULL;
   }
-  {
-    int i=0;
-    for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
-      LOG(5, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
-      OtherEndpoint[(i+1)%2] = runner->second->GetThirdEndpoint(Base)->node->getPosition();
-      BaseLineNormal[i++] = (runner->second->NormalVector);
-      ASSERT( i <= 2,
-          "Tesselation::FlipBaseline() - not exactly two triangles found.");
-    }
-  }
+  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
+    LOG(1, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");
+    BaseLineNormal += (runner->second->NormalVector);
+  }
+  BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
 
   // get the two triangles
@@ -1972 +1941 @@
 
   // index OldLines and OldPoints
-  // note that oldlines[0,1] belong to first triangle and hence first normal
-  // vector and oldlines[2,3] belong to second triangle and its normal vector
   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
@@ -1983 +1950 @@
         OldPoints[m++] = runner->second->endpoints[j];
 
-  Vector BasePoints[2]; // endpoints of Base
-  if (OldLines[0]->ContainsBoundaryPoint(Base->endpoints[0])) {
-    BasePoints[0] = Base->endpoints[0]->node->getPosition();
-    BasePoints[1] = Base->endpoints[1]->node->getPosition();
-  } else {
-    BasePoints[0] = Base->endpoints[1]->node->getPosition();
-    BasePoints[1] = Base->endpoints[0]->node->getPosition();
-  }
   // check whether everything is in place to create new lines and triangles
   if (i < 4) {
@@ -2006 +1965 @@
       return NULL;
     }
-
-  // construct plane of first triangle for calculating normal vectors later
-  const Plane triangleplane = Base->triangles.begin()->second->getPlane();
-  // get fourth point projected down onto this plane
-  const Vector Intersection = triangleplane.getClosestPoint(OtherEndpoint[0]);
 
   // remove triangles and baseline removes itself
@@ -2017 +1971 @@
   m = 0;
   // first obtain all triangle to delete ... (otherwise we pull the carpet (Base) from under the for-loop's feet)
-  list<BoundaryTriangleSet *> TrianglesOfBase;
+  list <BoundaryTriangleSet *> TrianglesOfBase;
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); ++runner)
     TrianglesOfBase.push_back(runner->second);
   // .. then delete each triangle (which deletes the line as well)
-  for (list<BoundaryTriangleSet *>::iterator runner = TrianglesOfBase.begin(); !TrianglesOfBase.empty(); runner = TrianglesOfBase.begin()) {
+  for (list <BoundaryTriangleSet *>::iterator runner = TrianglesOfBase.begin(); !TrianglesOfBase.empty(); runner = TrianglesOfBase.begin()) {
     LOG(3, "DEBUG: Deleting triangle " << *(*runner) << ".");
     OldTriangleNrs[m++] = (*runner)->Nr;
@@ -2034 +1988 @@
   LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
   LOG(3, "DEBUG: Created new baseline " << *NewLine << ".");
-
-  // Explanation for normal vector choice:
-  // Decisive for the normal vector of the new triangle is whether the fourth
-  // endpoint is with respect to the joining boundary line on one side or on
-  // the other with respect to the endpoint of the plane triangle that is not
-  // contained in the joining boundary line.
 
   // construct new triangles with flipped baseline
@@ -2052 +2000 @@
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
-    {
-      Line joiningline = makeLineThrough(
-          OldLines[0]->endpoints[0]->node->getPosition(),
-          OldLines[0]->endpoints[1]->node->getPosition());
-      // BasePoints[1] is not contained in OldLines[0], hence is third endpoint
-      // of plane triangle. BasePoints[0] is in the joining OldLines[0] and
-      // we check whether Intersection is on same side as BasePoints[1] or not.
-      const Vector pointinginside =
-          joiningline.getClosestPoint(BasePoints[1]) - BasePoints[1];
-      const Vector pointingtointersection =
-          joiningline.getClosestPoint(Intersection) - Intersection;
-      const double sign_of_intersection =
-          pointingtointersection.ScalarProduct(pointinginside);
-      LOG(4, "DEBUG: Sign of SKP between both lines w.r.t " << *OldLines[0]
-          << " is " << sign_of_intersection << ".");
-      const double sign_of_normal = (sign_of_intersection >= 0) ? -1. : +1.;
-      LOG(4, "DEBUG: Opposite normal direction is "
-          << sign_of_normal * BaseLineNormal[0] << ".");
-      BTS->GetNormalVector(sign_of_normal * BaseLineNormal[0]);
-    }
+    BTS->GetNormalVector(BaseLineNormal);
     AddTesselationTriangle(OldTriangleNrs[0]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");
@@ -2079 +2008 @@
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
-    {
-      Line joiningline = makeLineThrough(
-          OldLines[1]->endpoints[0]->node->getPosition(),
-          OldLines[1]->endpoints[1]->node->getPosition());
-      // BasePoints[0] is not contained in OldLines[1], hence is third endpoint
-      // of plane triangle. BasePoints[1] is in th1e joining OldLines[1] and
-      // we check whether Intersection is on same side as BasePoints[0] or not.
-      const Vector pointinginside =
-          joiningline.getClosestPoint(BasePoints[0]) - BasePoints[0];
-      const Vector pointingtointersection =
-          joiningline.getClosestPoint(Intersection) - Intersection;
-      const double sign_of_intersection =
-          pointingtointersection.ScalarProduct(pointinginside);
-      LOG(4, "DEBUG: Sign of SKP between both lines w.r.t " << *OldLines[1]
-          << " is " << sign_of_intersection << ".");
-      const double sign_of_normal = (sign_of_intersection >= 0) ? -1. : +1.;
-      LOG(4, "DEBUG: Opposite normal direction is "
-          << sign_of_normal * BaseLineNormal[0] << ".");
-      BTS->GetNormalVector(sign_of_normal * BaseLineNormal[0]);
-    }
+    BTS->GetNormalVector(BaseLineNormal);
     AddTesselationTriangle(OldTriangleNrs[1]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");
@@ -2245 +2155 @@
   TesselPoint *Candidate = NULL;
 
-  LOG(4, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << ".");
+  LOG(3, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << ".");
 
   // copy old center
@@ -2253 +2163 @@
 
   // construct center of circle
-  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
 
   // construct normal vector of circle
-  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
 
   RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
@@ -2265 +2177 @@
     CircleRadius = RADIUS * RADIUS - radius;
     CirclePlaneNormal.Normalize();
-    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+    LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
 
     // test whether old center is on the band's plane
@@ -2274 +2186 @@
     radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      LOG(4, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << ".");
+      LOG(3, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << ".");
 
       // check SearchDirection
-      LOG(4, "DEBUG: SearchDirection is " << SearchDirection << ".");
+      LOG(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
       if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
         ELOG(1, "SearchDirection and RelativeOldSphereCenter are not orthogonal!");
@@ -2318 +2230 @@
                   // find center on the plane
                   GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
-                  LOG(5, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << ".");
+                  LOG(3, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << ".");
 
                   try {
-                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()), (CandidateLine.BaseLine->endpoints[1]->node->getPosition()), (Candidate->getPosition())).getNormal();
-                    LOG(5, "DEBUG: NewNormalVector is " << NewNormalVector << ".");
+                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()),
+                                            (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
+                                            (Candidate->getPosition())).getNormal();
+                    LOG(3, "DEBUG: NewNormalVector is " << NewNormalVector << ".");
                     radius = CandidateLine.BaseLine->endpoints[0]->node->DistanceSquared(NewPlaneCenter);
-                    LOG(5, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
-                    LOG(5, "DEBUG: SearchDirection is " << SearchDirection << ".");
-                    LOG(5, "DEBUG: Radius of CircumCenterCircle is " << radius << ".");
+                    LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+                    LOG(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
+                    LOG(3, "DEBUG: Radius of CircumCenterCircle is " << radius << ".");
                     if (radius < RADIUS * RADIUS) {
                       otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter);
@@ -2332 +2246 @@
                         // construct both new centers
                         NewSphereCenter = NewPlaneCenter;
-                        OtherNewSphereCenter = NewPlaneCenter;
+                        OtherNewSphereCenter= NewPlaneCenter;
                         helper = NewNormalVector;
                         helper.Scale(sqrt(RADIUS * RADIUS - radius));
-                        LOG(5, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << ".");
+                        LOG(4, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << ".");
                         NewSphereCenter += helper;
-                        LOG(5, "DEBUG: NewSphereCenter is at " << NewSphereCenter << ".");
+                        LOG(4, "DEBUG: NewSphereCenter is at " << NewSphereCenter << ".");
                         // OtherNewSphereCenter is created by the same vector just in the other direction
                         helper.Scale(-1.);
                         OtherNewSphereCenter += helper;
-                        LOG(5, "DEBUG: OtherNewSphereCenter is at " << OtherNewSphereCenter << ".");
+                        LOG(4, "DEBUG: OtherNewSphereCenter is at " << OtherNewSphereCenter << ".");
                         alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
                         Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
@@ -2362 +2276 @@
                           if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
                             CandidateLine.pointlist.push_back(Candidate);
-                            LOG(3, "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
+                            LOG(2, "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
                           } else {
                             // remove all candidates from the list and then the list itself
                             CandidateLine.pointlist.clear();
                             CandidateLine.pointlist.push_back(Candidate);
-                            LOG(3, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
+                            LOG(2, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
                           }
                           CandidateLine.ShortestAngle = alpha;
-                          LOG(4, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now.");
+                          LOG(2, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now.");
                         } else {
                           if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
-                            LOG(4, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " .");
+                            LOG(3, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " .");
                           } else {
-                            LOG(4, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected.");
+                            LOG(3, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected.");
                           }
                         }
                       } else {
-                        ASSERT(0, std::string("FindThirdPointForTesselation() - ") + std::string("REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: ") + toString(fabs(radius - otherradius)));
+                        ELOG(0, "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius));
                       }
                     } else {
-                      LOG(4, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << ".");
+                      LOG(3, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << ".");
                     }
-                  } catch (LinearDependenceException &excp) {
-                    LOG(4, boost::diagnostic_information(excp));
-                    LOG(4, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent.");
+                  }
+                  catch (LinearDependenceException &excp){
+                    LOG(3, boost::diagnostic_information(excp));
+                    LOG(3, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent.");
                   }
                 } else {
                   if (ThirdPoint != NULL) {
-                    LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << ".");
+                    LOG(3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << ".");
                   } else {
-                    LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << ".");
+                    LOG(3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << ".");
                   }
                 }
@@ -2403 +2318 @@
   } else {
     if (ThirdPoint != NULL)
-      LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!");
+      LOG(3, "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!");
     else
-      LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!");
-  }
-
-  LOG(5, "DEBUG: Sorting candidate list ...");
+      LOG(3, "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!");
+  }
+
+  LOG(2, "DEBUG: Sorting candidate list ...");
   if (CandidateLine.pointlist.size() > 1) {
     CandidateLine.pointlist.unique();
@@ -2414 +2329 @@
   }
 
-  ASSERT(CandidateLine.pointlist.empty() || (CandidateLine.CheckValidity(RADIUS, LC)), std::string("Tesselation::FindThirdPointForTesselation()") + std::string("There were other points contained in the rolling sphere as well!"));
+  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
+    ELOG(0, "There were other points contained in the rolling sphere as well!");
+    performCriticalExit();
+  }
 }
 ;
@@ -2426 +2344 @@
 {
   //Info FunctionInfo(__func__);
-  const BoundaryLineSet * lines[2] = {line1, line2};
+  const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
   PointMap OrderMap;
@@ -2433 +2351 @@
     // 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]));
+      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;
@@ -2482 +2400 @@
               // we should make sure that both triangles end up in the same entry
               // in the distance multimap. Hence, we round to 6 digit precision.
-              const double distance = 1e-6 * floor(FindPoint->second->node->DistanceSquared(x) * 1e+6);
+              const double distance =
+                  1e-6*floor(FindPoint->second->node->DistanceSquared(x)*1e+6);
               points->insert(DistanceToPointPair(distance, FindPoint->second));
               LOG(3, "DEBUG: Putting " << *FindPoint->second << " into the list.");
@@ -2527 +2446 @@
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
       // calculate closest point on line to desired point
-      helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) + ((LineRunner->second)->endpoints[1]->node->getPosition()));
+      helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) +
+                      ((LineRunner->second)->endpoints[1]->node->getPosition()));
       Center = (x) - helper;
-      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - ((LineRunner->second)->endpoints[1]->node->getPosition());
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
       Center.ProjectOntoPlane(BaseLine);
       const double distance = Center.NormSquared();
@@ -2586 +2507 @@
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
 
-      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - ((LineRunner->second)->endpoints[1]->node->getPosition());
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
       const double lengthBase = BaseLine.NormSquared();
 
@@ -2602 +2524 @@
           MinDistance = lengthEnd;
           LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << ".");
-        } else
-          if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
-            ClosestLines.insert(LineRunner->second);
-            LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << ".");
-          } else { // line is worse
-            LOG(1, "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << ".");
-          }
+        } else if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
+          ClosestLines.insert(LineRunner->second);
+          LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << ".");
+        } else { // line is worse
+          LOG(1, "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << ".");
+        }
       } else { // intersection is closer, calculate
         // calculate closest point on line to desired point
@@ -2734 +2655 @@
   double distance = 0.;
 
-  if (triangle == NULL) { // is boundary point or only point in point cloud?
+  if (triangle == NULL) {// is boundary point or only point in point cloud?
     LOG(1, "No triangle given!");
     return -1.;
@@ -2843 +2764 @@
       takePoint = true;
       current = findLines->second->endpoints[1]->node;
-    } else
-      if (findLines->second->endpoints[1]->Nr == Point->getNr()) {
-        takePoint = true;
-        current = findLines->second->endpoints[0]->node;
-      }
+    } else if (findLines->second->endpoints[1]->Nr == Point->getNr()) {
+      takePoint = true;
+      current = findLines->second->endpoints[0]->node;
+    }
 
     if (takePoint) {
@@ -2887 +2807 @@
   Vector OrthogonalVector;
   Vector helper;
-  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
+  const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL };
   TriangleList *triangles = NULL;
 
@@ -2898 +2818 @@
   // calculate central point
   triangles = FindTriangles(TrianglePoints);
-  ASSERT((triangles == NULL) || (triangles->empty()), std::string("Tesselation::GetCircleOfConnectedTriangles()") + std::string("Could not find any triangles for point " + toString(*Point) + "."));
-  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-    PlaneNormal += (*Runner)->NormalVector;
+  if ((triangles != NULL) && (!triangles->empty())) {
+    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      PlaneNormal += (*Runner)->NormalVector;
+  } else {
+    ELOG(0, "Could not find any triangles for point " << *Point << ".");
+    performCriticalExit();
+  }
   PlaneNormal.Scale(1.0 / triangles->size());
   LOG(4, "DEBUG: Calculated PlaneNormal of all circle points is " << PlaneNormal << ".");
@@ -2912 +2836 @@
     AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
-    ASSERT(AngleZero.NormSquared() > MYEPSILON, std::string("Tesselation::GetCircleOfConnectedTriangles() - ") + std::string("AngleZero is 0 even with alternative reference.") + std::string("The algorithm has to be changed here!"));
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      ELOG(0, "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!");
+      performCriticalExit();
+    }
   }
   LOG(4, "DEBUG: Reference vector on this plane representing angle 0 is " << AngleZero << ".");
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal();
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
     OrthogonalVector.MakeNormalTo(PlaneNormal);
@@ -2927 +2854 @@
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     LOG(4, "DEBUG" << angle << " for point " << **listRunner << ".");
-    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
   }
 
@@ -2982 +2909 @@
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
-    helper = Plane(((*TesselA)->getPosition()), ((*TesselB)->getPosition()), ((*TesselC)->getPosition())).getNormal();
+    helper = Plane(((*TesselA)->getPosition()),
+                   ((*TesselB)->getPosition()),
+                   ((*TesselC)->getPosition())).getNormal();
     LOG(5, "DEBUG: Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper);
     counter++;
@@ -2991 +2920 @@
   }
   //LOG(0, "Summed vectors " << center << "; number of points " << connectedPoints.size() << "; scale factor " << counter);
-  PlaneNormal.Scale(1.0 / (double)counter);
+  PlaneNormal.Scale(1.0 / (double) counter);
   //  LOG(1, "INFO: Calculated center of all circle points is " << center << ".");
   //
@@ -3005 +2934 @@
     AngleZero.ProjectOntoPlane(PlaneNormal);
   }
-  if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON)) {
+  if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON )) {
     LOG(4, "DEBUG: Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer.");
     AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
-    ASSERT(AngleZero.NormSquared() > MYEPSILON, std::string("Tesselation::GetCircleOfSetOfPoints() - ") + std::string("AngleZero is 0 even with alternative reference.") + std::string("The algorithm has to be changed here!"));
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      ELOG(0, "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!");
+      performCriticalExit();
+    }
   }
   LOG(4, "DEBUG: Reference vector on this plane representing angle 0 is " << AngleZero << ".");
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal();
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
     OrthogonalVector.MakeNormalTo(PlaneNormal);
@@ -3027 +2959 @@
       angle = 2. * M_PI - angle;
     LOG(4, "DEBUG: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << ".");
-    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
-    ASSERT(InserterTest.second, std::string("Tesselation::GetCircleOfSetOfPoints() - ") + std::string("got two atoms with same angle " + toString(*((InserterTest.first)->second))) + std::string(" and " + toString((*listRunner))));
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
+    if (!InserterTest.second) {
+      ELOG(0, "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner));
+      performCriticalExit();
+    }
   }
 
@@ -3048 +2983 @@
   //Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPointList *> *ListOfPaths = new list<TesselPointList *>;
+  list<TesselPointList *> *ListOfPaths = new list<TesselPointList *> ;
   TesselPointList *connectedPath = NULL;
   Vector center;
@@ -3074 +3009 @@
   map<class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
   for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
-    LOG(4, "DEBUG: Adding " << *Runner->second << " to TouchedLine map.");
-    TouchedLine.insert(pair<class BoundaryLineSet *, bool>(Runner->second, false));
-    for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++) {
-      LOG(4, "DEBUG: Adding " << *Sprinter->second << " to TouchedTriangle map.");
-      TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool>(Sprinter->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));
   }
   if (!ReferencePoint->lines.empty()) {
@@ -3086 +3018 @@
       if (LineRunner == TouchedLine.end()) {
         ELOG(1, "I could not find " << *runner->second << " in the touched list.");
-      } else
-        if (!LineRunner->second) {
-          LineRunner->second = true;
-          connectedPath = new TesselPointList;
-          triangle = NULL;
-          CurrentLine = runner->second;
-          StartLine = CurrentLine;
-          CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
-          LOG(3, "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << ".");
-          do {
-            // push current one
-            LOG(3, "INFO: Putting " << *CurrentPoint << " at end of path.");
-            connectedPath->push_back(CurrentPoint->node);
-
-            // find next triangle
-            for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
-              LOG(4, "DEBUG: Inspecting triangle " << *Runner->second << ".");
-              if ((Runner->second != triangle)) { // look for first triangle not equal to old one
-                triangle = Runner->second;
-                TriangleRunner = TouchedTriangle.find(triangle);
-                if (TriangleRunner != TouchedTriangle.end()) {
-                  if (!TriangleRunner->second) {
-                    TriangleRunner->second = true;
-                    LOG(4, "DEBUG: Connecting triangle is " << *triangle << ".");
-                    break;
-                  } else {
-                    LOG(4, "DEBUG: Skipping " << *triangle << ", as we have already visited it.");
-                    triangle = NULL;
-                  }
+      } else if (!LineRunner->second) {
+        LineRunner->second = true;
+        connectedPath = new TesselPointList;
+        triangle = NULL;
+        CurrentLine = runner->second;
+        StartLine = CurrentLine;
+        CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+        LOG(1, "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << ".");
+        do {
+          // push current one
+          LOG(1, "INFO: Putting " << *CurrentPoint << " at end of path.");
+          connectedPath->push_back(CurrentPoint->node);
+
+          // find next triangle
+          for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
+            LOG(1, "INFO: Inspecting triangle " << *Runner->second << ".");
+            if ((Runner->second != triangle)) { // look for first triangle not equal to old one
+              triangle = Runner->second;
+              TriangleRunner = TouchedTriangle.find(triangle);
+              if (TriangleRunner != TouchedTriangle.end()) {
+                if (!TriangleRunner->second) {
+                  TriangleRunner->second = true;
+                  LOG(1, "INFO: Connecting triangle is " << *triangle << ".");
+                  break;
                 } else {
-                  ELOG(1, "I could not find " << *triangle << " in the touched list.");
+                  LOG(1, "INFO: Skipping " << *triangle << ", as we have already visited it.");
                   triangle = NULL;
                 }
               } else {
-                // as we have stumbled upon the same triangle, we don't need the check anymore
+                ELOG(1, "I could not find " << *triangle << " in the touched list.");
                 triangle = NULL;
               }
             }
-            if (triangle == NULL)
+          }
+          if (triangle == NULL)
+            break;
+          // find next line
+          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];
+              LOG(1, "INFO: Connecting line is " << *CurrentLine << ".");
               break;
-            // find next line
-            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];
-                LOG(3, "INFO: Connecting line is " << *CurrentLine << ".");
-                break;
-              }
             }
-            LineRunner = TouchedLine.find(CurrentLine);
-            if (LineRunner == TouchedLine.end())
-              ELOG(1, "I could not find " << *CurrentLine << " in the touched list.");
-            else
-              LineRunner->second = true;
-            // find next point
-            CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
-
-          } while (CurrentLine != StartLine);
-          // last point is missing, as it's on start line
-          if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back()) {
-            LOG(3, "INFO: Putting " << *CurrentPoint << " at end of path to close it.");
-            connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
           }
-
-          ListOfPaths->push_back(connectedPath);
-        } else {
-          LOG(3, "DEBUG: Skipping " << *runner->second << ", as we have already visited it.");
-        }
+          LineRunner = TouchedLine.find(CurrentLine);
+          if (LineRunner == TouchedLine.end())
+            ELOG(1, "I could not find " << *CurrentLine << " in the touched list.");
+          else
+            LineRunner->second = true;
+          // find next point
+          CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+
+        } while (CurrentLine != StartLine);
+        // last point is missing, as it's on start line
+        LOG(1, "INFO: Putting " << *CurrentPoint << " at end of path.");
+        if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
+          connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+
+        ListOfPaths->push_back(connectedPath);
+      } else {
+        LOG(1, "INFO: Skipping " << *runner->second << ", as we have already visited it.");
+      }
     }
   } else {
@@ -3171 +3098 @@
   //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;
@@ -3181 +3108 @@
     connectedPath = *ListRunner;
 
-    if (DoLog(2)) {
-      std::stringstream output;
-      output << "INFO: Current path is ";
-      BOOST_FOREACH(const TesselPoint * const item, *connectedPath) {
-        output << *item << " ";
-      }
-      LOG(1, output.str());
-    }
+    LOG(1, "INFO: Current path is " << connectedPath << ".");
 
     // go through list, look for reappearance of starting Point and count
@@ -3197 +3117 @@
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
-        if (DoLog(3)) {
+        if (DoLog(1)) {
           std::stringstream output;
-          output << "DEBUG: " << count + 1 << ". closed path consists of: ";
-          for (TesselPointList::iterator CircleSprinter = Marker; CircleSprinter != CircleRunner; CircleSprinter++)
+          output << count + 1 << ". closed path consists of: ";
+          for (TesselPointList::iterator CircleSprinter = Marker;
+              CircleSprinter != CircleRunner;
+              CircleSprinter++)
             output << (**CircleSprinter) << " <-> ";
           LOG(1, output.str());
@@ -3216 +3138 @@
     }
   }
-  LOG(2, "DEBUG: " << count << " closed additional path(s) have been created.");
+  LOG(1, "INFO: " << count << " closed additional path(s) have been created.");
 
   // delete list of paths
@@ -3254 +3176 @@
 }
 ;
-
-struct CloserToPiHalf
-{
-  bool operator()(double angle, double smallestangle)
-  {
-    return (fabs(angle - M_PI / 2.) < fabs(smallestangle - M_PI / 2.));
-  }
-};
-
-bool Tesselation::IsPointBelowSurroundingPolygon(const BoundaryPointSet *_point) const
-{
-  // check for NULL
-  if (_point == NULL) {
-    return false;
-  }
-
-  // get list of connected points
-  if (_point->lines.empty()) {
-    LOG(1, "INFO: point " << *_point << " is not connected to any lines.");
-    return false;
-  }
-  bool PointIsBelow = true;
-
-  // create Orthogonal vector as reference for angle (pointing into [pi,2pi) interval)
-  Vector OrthogonalVector;
-  for (LineMap::const_iterator lineiter = _point->lines.begin();
-      lineiter != _point->lines.end(); ++lineiter)
-    for (TriangleMap::const_iterator triangleiter = lineiter->second->triangles.begin();
-        triangleiter != lineiter->second->triangles.end();
-        ++triangleiter)
-        OrthogonalVector += triangleiter->second->NormalVector;
-  OrthogonalVector.Normalize();
-
-  // go through all closed paths for this point
-  typedef list<TesselPointList *> TPL_list_t;
-  const TPL_list_t *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(_point->node);
-  for (TPL_list_t::const_iterator closedpathsiter = ListOfClosedPaths->begin();
-      (closedpathsiter != ListOfClosedPaths->end()) && PointIsBelow;
-      ++closedpathsiter) {
-    const TesselPointList *connectedPath = *closedpathsiter;
-
-    TesselPointList::const_iterator ListAdvance = connectedPath->begin(); // gives angle direction
-    TesselPointList::const_iterator ListRunner = ListAdvance++;
-    for (; (ListAdvance != connectedPath->end()) && PointIsBelow;
-        ListRunner = ListAdvance++) { // go through all closed paths
-      LOG(2, "DEBUG: Current reference node is " << **ListRunner
-          << ", advanced node is " << **ListAdvance);
-
-      // reference vector to point to check for this point of connected path
-      const Vector Reference =
-          ((*ListAdvance)->getPosition()) - ((*ListRunner)->getPosition());
-
-      // go through all other points in this connected path
-      for (TesselPointList::const_iterator OtherRunner = ListAdvance;
-          (OtherRunner != ListRunner) && PointIsBelow;
-          ++OtherRunner == connectedPath->end() ?
-              OtherRunner = connectedPath->begin() :
-              OtherRunner) {
-        if (OtherRunner == ListAdvance)
-          continue;
-        LOG(3, "DEBUG: Current other node is " << **OtherRunner);
-
-        // build the plane with normal vector
-        const Vector onebeam =
-            ((*OtherRunner)->getPosition()) - ((*ListAdvance)->getPosition());
-        Vector NormalVector = Reference;
-        NormalVector.VectorProduct(onebeam);
-        // needs to point in same general direction as average NormalVector of all triangles
-        if (NormalVector.ScalarProduct(OrthogonalVector) < 0)
-          NormalVector *= -1.;
-        try {
-          Plane plane(NormalVector, ((*ListRunner)->getPosition()));
-
-          // check whether point is below
-          if (plane.SignedDistance(_point->node->getPosition()) > 0) {
-            LOG(2, "DEBUG: For plane " << plane << " point " << *_point
-                << " would remain above.");
-            PointIsBelow = false;
-          }
-        } catch (ZeroVectorException &e) {
-          ELOG(3, "Vectors are linear dependent, skipping.");
-        }
-      }
-    }
-  }
-  delete ListOfClosedPaths;
-
-  return PointIsBelow;
-}
-
-double Tesselation::RemovePointSurroundedByPolygon(
-    TesselPointList *connectedPath,
-    BoundaryPointSet *point)
-{
-  double volume = 0.;
-  const Vector OldPoint = point->node->getPosition();
-
-  TesselPoint *oldNode = point->node;
-  // remove present triangles for this connectedPath
-  unsigned int count = 0;
-  for (TesselPointList::const_iterator iter = connectedPath->begin();
-      iter != connectedPath->end(); ++iter)
-    LOG(4, "DEBUG: Node in connectedPath is " << **iter);
-
-  {
-    TesselPointList::iterator FirstNode, SecondNode;
-    SecondNode = connectedPath->begin();
-    FirstNode = SecondNode++;
-    for (;FirstNode != connectedPath->end(); ++SecondNode, ++FirstNode) {
-      LOG(3, "DEBUG: MiddleNode is " << **FirstNode << ".");
-      if (SecondNode == connectedPath->end())
-        SecondNode = connectedPath->begin();
-      TesselPoint *TriangleCandidates[3];
-      TriangleCandidates[0] = *FirstNode;
-      TriangleCandidates[1] = *SecondNode;
-      TriangleCandidates[2] = oldNode;
-      BoundaryTriangleSet *triangle = GetPresentTriangle(TriangleCandidates);
-      ASSERT( triangle != NULL,
-          "Tesselation::RemovePointSurroundedByPolygon() - triangle to points "
-          +toString((**SecondNode))+", "+toString((**FirstNode))+" and "
-          +toString(*oldNode)+" does not exist.");
-      LOG(3, "DEBUG: Removing triangle " << *triangle << ".");
-      RemoveTesselationTriangle(triangle);
-      ++count;
-    }
-    LOG(2, "INFO: " << count << " triangles were removed.");
-  }
-
-  // re-create all triangles by going through connected points list
-  LineList NewLines;
-  typedef std::vector<double> angles_t;
-  angles_t angles;
-  count = 0;
-  for (; !connectedPath->empty();) {
-    // search middle node with widest angle to next neighbours
-    TesselPointList::iterator StartNode, MiddleNode, EndNode;
-    for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
-      LOG(3, "INFO: MiddleNode is " << **MiddleNode << ".");
-      // construct vectors to next and previous neighbour
-      StartNode = MiddleNode;
-      if (StartNode == connectedPath->begin())
-        StartNode = connectedPath->end();
-      StartNode--;
-      //LOG(3, "INFO: StartNode is " << **StartNode << ".");
-      const Vector Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
-      EndNode = MiddleNode;
-      EndNode++;
-      if (EndNode == connectedPath->end())
-        EndNode = connectedPath->begin();
-      //LOG(3, "INFO: EndNode is " << **StartNode << ".");
-      const Vector Reference = ((*EndNode)->getPosition()) - ((*MiddleNode)->getPosition());
-      Vector OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
-      OrthogonalVector.MakeNormalTo(Reference);
-      angles.push_back(GetAngle(Point, Reference, OrthogonalVector));
-    }
-    ASSERT( !angles.empty(),
-        "Tesselation::RemovePointSurroundedByPolygon() - angles empty");
-    const angles_t::const_iterator maxiter = std::max_element(angles.begin(), angles.end());
-    angles_t::const_iterator miniter = angles.begin();
-    // distinguish between convex and nonconvex polygon
-    if (*maxiter > M_PI) {
-      // connectedPath is not convex: The idea is to fill any kinks pointing
-      // inside into the connectedPath close to this concave spot first, making
-      // it eventually become convex.
-      // Hence, use adjacent (and convex) fill-in point
-      miniter = (maxiter == angles.begin()) ? angles.end()-1 : maxiter-1;
-      if (*miniter > M_PI) {
-        miniter = (maxiter+1 == angles.end()) ? angles.begin() : maxiter+1;
-        if (*miniter > M_PI) {
-          miniter = std::min_element(angles.begin(), angles.end());
-        }
-      }
-    } else {
-      // is convex
-      miniter = std::min_element(angles.begin(), angles.end(), CloserToPiHalf());
-    }
-    MiddleNode = connectedPath->begin();
-    std::advance(MiddleNode, std::distance(const_cast<const angles_t &>(angles).begin(), miniter));
-
-    ASSERT(MiddleNode != connectedPath->end(),
-        "Tesselation::RemovePointSurroundedByPolygon() - Could not find a smallest angle!");
-    angles.clear();
-    StartNode = MiddleNode;
-    EndNode = MiddleNode;
-    if (StartNode == connectedPath->begin())
-      StartNode = connectedPath->end();
-    StartNode--;
-    EndNode++;
-    if (EndNode == connectedPath->end())
-      EndNode = connectedPath->begin();
-    LOG(2, "INFO: StartNode is " << **StartNode << ".");
-    LOG(2, "INFO: MiddleNode is " << **MiddleNode << ".");
-    LOG(2, "INFO: EndNode is " << **EndNode << ".");
-    LOG(1, "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << ".");
-    TesselPoint *TriangleCandidates[3];
-    TriangleCandidates[0] = *StartNode;
-    TriangleCandidates[1] = *MiddleNode;
-    TriangleCandidates[2] = *EndNode;
-    BoundaryTriangleSet *triangle = GetPresentTriangle(TriangleCandidates);
-    if (triangle != NULL) {
-      const Vector center = 1./3. * ((*StartNode)->getPosition()
-              + (*MiddleNode)->getPosition()
-              + (*EndNode)->getPosition());
-      const Vector NormalVector = OldPoint - center;
-      // check orientation of normal vector (that points inside)
-      ASSERT( triangle->NormalVector.ScalarProduct(NormalVector) > std::numeric_limits<double>::epsilon()*1e2,
-          "Tesselation::RemovePointSurroundedByPolygon() - New triangle with same orientation already present as "
-          +toString(*triangle)+"!");
-    }
-
-    LOG(3, "DEBUG: Adding new triangle points.");
-    AddTesselationPoint(*StartNode, 0);
-    AddTesselationPoint(*MiddleNode, 1);
-    AddTesselationPoint(*EndNode, 2);
-    LOG(3, "DEBUG: Adding new triangle lines.");
-    AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
-    // line between start and end must be new (except for very last triangle)
-    if (AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1))
-      NewLines.push_back(BLS[1]);
-    AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
-    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-    const Vector center = 1./3. * ((*StartNode)->getPosition()
-            + (*MiddleNode)->getPosition()
-            + (*EndNode)->getPosition());
-    const Vector NormalVector = OldPoint - center;
-    BTS->GetNormalVector(NormalVector);
-    AddTesselationTriangle();
-    // calculate volume summand as a general tetraeder
-    volume += CalculateVolumeofGeneralTetraeder(
-        TPS[0]->node->getPosition(),
-        TPS[1]->node->getPosition(),
-        TPS[2]->node->getPosition(),
-        OldPoint);
-    // advance number
-    ++count;
-
-    // prepare nodes for next triangle
-    LOG(3, "DEBUG: Removing " << **MiddleNode << " from closed path.");
-    connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
-    LOG(3, "DEBUG: Remaining points: " << connectedPath->size() << ".");
-    ASSERT(connectedPath->size() >= 2,
-        "Tesselation::RemovePointSurroundedByPolygon() - There are only two endpoints left!");
-    if (connectedPath->size() == 2) { // we are done
-      connectedPath->remove(*StartNode); // remove the start node
-      connectedPath->remove(*EndNode); // remove the end node
-    }
-  }
-  LOG(1, "INFO: " << count << " triangles were created.");
-
-  return volume;
-}
 
 /** Removes a boundary point from the envelope while keeping it closed.
@@ -3519 +3190 @@
 double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point)
 {
+  class BoundaryLineSet *line = NULL;
+  class BoundaryTriangleSet *triangle = NULL;
+  Vector OldPoint, NormalVector;
   double volume = 0;
+  int count = 0;
 
   if (point == NULL) {
@@ -3527 +3202 @@
     LOG(4, "DEBUG: Removing point " << *point << " from tesselated boundary ...");
 
+  // copy old location for the volume
+  OldPoint = (point->node->getPosition());
+
   // get list of connected points
   if (point->lines.empty()) {
@@ -3534 +3212 @@
 
   list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
-  list<TesselPointList *>::iterator ListRunner = ListAdvance;
-//  TriangleMap::iterator NumberRunner = Candidates.begin();
-  Vector Point, Reference, OrthogonalVector;
-  for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
-    if (ListAdvance != ListOfClosedPaths->end())
-      ListAdvance++;
-
-    TesselPointList *connectedPath = *ListRunner;
-
-    volume += RemovePointSurroundedByPolygon(connectedPath, point);
-
-    ListOfClosedPaths->remove(connectedPath);
-    delete (connectedPath);
-  }
-  delete (ListOfClosedPaths);
-
-  LOG(1, "INFO: Removed volume is " << volume << ".");
-
-  return volume;
-}
-;
-
-bool Tesselation::CheckAllConcaveInPolygon(
-    const TesselPointList *connectedPath,
-    const BoundaryPointSet *point
-    )
-{
-  // check whether lines in this path to point to remove are all concave
-  bool all_lines_concave = true;
-  if (connectedPath->size() >= 2) {
-    TesselPointList::const_iterator StartNode, MiddleNode, EndNode;
-    // have nearest neighbors to a middle node to know adjacent triangles
-    for (MiddleNode = connectedPath->begin();
-        all_lines_concave && (MiddleNode != connectedPath->end());
-        MiddleNode++) {
-      LOG(3, "INFO: MiddleNode is " << **MiddleNode << ".");
-      EndNode = MiddleNode;
-      if (EndNode == connectedPath->begin())
-        EndNode = connectedPath->end();
-      --EndNode;
-      StartNode = MiddleNode;
-      ++StartNode;
-      if (StartNode == connectedPath->end())
-        StartNode = connectedPath->begin();
-
-      AddTesselationPoint(point->node, 0);
-      AddTesselationPoint(*MiddleNode, 1);
-
-      ASSERT( point->lines.find((*MiddleNode)->getNr()) != point->lines.end(),
-          "Tesselation::CheckAllConcaveInPolygon() - MiddleNode "
-          +toString(**MiddleNode)+" not present in "
-          +toString(*point)+"'s lines.");
-
-      // get line between Node and point and check
-      std::pair<LineMap::const_iterator, LineMap::const_iterator> FindPair =
-          point->lines.equal_range((*MiddleNode)->getNr());
-      LineMap::const_iterator FindLine = FindPair.first;
-      for (; FindLine != FindPair.second; ++FindLine) {
-        // line  has got two triangles, check whether they resemble those
-        // with start and endnode
-        const BoundaryLineSet *currentline = FindLine->second;
-        unsigned int matching_triangles = 0;
-        for (TriangleMap::const_iterator triangleiter = currentline->triangles.begin();
-            triangleiter != currentline->triangles.end(); ++triangleiter) {
-          const BoundaryTriangleSet *triangle = triangleiter->second;
-          AddTesselationPoint(*StartNode, 2);
-          if (triangle->IsPresentTupel(TPS))
-            ++matching_triangles;
-          AddTesselationPoint(*EndNode, 2);
-          if (triangle->IsPresentTupel(TPS))
-            ++matching_triangles;
-        }
-        if (matching_triangles == 2)
-          break;
-      }
-      if (FindLine != FindPair.second) {
-        LOG(3, "INFO: Current line of point " << *point << " is " << *FindLine->second << ".");
-        all_lines_concave &= !FindLine->second->CheckConvexityCriterion();
-      }
-    }
-  } else {
-    // check the single line
-    if (connectedPath->empty())
-      return false;
-    LineMap::const_iterator FindLine = point->lines.find((*connectedPath->begin())->getNr());
-    ASSERT( FindLine != point->lines.end(),
-        "Tesselation::CheckAllConcaveInPolygon() - point "
-        +toString((*connectedPath->begin())->getNr())+" not present in "
-        +toString(*point)+"'s lines.");
-    return !FindLine->second->CheckConvexityCriterion();
-  }
-
-  return all_lines_concave;
-}
-
-/** Removes a boundary point from the envelope while keeping it closed.
- * We remove the old triangles connected to the point and re-create new triangles to close the surface following this ansatz:
- *  -# a closed path(s) of boundary points surrounding the point to be removed is constructed
- *  -# on each closed path, we pick three adjacent points, create a triangle with them and subtract the middle point from the path
- *  -# we advance two points (i.e. the next triangle will start at the ending point of the last triangle) and continue as before
- *  -# the surface is closed, when the path is empty
- * Thereby, we (hopefully) make sure that the removed points remains beneath the surface (this is checked via IsInnerPoint eventually).
- * \param *out output stream for debugging
- * \param *point point to be removed
- * \return volume added to the volume inside the tesselated surface by the removal
- */
-double Tesselation::RemoveFullConcavePointFromTesselatedSurface(class BoundaryPointSet *point)
-{
-  double volume = 0;
-
-  if (point == NULL) {
-    ELOG(1, "Cannot remove the point " << point << ", it's NULL!");
-    return 0.;
-  } else
-    LOG(4, "DEBUG: Removing point " << *point << " from tesselated boundary ...");
-
-  // get list of connected points
-  if (point->lines.empty()) {
-    ELOG(1, "Cannot remove the point " << *point << ", it's connected to no lines!");
-    return 0.;
-  }
-
-  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
+
+  // gather all triangles
+  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
+    count += LineRunner->second->triangles.size();
+  TriangleMap Candidates;
+  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+    line = LineRunner->second;
+    for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
+      triangle = TriangleRunner->second;
+      Candidates.insert(TrianglePair(triangle->Nr, triangle));
+    }
+  }
+
+  // remove all triangles
+  count = 0;
+  NormalVector.Zero();
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    LOG(1, "INFO: Removing triangle " << *(Runner->second) << ".");
+    NormalVector -= Runner->second->NormalVector; // has to point inward
+    RemoveTesselationTriangle(Runner->second);
+    count++;
+  }
+  LOG(1, count << " triangles were removed.");
+
   list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
   list<TesselPointList *>::iterator ListRunner = ListAdvance;
 //  TriangleMap::iterator NumberRunner = Candidates.begin();
   TesselPointList::iterator StartNode, MiddleNode, EndNode;
+  double angle;
+  double smallestangle;
   Vector Point, Reference, OrthogonalVector;
-  for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
-    if (ListAdvance != ListOfClosedPaths->end())
-      ListAdvance++;
-
-    TesselPointList *connectedPath = *ListRunner;
-
-    if (CheckAllConcaveInPolygon(connectedPath, point)) {
-      LOG(1, "INFO: ... point " << *point << " cannot be on convex envelope, all lines concave.");
-      volume += RemovePointSurroundedByPolygon(connectedPath, point);
-    }
-
-    ListOfClosedPaths->remove(connectedPath);
-    delete (connectedPath);
+  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
+      if (ListAdvance != ListOfClosedPaths->end())
+        ListAdvance++;
+
+      connectedPath = *ListRunner;
+      // re-create all triangles by going through connected points list
+      LineList NewLines;
+      for (; !connectedPath->empty();) {
+        // search middle node with widest angle to next neighbours
+        EndNode = connectedPath->end();
+        smallestangle = 0.;
+        for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
+          LOG(1, "INFO: MiddleNode is " << **MiddleNode << ".");
+          // construct vectors to next and previous neighbour
+          StartNode = MiddleNode;
+          if (StartNode == connectedPath->begin())
+            StartNode = connectedPath->end();
+          StartNode--;
+          //LOG(3, "INFO: StartNode is " << **StartNode << ".");
+          Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+          StartNode = MiddleNode;
+          StartNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          //LOG(3, "INFO: EndNode is " << **StartNode << ".");
+          Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+          OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
+          OrthogonalVector.MakeNormalTo(Reference);
+          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;
+          }
+        }
+        MiddleNode = EndNode;
+        if (MiddleNode == connectedPath->end()) {
+          ELOG(0, "CRITICAL: Could not find a smallest angle!");
+          performCriticalExit();
+        }
+        StartNode = MiddleNode;
+        if (StartNode == connectedPath->begin())
+          StartNode = connectedPath->end();
+        StartNode--;
+        EndNode++;
+        if (EndNode == connectedPath->end())
+          EndNode = connectedPath->begin();
+        LOG(2, "INFO: StartNode is " << **StartNode << ".");
+        LOG(2, "INFO: MiddleNode is " << **MiddleNode << ".");
+        LOG(2, "INFO: EndNode is " << **EndNode << ".");
+        LOG(1, "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << ".");
+        TriangleCandidates[0] = *StartNode;
+        TriangleCandidates[1] = *MiddleNode;
+        TriangleCandidates[2] = *EndNode;
+        triangle = GetPresentTriangle(TriangleCandidates);
+        if (triangle != NULL) {
+          ELOG(0, "New triangle already present, skipping!");
+          StartNode++;
+          MiddleNode++;
+          EndNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+          continue;
+        }
+        LOG(3, "Adding new triangle points.");
+        AddTesselationPoint(*StartNode, 0);
+        AddTesselationPoint(*MiddleNode, 1);
+        AddTesselationPoint(*EndNode, 2);
+        LOG(3, "Adding new triangle lines.");
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+        NewLines.push_back(BLS[1]);
+        AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
+        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+        BTS->GetNormalVector(NormalVector);
+        AddTesselationTriangle();
+        // calculate volume summand as a general tetraeder
+        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
+        // advance number
+        count++;
+
+        // prepare nodes for next triangle
+        StartNode = EndNode;
+        LOG(2, "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << ".");
+        connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
+        if (connectedPath->size() == 2) { // we are done
+          connectedPath->remove(*StartNode); // remove the start node
+          connectedPath->remove(*EndNode); // remove the end node
+          break;
+        } else if (connectedPath->size() < 2) { // something's gone wrong!
+          ELOG(0, "CRITICAL: There are only two endpoints left!");
+          performCriticalExit();
+        } else {
+          MiddleNode = StartNode;
+          MiddleNode++;
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          EndNode = MiddleNode;
+          EndNode++;
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+        }
+      }
+      // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
+      if (NewLines.size() > 1) {
+        LineList::iterator Candidate;
+        class BoundaryLineSet *OtherBase = NULL;
+        double tmp, maxgain;
+        do {
+          maxgain = 0;
+          for (LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+            tmp = PickFarthestofTwoBaselines(*Runner);
+            if (maxgain < tmp) {
+              maxgain = tmp;
+              Candidate = Runner;
+            }
+          }
+          if (maxgain != 0) {
+            volume += maxgain;
+            LOG(1, "Flipping baseline with highest volume" << **Candidate << ".");
+            OtherBase = FlipBaseline(*Candidate);
+            NewLines.erase(Candidate);
+            NewLines.push_back(OtherBase);
+          }
+        } while (maxgain != 0.);
+      }
+
+      ListOfClosedPaths->remove(connectedPath);
+      delete (connectedPath);
+    }
+    LOG(1, "INFO: " << count << " triangles were created.");
+  } else {
+    while (!ListOfClosedPaths->empty()) {
+      ListRunner = ListOfClosedPaths->begin();
+      connectedPath = *ListRunner;
+      ListOfClosedPaths->remove(connectedPath);
+      delete (connectedPath);
+    }
+    LOG(3, "DEBUG: No need to create any triangles.");
   }
   delete (ListOfClosedPaths);
 
-  if (volume > 0.)
-    LOG(1, "INFO: Removed volume is " << volume << ".");
+  LOG(1, "INFO: Removed volume is " << volume << ".");
 
   return volume;
@@ -3724 +3437 @@
             if (TrianglePoints[j] != NULL) {
               for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->getNr()); // is a multimap!
-                  (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->getNr()); FindLine++) {
+              (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->getNr()); FindLine++) {
                 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
                   if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
@@ -3745 +3458 @@
           break;
       for (FindLine = TrianglePoints[(i + 1) % 3]->lines.find(TrianglePoints[(i + 2) % 3]->node->getNr()); // is a multimap!
-          (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->getNr()); FindLine++) {
+      (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->getNr()); FindLine++) {
         for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
           if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
@@ -3774 +3487 @@
     }
     default:
-      ASSERT(0, "Tesselation::FindTriangles() - Number of wildcards is greater than 3, cannot happen!");
+      ELOG(0, "Number of wildcards is greater than 3, cannot happen!");
+      performCriticalExit();
       break;
   }
@@ -3800 +3514 @@
     if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])
       return true;
-    else
-      if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
+    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;
   }
@@ -3830 +3542 @@
 
   // sanity check
-      if (LinesOnBoundary.empty()) {
-        ELOG(2, "FindAllDegeneratedTriangles() was called without any tesselation structure.");
-        return DegeneratedLines;
-      }
-      LineMap::iterator LineRunner1;
-      pair<UniqueLines::iterator, bool> tester;
+  if (LinesOnBoundary.empty()) {
+    ELOG(2, "FindAllDegeneratedTriangles() was called without any tesselation structure.");
+    return DegeneratedLines;
+  }
+  LineMap::iterator LineRunner1;
+  pair<UniqueLines::iterator, bool> tester;
   for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
     tester = AllLines.insert(LineRunner1->second);
@@ -3852 +3564 @@
     const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second);
     if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end())
-    LOG(3, "DEBUG: " << *Line1->second << " => " << *Line2->second);
+      LOG(3, "DEBUG: " << *Line1->second << " => " << *Line2->second);
     else
-    ELOG(1, "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!");
+      ELOG(1, "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!");
   }
 
@@ -3886 +3598 @@
       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));
+          DegeneratedTriangles->insert(pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr));
+          DegeneratedTriangles->insert(pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr));
         }
       }
@@ -3914 +3626 @@
 
   // iterate over all degenerated triangles
-  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
-      !DegeneratedTriangles->empty();
-      TriangleKeyRunner = DegeneratedTriangles->begin()) {
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); !DegeneratedTriangles->empty(); TriangleKeyRunner = DegeneratedTriangles->begin()) {
     LOG(3, "DEBUG: Checking presence of triangles " << TriangleKeyRunner->first << " and " << TriangleKeyRunner->second << ".");
     // both ways are stored in the map, only use one
@@ -3957 +3667 @@
 //                OtherpartnerTriangle = TriangleRunner->second;
 //              }
-              /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]
-              // the line of triangle receives the degenerated ones
-              triangle->lines[i]->triangles.erase(Othertriangle->Nr);
+            /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]
+            // 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++)
@@ -3973 +3683 @@
 
       // erase the pair
-      count += (int)DegeneratedTriangles->erase(triangle->Nr);
+      count += (int) DegeneratedTriangles->erase(triangle->Nr);
       LOG(4, "DEBUG: RemoveDegeneratedTriangles() removes triangle " << *triangle << ".");
       RemoveTesselationTriangle(triangle);
-      count += (int)DegeneratedTriangles->erase(partnerTriangle->Nr);
+      count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
       LOG(4, "DEBUG: RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << ".");
       RemoveTesselationTriangle(partnerTriangle);
     } else {
       LOG(4, "DEBUG: 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.");
-      // we need to remove them from the list nonetheless
-      DegeneratedTriangles->erase(triangle->Nr);
-      DegeneratedTriangles->erase(partnerTriangle->Nr);
     }
   }
@@ -4027 +3734 @@
   class BoundaryLineSet *BestLine = NULL;
   for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
-    BaseLine = (Runner->second->endpoints[0]->node->getPosition()) - (Runner->second->endpoints[1]->node->getPosition());
-    CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) + (Runner->second->endpoints[1]->node->getPosition()));
+    BaseLine = (Runner->second->endpoints[0]->node->getPosition()) -
+               (Runner->second->endpoints[1]->node->getPosition());
+    CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) +
+                           (Runner->second->endpoints[1]->node->getPosition()));
     CenterToPoint -= (point->getPosition());
     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;
@@ -4080 +3789 @@
   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]))) {
-      ASSERT(BestLine != BTS->lines[i], std::string("Tesselation::AddBoundaryPointByDegeneratedTriangle() - ") + std::string("BestLine is same as found line, something's wrong here!"));
-      BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *>(TempTriangle->Nr, TempTriangle));
+      if (BestLine == BTS->lines[i]) {
+        ELOG(0, "BestLine is same as found line, something's wrong here!");
+        performCriticalExit();
+      }
+      BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle));
       TempTriangle->lines[nr] = BTS->lines[i];
       break;
@@ -4103 +3815 @@
   if (LastTriangle != NULL) {
     stringstream sstr;
-    sstr << "-" << TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
+    sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
     NumberName = sstr.str();
     if (DoTecplotOutput) {
@@ -4145 +3857 @@
     if (s1->endpoints.size() < s2->endpoints.size())
       return true;
-    else
-      if (s1->endpoints.size() > s2->endpoints.size())
-        return false;
-      else { // equality of number of endpoints
-        PointSet::const_iterator Walker1 = s1->endpoints.begin();
-        PointSet::const_iterator Walker2 = s2->endpoints.begin();
-        while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
-          if ((*Walker1)->Nr < (*Walker2)->Nr)
-            return true;
-          else
-            if ((*Walker1)->Nr > (*Walker2)->Nr)
-              return false;
-          Walker1++;
-          Walker2++;
-        }
-        return false;
+    else if (s1->endpoints.size() > s2->endpoints.size())
+      return false;
+    else { // equality of number of endpoints
+      PointSet::const_iterator Walker1 = s1->endpoints.begin();
+      PointSet::const_iterator Walker2 = s2->endpoints.begin();
+      while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
+        if ((*Walker1)->Nr < (*Walker2)->Nr)
+          return true;
+        else if ((*Walker1)->Nr > (*Walker2)->Nr)
+          return false;
+        Walker1++;
+        Walker2++;
       }
+      return false;
+    }
   }
 };
@@ -4186 +3896 @@
       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(5, "DEBUG:  Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list.");
@@ -4199 +3909 @@
         if (VectorWalker != VectorRunner) { // skip equals
           const double SCP = VectorWalker->second->ScalarProduct(*VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles
-          LOG(4, "DEBUG: Checking " << *(VectorWalker->second) << " against " << *(VectorRunner->second) << ": " << SCP);
+          LOG(4, "DEBUG: Checking " << *VectorWalker->second << " against " << *VectorRunner->second << ": " << SCP);
           if (fabs(SCP + 1.) < ParallelEpsilon) {
             InsertionTester = EndpointCandidateList.insert((Runner->second));
@@ -4281 +3991 @@
     // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
     // connections to either polygon ...
-    ASSERT(T->size() % 2 == 0, std::string("Tesselation::CorrectAllDegeneratedPolygons() - ") + std::string(" degenerated polygon contains an odd number of triangles,") + std::string(" probably contains bridging non-degenerated ones, too!"));
+    if (T->size() % 2 != 0) {
+      ELOG(0, " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!");
+      performCriticalExit();
+    }
     TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator
     /// 4a. Get NormalVector for one side (this is "front")