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

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>
 
@@ -66 +69 @@
 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;
 
@@ -76 +79 @@
  */
 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__);
@@ -113 +111 @@
 {
   // 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();
@@ -126 +124 @@
       return;
     }
-        }
-
-        // find a starting triangle
+  }
+
+  // find a starting triangle
   FindStartingTriangle(SPHERERADIUS, LinkedList);
 
@@ -142 +140 @@
         //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) {
@@ -155 +152 @@
       }
     }
-    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+    if ((ShortestAngle == 4. * M_PI) || (baseline->pointlist.empty()))
       OneLoopWithoutSuccessFlag = false;
     else {
@@ -161 +158 @@
     }
   }
+
+  delete LinkedList;
 }
 
@@ -170 +169 @@
 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();
+  }
+
+  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
   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
-  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.");
+  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.");
 
   return volume;
@@ -207 +208 @@
 
   // 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.
@@ -251 +248 @@
         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)));
       }
     }
@@ -271 +268 @@
     // 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
@@ -388 +383 @@
         // 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());
 
@@ -407 +401 @@
         // 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) << ".");
@@ -460 +454 @@
             // 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
@@ -477 +467 @@
               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
+              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 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
@@ -562 +553 @@
 
   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
@@ -729 +720 @@
  * @param *b second endpoint
  * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
+ * @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)
 {
   bool insertNewLine = true;
@@ -746 +738 @@
       if (FindLine->second->triangles.size() == 1) {
         CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
-        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 << ".");
+        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.");
             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;
         }
       }
@@ -770 +772 @@
     AddExistingTesselationTriangleLine(WinningLine, n);
   }
+
+  return insertNewLine;
 }
 ;
@@ -794 +798 @@
   // 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));
 }
 ;
@@ -883 +887 @@
       } 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());
@@ -919 +923 @@
   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++) {
@@ -937 +948 @@
         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
@@ -985 +997 @@
   //Info FunctionInfo(__func__);
 
-  LOG(3, "DEBUG: Checking whether sphere contains no others points ...");
+  LOG(3, "DEBUG: Checking degeneracy by whether sphere contains no others points ...");
   bool flag = true;
 
@@ -992 +1004 @@
   TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);
 
-  LOG(3, "DEBUG: The following atoms are inside sphere at " << 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 << ":");
   for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-    LOG(3, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
+    LOG(4, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
 
   // remove triangles's endpoints
@@ -1008 +1021 @@
     LOG(3, "DEBUG: CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere.");
     flag = false;
-    LOG(3, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":");
+    LOG(4, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":");
     for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      LOG(3, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
-  }
-  delete (ListofPoints);
+      LOG(4, "DEBUG:   " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");
+  }
+  delete ListofPoints;
 
   return flag;
@@ -1155 +1168 @@
 
   // 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
@@ -1176 +1190 @@
   }
 
-  if (DoLog(1)) {
+  if (DoLog(3)) {
     std::stringstream output;
     output << "Found maximum coordinates: ";
@@ -1190 +1204 @@
     BaseLine = new BoundaryLineSet();
     BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
-    LOG(2, "DEBUG: Coordinates of start node at " << *BaseLine->endpoints[0]->node << ".");
+    LOG(2, "INFO: Coordinates of start node at " << *BaseLine->endpoints[0]->node << ".");
 
     double ShortestAngle;
@@ -1203 +1217 @@
     }
     BaseLine->endpoints[1] = new BoundaryPointSet(Temporary);
-    LOG(1, "INFO: Second node is at " << *Temporary << ".");
+    LOG(2, "INFO: Second node is at " << *Temporary << ".");
 
     // construct center of circle
     CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
-    LOG(1, "INFO: CircleCenter is at " << CircleCenter << ".");
+    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ".");
 
     // construct normal vector of circle
     CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
-    LOG(1, "INFO: CirclePlaneNormal is at " << CirclePlaneNormal << ".");
+    LOG(4, "DEBUG: CirclePlaneNormal is at " << CirclePlaneNormal << ".");
 
     double radius = CirclePlaneNormal.NormSquared();
@@ -1218 +1232 @@
     NormalVector.ProjectOntoPlane(CirclePlaneNormal);
     NormalVector.Normalize();
-    LOG(1, "INFO: NormalVector is at " << NormalVector << ".");
+    LOG(4, "DEBUG: 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) {
-      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 &e) {
+      ELOG(1, "Vectors are linear dependent: " << CirclePlaneNormal << ", " << NormalVector << ".");
       delete BaseLine;
       continue;
@@ -1235 +1248 @@
 
     // adding point 1 and point 2 and add the line between them
-    LOG(2, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << ".");
+    LOG(3, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << ".");
 
     //LOG(1, "INFO: OldSphereCenter is at " << helper << ".");
@@ -1244 +1257 @@
       for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++)
         output << *(*it);
-      LOG(2, "DEBUG: List of third Points is: " << output.str());
+      LOG(3, "DEBUG: List of third Points is: " << output.str());
     }
     if (!OptCandidates.pointlist.empty()) {
@@ -1406 +1419 @@
 //  return result;
 //};
-
 /** This function finds a triangle to a line, adjacent to an existing one.
  * @param out output stream for debugging
@@ -1438 +1450 @@
 
   // 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
@@ -1453 +1463 @@
     CircleRadius = RADIUS * RADIUS - radius / 4.;
     CirclePlaneNormal.Normalize();
-    LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
-
-    LOG(3, "DEBUG: OldSphereCenter is at " << T.SphereCenter << ".");
+    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+
+    LOG(4, "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(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
+    LOG(4, "DEBUG: SearchDirection is " << SearchDirection << ".");
     if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
@@ -1472 +1482 @@
 
   } else {
-    LOG(3, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!");
+    LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!");
   }
 
@@ -1505 +1515 @@
     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);
@@ -1538 +1548 @@
   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 << ":");
+    LOG(3, "DEBUG: List of Candidates for Turning Point " << *TurningPoint << ":" << output.str());
   }
 
@@ -1583 +1593 @@
   }
   delete (connectedClosestPoints);
-};
+}
+;
 
 /** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate.
@@ -1603 +1614 @@
       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.;
@@ -1610 +1621 @@
     }
   }
-};
+}
+;
 
 /** If a given \a *triangle is degenerated, this adds both sides.
@@ -1646 +1658 @@
   // give some verbose output about the whole procedure
   if (CandidateLine.T != NULL)
-    LOG(2, "DEBUG: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << ".");
+    LOG(2, "INFO: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << ".");
   else
-    LOG(2, "DEBUG: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle.");
+    LOG(2, "INFO: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle.");
   triangle = BTS;
 
@@ -1732 +1744 @@
 ;
 
+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
@@ -1868 +1894 @@
     }
     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);
     }
@@ -1897 +1923 @@
   class BoundaryLineSet *OldLines[4], *NewLine;
   class BoundaryPointSet *OldPoints[2];
-  Vector BaseLineNormal;
+  Vector BaseLineNormal[2];
+  Vector OtherEndpoint[2]; // fourth point to either triangle
   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;
   }
-  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()
+  {
+    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.");
+    }
+  }
 
   // get the two triangles
@@ -1941 +1972 @@
 
   // 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
@@ -1950 +1983 @@
         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) {
@@ -1965 +2006 @@
       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
@@ -1971 +2017 @@
   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;
@@ -1988 +2034 @@
   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
@@ -2000 +2052 @@
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
-    BTS->GetNormalVector(BaseLineNormal);
+    {
+      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]);
+    }
     AddTesselationTriangle(OldTriangleNrs[0]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");
@@ -2008 +2079 @@
     BLS[2] = NewLine;
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
-    BTS->GetNormalVector(BaseLineNormal);
+    {
+      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]);
+    }
     AddTesselationTriangle(OldTriangleNrs[1]);
     LOG(3, "DEBUG: Created new triangle " << *BTS << ".");
@@ -2155 +2245 @@
   TesselPoint *Candidate = NULL;
 
-  LOG(3, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << ".");
+  LOG(4, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << ".");
 
   // copy old center
@@ -2163 +2253 @@
 
   // 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;
@@ -2177 +2265 @@
     CircleRadius = RADIUS * RADIUS - radius;
     CirclePlaneNormal.Normalize();
-    LOG(3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
+    LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");
 
     // test whether old center is on the band's plane
@@ -2186 +2274 @@
     radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      LOG(3, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << ".");
+      LOG(4, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << ".");
 
       // check SearchDirection
-      LOG(3, "DEBUG: SearchDirection is " << SearchDirection << ".");
+      LOG(4, "DEBUG: SearchDirection is " << SearchDirection << ".");
       if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
         ELOG(1, "SearchDirection and RelativeOldSphereCenter are not orthogonal!");
@@ -2230 +2318 @@
                   // find center on the plane
                   GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
-                  LOG(3, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << ".");
+                  LOG(5, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << ".");
 
                   try {
-                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()),
-                                            (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
-                                            (Candidate->getPosition())).getNormal();
-                    LOG(3, "DEBUG: NewNormalVector is " << NewNormalVector << ".");
+                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()), (CandidateLine.BaseLine->endpoints[1]->node->getPosition()), (Candidate->getPosition())).getNormal();
+                    LOG(5, "DEBUG: NewNormalVector is " << NewNormalVector << ".");
                     radius = CandidateLine.BaseLine->endpoints[0]->node->DistanceSquared(NewPlaneCenter);
-                    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 << ".");
+                    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 << ".");
                     if (radius < RADIUS * RADIUS) {
                       otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter);
@@ -2246 +2332 @@
                         // construct both new centers
                         NewSphereCenter = NewPlaneCenter;
-                        OtherNewSphereCenter= NewPlaneCenter;
+                        OtherNewSphereCenter = NewPlaneCenter;
                         helper = NewNormalVector;
                         helper.Scale(sqrt(RADIUS * RADIUS - radius));
-                        LOG(4, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << ".");
+                        LOG(5, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << ".");
                         NewSphereCenter += helper;
-                        LOG(4, "DEBUG: NewSphereCenter is at " << NewSphereCenter << ".");
+                        LOG(5, "DEBUG: NewSphereCenter is at " << NewSphereCenter << ".");
                         // OtherNewSphereCenter is created by the same vector just in the other direction
                         helper.Scale(-1.);
                         OtherNewSphereCenter += helper;
-                        LOG(4, "DEBUG: OtherNewSphereCenter is at " << OtherNewSphereCenter << ".");
+                        LOG(5, "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);
@@ -2276 +2362 @@
                           if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
                             CandidateLine.pointlist.push_back(Candidate);
-                            LOG(2, "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
+                            LOG(3, "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(2, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
+                            LOG(3, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");
                           }
                           CandidateLine.ShortestAngle = alpha;
-                          LOG(2, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now.");
+                          LOG(4, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now.");
                         } else {
                           if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
-                            LOG(3, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " .");
+                            LOG(4, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " .");
                           } else {
-                            LOG(3, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected.");
+                            LOG(4, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected.");
                           }
                         }
                       } else {
-                        ELOG(0, "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius));
+                        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)));
                       }
                     } else {
-                      LOG(3, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << ".");
+                      LOG(4, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << ".");
                     }
-                  }
-                  catch (LinearDependenceException &excp){
-                    LOG(3, boost::diagnostic_information(excp));
-                    LOG(3, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent.");
+                  } catch (LinearDependenceException &excp) {
+                    LOG(4, boost::diagnostic_information(excp));
+                    LOG(4, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent.");
                   }
                 } else {
                   if (ThirdPoint != NULL) {
-                    LOG(3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << ".");
+                    LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << ".");
                   } else {
-                    LOG(3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << ".");
+                    LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << ".");
                   }
                 }
@@ -2318 +2403 @@
   } else {
     if (ThirdPoint != NULL)
-      LOG(3, "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!");
+      LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!");
     else
-      LOG(3, "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!");
-  }
-
-  LOG(2, "DEBUG: Sorting candidate list ...");
+      LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!");
+  }
+
+  LOG(5, "DEBUG: Sorting candidate list ...");
   if (CandidateLine.pointlist.size() > 1) {
     CandidateLine.pointlist.unique();
@@ -2329 +2414 @@
   }
 
-  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
-    ELOG(0, "There were other points contained in the rolling sphere as well!");
-    performCriticalExit();
-  }
+  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!"));
 }
 ;
@@ -2344 +2426 @@
 {
   //Info FunctionInfo(__func__);
-  const BoundaryLineSet * lines[2] = { line1, line2 };
+  const BoundaryLineSet * lines[2] = {line1, line2};
   class BoundaryPointSet *node = NULL;
   PointMap OrderMap;
@@ -2351 +2433 @@
     // 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;
@@ -2400 +2482 @@
               // 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.");
@@ -2446 +2527 @@
     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();
@@ -2507 +2586 @@
     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();
 
@@ -2524 +2602 @@
           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
@@ -2655 +2734 @@
   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.;
@@ -2764 +2843 @@
       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) {
@@ -2807 +2887 @@
   Vector OrthogonalVector;
   Vector helper;
-  const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL };
+  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
   TriangleList *triangles = NULL;
 
@@ -2818 +2898 @@
   // calculate central point
   triangles = FindTriangles(TrianglePoints);
-  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();
-  }
+  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;
   PlaneNormal.Scale(1.0 / triangles->size());
   LOG(4, "DEBUG: Calculated PlaneNormal of all circle points is " << PlaneNormal << ".");
@@ -2836 +2912 @@
     AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
-    if (AngleZero.NormSquared() < MYEPSILON) {
-      ELOG(0, "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!");
-      performCriticalExit();
-    }
+    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!"));
   }
   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);
@@ -2854 +2927 @@
     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)));
   }
 
@@ -2909 +2982 @@
   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++;
@@ -2920 +2991 @@
   }
   //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 << ".");
   //
@@ -2934 +3005 @@
     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);
-    if (AngleZero.NormSquared() < MYEPSILON) {
-      ELOG(0, "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!");
-      performCriticalExit();
-    }
+    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!"));
   }
   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);
@@ -2959 +3027 @@
       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)));
-    if (!InserterTest.second) {
-      ELOG(0, "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner));
-      performCriticalExit();
-    }
+    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))));
   }
 
@@ -2983 +3048 @@
   //Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPointList *> *ListOfPaths = new list<TesselPointList *> ;
+  list<TesselPointList *> *ListOfPaths = new list<TesselPointList *>;
   TesselPointList *connectedPath = NULL;
   Vector center;
@@ -3009 +3074 @@
   map<class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
   for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
-    TouchedLine.insert(pair<class BoundaryLineSet *, bool> (Runner->second, false));
-    for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++)
-      TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool> (Sprinter->second, false));
+    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));
+    }
   }
   if (!ReferencePoint->lines.empty()) {
@@ -3018 +3086 @@
       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(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
+        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 {
-                  LOG(1, "INFO: Skipping " << *triangle << ", as we have already visited it.");
+                  ELOG(1, "I could not find " << *triangle << " in the touched list.");
                   triangle = NULL;
                 }
               } else {
-                ELOG(1, "I could not find " << *triangle << " in the touched list.");
+                // as we have stumbled upon the same triangle, we don't need the check anymore
                 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(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);
           }
-          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;
-            }
-          }
-          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.");
-      }
+
+          ListOfPaths->push_back(connectedPath);
+        } else {
+          LOG(3, "DEBUG: Skipping " << *runner->second << ", as we have already visited it.");
+        }
     }
   } else {
@@ -3098 +3171 @@
   //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;
@@ -3108 +3181 @@
     connectedPath = *ListRunner;
 
-    LOG(1, "INFO: Current path is " << connectedPath << ".");
+    if (DoLog(2)) {
+      std::stringstream output;
+      output << "INFO: Current path is ";
+      BOOST_FOREACH(const TesselPoint * const item, *connectedPath) {
+        output << *item << " ";
+      }
+      LOG(1, output.str());
+    }
 
     // go through list, look for reappearance of starting Point and count
@@ -3117 +3197 @@
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
-        if (DoLog(1)) {
+        if (DoLog(3)) {
           std::stringstream output;
-          output << count + 1 << ". closed path consists of: ";
-          for (TesselPointList::iterator CircleSprinter = Marker;
-              CircleSprinter != CircleRunner;
-              CircleSprinter++)
+          output << "DEBUG: " << count + 1 << ". closed path consists of: ";
+          for (TesselPointList::iterator CircleSprinter = Marker; CircleSprinter != CircleRunner; CircleSprinter++)
             output << (**CircleSprinter) << " <-> ";
           LOG(1, output.str());
@@ -3138 +3216 @@
     }
   }
-  LOG(1, "INFO: " << count << " closed additional path(s) have been created.");
+  LOG(2, "DEBUG: " << count << " closed additional path(s) have been created.");
 
   // delete list of paths
@@ -3176 +3254 @@
 }
 ;
+
+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.
@@ -3190 +3519 @@
 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) {
@@ -3202 +3527 @@
     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()) {
@@ -3212 +3534 @@
 
   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();
+  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);
   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;
-  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.");
+  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);
   }
   delete (ListOfClosedPaths);
 
-  LOG(1, "INFO: Removed volume is " << volume << ".");
+  if (volume > 0.)
+    LOG(1, "INFO: Removed volume is " << volume << ".");
 
   return volume;
@@ -3437 +3724 @@
             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)) {
@@ -3458 +3745 @@
           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)) {
@@ -3487 +3774 @@
     }
     default:
-      ELOG(0, "Number of wildcards is greater than 3, cannot happen!");
-      performCriticalExit();
+      ASSERT(0, "Tesselation::FindTriangles() - Number of wildcards is greater than 3, cannot happen!");
       break;
   }
@@ -3514 +3800 @@
     if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])
       return true;
-    else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
-      return false;
-    else { // both lower-numbered endpoints are the same ...
-      if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2])
-        return true;
-      else if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2])
+    else
+      if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
         return false;
-    }
+      else { // both lower-numbered endpoints are the same ...
+        if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2])
+          return true;
+        else
+          if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2])
+            return false;
+      }
     return false;
   }
@@ -3542 +3830 @@
 
   // 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);
@@ -3564 +3852 @@
     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!");
   }
 
@@ -3598 +3886 @@
       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));
         }
       }
@@ -3626 +3914 @@
 
   // 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
@@ -3667 +3957 @@
 //                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++)
@@ -3683 +3973 @@
 
       // 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);
     }
   }
@@ -3734 +4027 @@
   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;
@@ -3789 +4080 @@
   for (int i = 0; i < 3; i++) { // look for the same line as BestLine (only it's its degenerated companion)
     if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
-      if (BestLine == BTS->lines[i]) {
-        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));
+      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));
       TempTriangle->lines[nr] = BTS->lines[i];
       break;
@@ -3815 +4103 @@
   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) {
@@ -3857 +4145 @@
     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++;
+    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;
       }
-      return false;
-    }
   }
 };
@@ -3896 +4186 @@
       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.");
@@ -3909 +4199 @@
         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));
@@ -3991 +4281 @@
     // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
     // connections to either polygon ...
-    if (T->size() % 2 != 0) {
-      ELOG(0, " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!");
-      performCriticalExit();
-    }
+    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!"));
     TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator
     /// 4a. Get NormalVector for one side (this is "front")