Changes in src/tesselation.cpp [29812d:7dea7c]

File:

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

src/tesselation.cpp

@@ -6 +6 @@
  */
 
+#include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "tesselation.hpp"
-#include "memoryallocator.hpp"
+#include "tesselationhelpers.hpp"
+#include "vector.hpp"
+
+class molecule;
 
 // ======================================== Points on Boundary =================================
@@ -37 +42 @@
 BoundaryPointSet::~BoundaryPointSet()
 {
-  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  //cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
   if (!lines.empty())
     cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
@@ -67 +72 @@
 ostream & operator <<(ostream &ost, BoundaryPointSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.node->Name << " at " << *a.node->node << "]";
   return ost;
 }
@@ -125 +130 @@
         for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
           if ((*Runner).second == this) {
-            cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+            //cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
             endpoints[i]->lines.erase(Runner);
             break;
           }
       } else { // there's just a single line left
-        if (endpoints[i]->lines.erase(Nr))
-          cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+        if (endpoints[i]->lines.erase(Nr)) {
+          //cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+        }
       }
       if (endpoints[i]->lines.empty()) {
-        cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+        //cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         if (endpoints[i] != NULL) {
           delete(endpoints[i]);
@@ -168 +174 @@
 
 /** Checks whether the adjacent triangles of a baseline are convex or not.
- * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
+ * We sum the two angles of each height vector with respect to the center of the baseline.
  * If greater/equal M_PI than we are convex.
  * \param *out output stream for debugging
@@ -178 +184 @@
   // get the two triangles
   if (triangles.size() != 2) {
-    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
+    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;
     return true;
   }
@@ -201 +207 @@
     NormalCheck.Scale(sign);
     sign = -sign;
-    BaseLineNormal.SubtractVector(&runner->second->NormalVector);   // we subtract as BaseLineNormal has to point inward in direction of [pi,2pi]
+    if (runner->second->NormalVector.NormSquared() > MYEPSILON)
+      BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
+    else {
+      *out << Verbose(1) << "CRITICAL: Triangle " << *runner->second << " has zero normal vector!" << endl;
+      exit(255);
+    }
     node = runner->second->GetThirdEndpoint(this);
     if (node != NULL) {
@@ -211 +222 @@
       i++;
     } else {
-      //*out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
+      //*out << Verbose(2) << "ERROR: I cannot find third node in triangle, something's wrong." << endl;
       return true;
     }
@@ -217 +228 @@
   //*out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
   if (NormalCheck.NormSquared() < MYEPSILON) {
-    *out << Verbose(2) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
+    *out << Verbose(3) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
     return true;
   }
+  BaseLineNormal.Scale(-1.);
   double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
   if ((angle - M_PI) > -MYEPSILON) {
-    *out << Verbose(2) << "ACCEPT: Angle is greater than pi: convex." << endl;
+    *out << Verbose(3) << "ACCEPT: Angle is greater than pi: convex." << endl;
     return true;
   } else {
-    *out << Verbose(2) << "REJECT: Angle is less than pi: concave." << endl;
+    *out << Verbose(3) << "REJECT: Angle is less than pi: concave." << endl;
     return false;
   }
@@ -262 +274 @@
 ostream & operator <<(ostream &ost, BoundaryLineSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]";
   return ost;
 };
@@ -332 +344 @@
   for (int i = 0; i < 3; i++) {
     if (lines[i] != NULL) {
-      if (lines[i]->triangles.erase(Nr))
-        cout << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      if (lines[i]->triangles.erase(Nr)) {
+        //cout << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      }
       if (lines[i]->triangles.empty()) {
-          cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+          //cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
           delete (lines[i]);
           lines[i] = NULL;
@@ -341 +354 @@
     }
   }
-  cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+  //cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
 };
 
@@ -361 +374 @@
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
  * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
- * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
- * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
- * smaller than the first line.
+ * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
+ * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
+ * the first two basepoints) or not.
  * \param *out output stream for debugging
  * \param *MolCenter offset vector of line
@@ -395 +408 @@
     exit(255);
   }
-  CrossPoint.SubtractVector(endpoints[i%3]->node->node);
+  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
 
   // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
     return true;
   } else { // outside!
@@ -426 +439 @@
   for(int i=0;i<3;i++)
     if (point == endpoints[i])
+      return true;
+  return false;
+};
+
+/** Checks whether point is any of the three endpoints this triangle contains.
+ * \param *point TesselPoint to test
+ * \return true - point is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
+{
+  for(int i=0;i<3;i++)
+    if (point == endpoints[i]->node)
       return true;
   return false;
@@ -451 +476 @@
 };
 
+/** Checks whether three given \a *Points coincide with triangle's endpoints.
+ * \param *Points[3] pointer to BoundaryPointSet
+ * \return true - is the very triangle, false - is not
+ */
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
+{
+  return (((endpoints[0] == T->endpoints[0])
+            || (endpoints[0] == T->endpoints[1])
+            || (endpoints[0] == T->endpoints[2])
+          ) && (
+            (endpoints[1] == T->endpoints[0])
+            || (endpoints[1] == T->endpoints[1])
+            || (endpoints[1] == T->endpoints[2])
+          ) && (
+            (endpoints[2] == T->endpoints[0])
+            || (endpoints[2] == T->endpoints[1])
+            || (endpoints[2] == T->endpoints[2])
+
+          ));
+};
+
 /** Returns the endpoint which is not contained in the given \a *line.
  * \param *line baseline defining two endpoints
@@ -485 +531 @@
 ostream &operator <<(ostream &ost, BoundaryTriangleSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
-      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
+      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
   return ost;
 };
@@ -505 +551 @@
 TesselPoint::~TesselPoint()
 {
-  Free(&Name);
 };
 
@@ -512 +557 @@
 ostream & operator << (ostream &ost, const TesselPoint &a)
 {
-  ost << "[" << (a.Name) << "|" << &a << "]";
+  ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
   return ost;
 };
@@ -569 +614 @@
   TrianglesOnBoundaryCount = 0;
   InternalPointer = PointsOnBoundary.begin();
+  LastTriangle = NULL;
+  TriangleFilesWritten = 0;
 }
 ;
@@ -585 +632 @@
       cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
   }
+  cout << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
 }
 ;
@@ -1052 +1100 @@
   Vector *Center = cloud->GetCenter(out);
   list<BoundaryTriangleSet*> *triangles = NULL;
+  bool AddFlag = false;
+  LinkedCell *BoundaryPoints = NULL;
 
   *out << Verbose(1) << "Begin of InsertStraddlingPoints" << endl;
 
   cloud->GoToFirst();
+  BoundaryPoints = new LinkedCell(this, 5.);
   while (!cloud->IsEnd()) {  // we only have to go once through all points, as boundary can become only bigger
-    LinkedCell BoundaryPoints(this, 5.);
+    if (AddFlag) {
+      delete(BoundaryPoints);
+      BoundaryPoints = new LinkedCell(this, 5.);
+      AddFlag = false;
+    }
     Walker = cloud->GetPoint();
     *out << Verbose(2) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToPoint(out, Walker->node, &BoundaryPoints);
-    if (triangles == NULL) {
-      *out << Verbose(1) << "No triangles found, probably a tesselation point itself." << endl;
+    triangles = FindClosestTrianglesToPoint(out, Walker->node, BoundaryPoints);
+    BTS = triangles->front();
+    if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
+      *out << Verbose(2) << "No triangles found, probably a tesselation point itself." << endl;
       cloud->GoToNext();
       continue;
     } else {
-      BTS = triangles->front();
     }
     *out << Verbose(2) << "Closest triangle is " << *BTS << "." << endl;
@@ -1090 +1145 @@
         // add Walker to boundary points
         *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
+        AddFlag = true;
         if (AddBoundaryPoint(Walker,0))
           NewPoint = BPS[0];
@@ -1180 +1236 @@
   } else {
     delete TPS[n];
-    cout << Verbose(3) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
+    cout << Verbose(4) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
     TPS[n] = (InsertUnique.first)->second;
   }
@@ -1197 +1253 @@
 
   if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine;
+    LineMap::iterator FindLine = a->lines.find(b->node->nr);
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-
-    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+    cout << Verbose(5) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
+
+    for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
       // If there is a line with less than two attached triangles, we don't need a new line.
       if (FindLine->second->triangles.size() < 2) {
         insertNewLine = false;
-        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
+        cout << Verbose(4) << "Using existing line " << *FindLine->second << endl;
 
         BPS[0] = FindLine->second->endpoints[0];
@@ -1232 +1289 @@
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
 {
-  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
+  cout << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
@@ -1242 +1299 @@
 };
 
-/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
- * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
+/** Function adds triangle to global list.
+ * Furthermore, the triangle receives the next free id and id counter \a TrianglesOnBoundaryCount is increased.
  */
 void Tesselation::AddTesselationTriangle()
@@ -1252 +1309 @@
   TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
   TrianglesOnBoundaryCount++;
+
+  // set as last new triangle
+  LastTriangle = BTS;
+
+  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
+};
+
+/** Function adds triangle to global list.
+ * Furthermore, the triangle number is set to \a nr.
+ * \param nr triangle number
+ */
+void Tesselation::AddTesselationTriangle(int nr)
+{
+  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
+
+  // add triangle to global map
+  TrianglesOnBoundary.insert(TrianglePair(nr, BTS));
+
+  // set as last new triangle
+  LastTriangle = BTS;
 
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
@@ -1272 +1349 @@
           cout << Verbose(5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
           RemoveTesselationLine(triangle->lines[i]);
-          triangle->lines[i] = NULL;
-      } else
-        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle." << endl;
+      } else {
+        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle: ";
+        for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
+          cout << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t";
+        cout << endl;
+//        for (int j=0;j<2;j++) {
+//          cout << Verbose(5) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
+//          for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
+//            cout << "[" << *(LineRunner->second) << "] \t";
+//          cout << endl;
+//        }
+      }
+      triangle->lines[i] = NULL;  // free'd or not: disconnect
     } else
       cerr << "ERROR: This line " << i << " has already been free'd." << endl;
@@ -1294 +1381 @@
   if (line == NULL)
     return;
-  // get other endpoint number of finding copies of same line
+  // get other endpoint number for finding copies of same line
   if (line->endpoints[1] != NULL)
     Numbers[0] = line->endpoints[1]->Nr;
@@ -1321 +1408 @@
         cout << Verbose(5) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         RemoveTesselationPoint(line->endpoints[i]);
-        line->endpoints[i] = NULL;
-      } else
-        cout << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to." << endl;
+      } else {
+        cout << Verbose(5) << *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++)
+          cout << "[" << *(LineRunner->second) << "] \t";
+        cout << endl;
+      }
+      line->endpoints[i] = NULL;  // free'd or not: disconnect
     } else
       cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
@@ -1390 +1481 @@
           }
           // Only one of the triangle lines must be considered for the triangle count.
-          *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-          return adjacentTriangleCount;
+          //*out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          //return adjacentTriangleCount;
         }
       }
@@ -1400 +1491 @@
   *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
   return adjacentTriangleCount;
+};
+
+/** Checks whether the triangle consisting of the three points is already present.
+ * Searches for the points in Tesselation::PointsOnBoundary and checks their
+ * lines. If any of the three edges already has two triangles attached, false is
+ * returned.
+ * \param *out output stream for debugging
+ * \param *Candidates endpoints of the triangle candidate
+ * \return NULL - none found or pointer to triangle
+ */
+class BoundaryTriangleSet * Tesselation::GetPresentTriangle(ofstream *out, TesselPoint *Candidates[3])
+{
+  class BoundaryTriangleSet *triangle = NULL;
+  class BoundaryPointSet *Points[3];
+
+  // builds a triangle point set (Points) of the end points
+  for (int i = 0; i < 3; i++) {
+    PointMap::iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      Points[i] = FindPoint->second;
+    } else {
+      Points[i] = NULL;
+    }
+  }
+
+  // checks lines between the points in the Points for their adjacent triangles
+  for (int i = 0; i < 3; i++) {
+    if (Points[i] != NULL) {
+      for (int j = i; j < 3; j++) {
+        if (Points[j] != NULL) {
+          LineMap::iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
+          for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
+            TriangleMap *triangles = &FindLine->second->triangles;
+            for (TriangleMap::iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
+              if (FindTriangle->second->IsPresentTupel(Points)) {
+                if ((triangle == NULL) || (triangle->Nr > FindTriangle->second->Nr))
+                  triangle = FindTriangle->second;
+              }
+            }
+          }
+          // Only one of the triangle lines must be considered for the triangle count.
+          //*out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          //return adjacentTriangleCount;
+        }
+      }
+    }
+  }
+
+  return triangle;
 };
 
@@ -1461 +1601 @@
   CandidateList *OptCandidates = new CandidateList();
   for (int k=0;k<NDIM;k++) {
+    Oben.Zero();
     Oben.x[k] = 1.;
     FirstPoint = MaxPoint[k];
@@ -1469 +1610 @@
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
 
-    FindSecondPointForTesselation(FirstPoint, NULL, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(FirstPoint, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     SecondPoint = OptCandidate;
     if (SecondPoint == NULL)  // have we found a second point?
       continue;
-    else
-      cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
 
     helper.CopyVector(FirstPoint->node);
@@ -1496 +1635 @@
 
     // adding point 1 and point 2 and add the line between them
+    cout << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;
     AddTesselationPoint(FirstPoint, 0);
+    cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
     AddTesselationPoint(SecondPoint, 1);
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -1569 +1710 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC)
 {
   cout << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
@@ -1604 +1745 @@
     CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
-    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
     // construct old center
@@ -1691 +1832 @@
         result = false;
       }
-    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
+    } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
         AddTesselationPoint((*it)->point, 0);
         AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
@@ -1732 +1873 @@
     // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
     BaseRay = BLS[0];
+    if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) {
+      *out << Verbose(1) << "CRITICAL: Triangle " << *BTS << " has zero normal vector!" << endl;
+      exit(255);
+    }
+
   }
 
@@ -1750 +1896 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return NULL - concave, otherwise endpoint that makes it concave
+ * \return NULL - convex, otherwise endpoint that makes it concave
  */
 class BoundaryPointSet *Tesselation::IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base)
@@ -1827 +1973 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return true - line was changed, false - same line as before
- */
-bool Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
+ * \return volume change due to flipping (0 - then no flipped occured)
+ */
+double Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
 {
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint[2];
+  double volume;
 
   int m=0;
@@ -1852 +1999 @@
   Distance.CopyVector(ClosestPoint[1]);
   Distance.SubtractVector(ClosestPoint[0]);
+
+  // calculate volume
+  volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->node, OtherBase->endpoints[0]->node->node, OtherBase->endpoints[1]->node->node, Base->endpoints[0]->node->node);
 
   // delete the temporary other base line and the closest points
@@ -1866 +2016 @@
     if (Base->triangles.size() < 2) {
       *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
-      return false;
+      return 0.;
     }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
@@ -1876 +2026 @@
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
       *out << Verbose(2) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;
-      FlipBaseline(out, Base);
-      return true;
+      // calculate volume summand as a general tetraeder
+      return volume;
     } else {  // Base higher than OtherBase -> do nothing
       *out << Verbose(2) << "REJECT: Base line is higher: Nothing to do." << endl;
-      return false;
-    }
-  }
-};
-
-/** Returns the closest point on \a *Base with respect to \a *OtherBase.
- * \param *out output stream for debugging
- * \param *Base reference line
- * \param *OtherBase other base line
- * \return Vector on reference line that has closest distance
- */
-Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
-{
-  // construct the plane of the two baselines (i.e. take both their directional vectors)
-  Vector Normal;
-  Vector Baseline, OtherBaseline;
-  Baseline.CopyVector(Base->endpoints[1]->node->node);
-  Baseline.SubtractVector(Base->endpoints[0]->node->node);
-  OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
-  OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
-  Normal.CopyVector(&Baseline);
-  Normal.VectorProduct(&OtherBaseline);
-  Normal.Normalize();
-  *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
-
-  // project one offset point of OtherBase onto this plane (and add plane offset vector)
-  Vector NewOffset;
-  NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
-  NewOffset.SubtractVector(Base->endpoints[0]->node->node);
-  NewOffset.ProjectOntoPlane(&Normal);
-  NewOffset.AddVector(Base->endpoints[0]->node->node);
-  Vector NewDirection;
-  NewDirection.CopyVector(&NewOffset);
-  NewDirection.AddVector(&OtherBaseline);
-
-  // calculate the intersection between this projected baseline and Base
-  Vector *Intersection = new Vector;
-  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
-  Normal.CopyVector(Intersection);
-  Normal.SubtractVector(Base->endpoints[0]->node->node);
-  *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
-
-  return Intersection;
+      return 0.;
+    }
+  }
 };
 
@@ -1930 +2040 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return true - flipping successful, false - something went awry
- */
-bool Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
+ * \return pointer to allocated new baseline - flipping successful, NULL - something went awry
+ */
+class BoundaryLineSet * Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
 {
   class BoundaryLineSet *OldLines[4], *NewLine;
@@ -1946 +2056 @@
   if (Base->triangles.size() < 2) {
     *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
-    return false;
+    return NULL;
   }
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
@@ -1982 +2092 @@
   if (i<4) {
     *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-    return false;
+    return NULL;
   }
   for (int j=0;j<4;j++)
     if (OldLines[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-      return false;
+      return NULL;
     }
   for (int j=0;j<2;j++)
     if (OldPoints[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
-      return false;
+      return NULL;
     }
 
@@ -2024 +2134 @@
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
     BTS->GetNormalVector(BaseLineNormal);
-    TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[0], BTS));
+    AddTesselationTriangle(OldTriangleNrs[0]);
     *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
 
@@ -2032 +2142 @@
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
     BTS->GetNormalVector(BaseLineNormal);
-    TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[1], BTS));
+    AddTesselationTriangle(OldTriangleNrs[1]);
     *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
   } else {
     *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
-    return false;
+    return NULL;
   }
 
   *out << Verbose(1) << "End of FlipBaseline" << endl;
-  return true;
+  return NewLine;
 };
 
@@ -2046 +2156 @@
 /** Finds the second point of starting triangle.
  * \param *a first node
- * \param *Candidate pointer to candidate node on return
  * \param Oben vector indicating the outside
  * \param OptCandidate reference to recommended candidate on return
@@ -2053 +2162 @@
  * \param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindSecondPointForTesselation(TesselPoint* a, TesselPoint* Candidate, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
+void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
 {
   cout << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl;
   Vector AngleCheck;
+  class TesselPoint* Candidate = NULL;
   double norm = -1., angle;
   LinkedNodes *List = NULL;
@@ -2191 +2301 @@
   cout << Verbose(1) << "Begin of FindThirdPointForTesselation" << endl;
 
-  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+  cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
 
   // construct center of circle
@@ -2216 +2326 @@
     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
+      //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 
       // check SearchDirection
@@ -2392 +2503 @@
   TesselPoint *trianglePoints[3];
   TesselPoint *SecondPoint = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
 
   if (LinesOnBoundary.empty()) {
@@ -2402 +2514 @@
   // check whether closest point is "too close" :), then it's inside
   if (trianglePoints[0] == NULL) {
-    *out << Verbose(1) << "Is the only point, no one else is closeby." << endl;
+    *out << Verbose(2) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
   if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
-    *out << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
-    return NULL;
-  }
-  list<TesselPoint*> *connectedPoints = GetCircleOfConnectedPoints(out, trianglePoints[0]);
-  list<TesselPoint*> *connectedClosestPoints = GetNeighboursOnCircleOfConnectedPoints(out, connectedPoints, trianglePoints[0], x);
-  delete(connectedPoints);
-  trianglePoints[1] = connectedClosestPoints->front();
-  trianglePoints[2] = connectedClosestPoints->back();
-  for (int i=0;i<3;i++) {
-    if (trianglePoints[i] == NULL) {
-      *out << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
-    }
-    //*out << Verbose(1) << "List of possible points:" << endl;
-    //*out << Verbose(2) << *trianglePoints[i] << endl;
-  }
-
-  list<BoundaryTriangleSet*> *triangles = FindTriangles(trianglePoints);
-
-  delete(connectedClosestPoints);
+    *out << Verbose(3) << "Point is right on a tesselation point, no nearest triangle." << endl;
+    PointMap::iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
+    triangles = new list<BoundaryTriangleSet*>;
+    if (PointRunner != PointsOnBoundary.end()) {
+      for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
+        for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
+          triangles->push_back(TriangleRunner->second);
+      triangles->sort();
+      triangles->unique();
+    } else {
+      PointRunner = PointsOnBoundary.find(SecondPoint->nr);
+      trianglePoints[0] = SecondPoint;
+      if (PointRunner != PointsOnBoundary.end()) {
+        for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
+          for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
+            triangles->push_back(TriangleRunner->second);
+        triangles->sort();
+        triangles->unique();
+      } else {
+        *out << Verbose(1) << "ERROR: I cannot find a boundary point to the tessel point " << *trianglePoints[0] << "." << endl;
+        return NULL;
+      }
+    }
+  } else {
+    list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(out, trianglePoints[0], x);
+    trianglePoints[1] = connectedClosestPoints->front();
+    trianglePoints[2] = connectedClosestPoints->back();
+    for (int i=0;i<3;i++) {
+      if (trianglePoints[i] == NULL) {
+        *out << Verbose(1) << "ERROR: IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+      }
+      //*out << Verbose(2) << "List of triangle points:" << endl;
+      //*out << Verbose(3) << *trianglePoints[i] << endl;
+    }
+
+    triangles = FindTriangles(trianglePoints);
+    *out << Verbose(2) << "List of possible triangles:" << endl;
+    for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      *out << Verbose(3) << **Runner << endl;
+
+    delete(connectedClosestPoints);
+  }
   
   if (triangles->empty()) {
-    *out << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
+    *out << Verbose(0) << "ERROR: There is no nearest triangle. Please check the tesselation structure.";
+    delete(triangles);
     return NULL;
   } else
@@ -2443 +2579 @@
   class BoundaryTriangleSet *result = NULL;
   list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(out, x, LC);
+  Vector Center;
 
   if (triangles == NULL)
     return NULL;
 
-  if (x->ScalarProduct(&triangles->front()->NormalVector) < 0)
-    result = triangles->back();
-  else
+  if (triangles->size() == 1) { // there is no degenerate case
     result = triangles->front();
-
+    *out << Verbose(2) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
+  } else {
+    result = triangles->front();
+    result->GetCenter(&Center);
+    Center.SubtractVector(x);
+    *out << Verbose(2) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
+    if (Center.ScalarProduct(&result->NormalVector) < 0) {
+      result = triangles->back();
+      *out << Verbose(2) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
+      if (Center.ScalarProduct(&result->NormalVector) < 0) {
+        *out << Verbose(1) << "ERROR: Front and back side yield NormalVector in wrong direction!" << endl;
+      }
+    }
+  }
   delete(triangles);
   return result;
@@ -2466 +2614 @@
 {
   class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, &Point, LC);
-  if (result == NULL)
+  Vector Center;
+
+  if (result == NULL) {// is boundary point or only point in point cloud?
+    *out << Verbose(1) << Point << " is the only point in vicinity." << endl;
+    return false;
+  }
+
+  result->GetCenter(&Center);
+  *out << Verbose(3) << "INFO: Central point of the triangle is " << Center << "." << endl;
+  Center.SubtractVector(&Point);
+  *out << Verbose(3) << "INFO: Vector from center to point to test is " << Center << "." << endl;
+  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
+    *out << Verbose(1) << Point << " is an inner point." << endl;
     return true;
-  if (Point.ScalarProduct(&result->NormalVector) < 0)
-    return true;
-  else
+  } else {
+    *out << Verbose(1) << Point << " is NOT an inner point." << endl;
     return false;
+  }
 }
 
@@ -2483 +2643 @@
 bool Tesselation::IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC)
 {
-  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, Point->node, LC);
-  if (result == NULL)
-    return true;
-  if (Point->node->ScalarProduct(&result->NormalVector) < 0)
-    return true;
-  else
-    return false;
+  return IsInnerPoint(out, *(Point->node), LC);
 }
 
@@ -2496 +2650 @@
  * @param *Point of which get all connected points
  *
- * @return list of the all points linked to the provided one
- */
-list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point)
-{
-  list<TesselPoint*> *connectedPoints = new list<TesselPoint*>;
+ * @return set of the all points linked to the provided one
+ */
+set<TesselPoint*> * Tesselation::GetAllConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -2510 +2664 @@
     ReferencePoint = PointRunner->second;
   } else {
-    *out << Verbose(2) << "getCircleOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    *out << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
     ReferencePoint = NULL;
   }
 
-  // little trick so that we look just through lines connect to the BoundaryPoint 
+  // little trick so that we look just through lines connect to the BoundaryPoint
   // OR fall-back to look through all lines if there is no such BoundaryPoint
   LineMap *Lines = &LinesOnBoundary;
@@ -2521 +2675 @@
   LineMap::iterator findLines = Lines->begin();
   while (findLines != Lines->end()) {
-    takePoint = false;
-
-    if (findLines->second->endpoints[0]->Nr == Point->nr) {
-      takePoint = true;
-      current = findLines->second->endpoints[1]->node;
-    } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
-      takePoint = true;
-      current = findLines->second->endpoints[0]->node;
-    }
-    
-    if (takePoint) {
-      *out << Verbose(3) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is taken into the circle." << endl;
-      connectedPoints->push_back(current);
-    }
-
-    findLines++;
+   takePoint = false;
+
+   if (findLines->second->endpoints[0]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[1]->node;
+   } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[0]->node;
+   }
+
+   if (takePoint) {
+     *out << Verbose(5) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
+     connectedPoints->insert(current);
+   }
+
+   findLines++;
   }
 
@@ -2543 +2697 @@
     return NULL;
   }
+
   return connectedPoints;
-}
-
-/** Gets the two neighbouring points with respect to a reference line to the provided point.
+};
+
+
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
  * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
  * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
@@ -2553 +2709 @@
  *
  * @param *out output stream for debugging
- * @param *connectedPoints list of connected points to the central \a *Point
  * @param *Point of which get all connected points
- * @param *Reference Vector to be checked whether it is an inner point
- *
- * @return list of the two points linked to the provided one and closest to the point to be checked,
- */
-list<TesselPoint*> * Tesselation::GetNeighboursOnCircleOfConnectedPoints(ofstream *out, list<TesselPoint*> *connectedPoints, TesselPoint* Point, Vector* Reference)
+ * @param *Reference Reference vector for zero angle or NULL for no preference
+ * @return list of the all points linked to the provided one
+ */
+list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector *Reference)
 {
   map<double, TesselPoint*> anglesOfPoints;
-  map<double, TesselPoint*>::iterator runner;
-  ;
-  Vector center, PlaneNormal, OrthogonalVector, helper, AngleZero;
-
-  if (connectedPoints->size() == 0) { // if have not found any points
-    *out << Verbose(1) << "ERROR: We have not found any connected points to " << *Point<< "." << endl;
-    return NULL;
-  }
+  set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(out, Point);
+  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
 
   // calculate central point
-  for (list<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
+  for (set<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //*out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -2586 +2739 @@
 
   // construct one orthogonal vector
-  AngleZero.CopyVector(Reference);
+  if (Reference != NULL)
+    AngleZero.CopyVector(Reference);
+  else
+    AngleZero.CopyVector((*connectedPoints->begin())->node);
   AngleZero.SubtractVector(Point->node);
   AngleZero.ProjectOntoPlane(&PlaneNormal);
@@ -2594 +2750 @@
 
   // go through all connected points and calculate angle
-  for (list<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
+  for (set<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
     helper.ProjectOntoPlane(&PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
-    *out << Verbose(2) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
+    *out << Verbose(3) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
     anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
   }
 
-  list<TesselPoint*> *result = new list<TesselPoint*>;
-  runner = anglesOfPoints.begin();
-  result->push_back(runner->second);
-  runner = anglesOfPoints.end();
-  runner--;
-  result->push_back(runner->second);
-
-  *out << Verbose(2) << "List of closest points has " << result->size() << " elements, which are "
-    << *(result->front()) << " and " << *(result->back()) << endl;
-
-  return result;
+  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
+    connectedCircle->push_back(AngleRunner->second);
+  }
+
+  delete(connectedPoints);
+  return connectedCircle;
 }
 
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we walk along a closed path.
+ *
+ * @param *out output stream for debugging
+ * @param *Point of which get all connected points
+ * @return list of the all points linked to the provided one
+ */
+list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  map<double, TesselPoint*> anglesOfPoints;
+  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+  class BoundaryPointSet *ReferencePoint = NULL;
+  class BoundaryPointSet *CurrentPoint = NULL;
+  class BoundaryTriangleSet *triangle = NULL;
+  class BoundaryLineSet *CurrentLine = NULL;
+  class BoundaryLineSet *StartLine = NULL;
+
+  // find the respective boundary point
+  PointMap::iterator PointRunner = PointsOnBoundary.find(Point->nr);
+  if (PointRunner != PointsOnBoundary.end()) {
+    ReferencePoint = PointRunner->second;
+  } else {
+    *out << Verbose(2) << "ERROR: GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    return NULL;
+  }
+
+  map <class BoundaryLineSet *, bool> TouchedLine;
+  map <class BoundaryTriangleSet *, bool> TouchedTriangle;
+  map <class BoundaryLineSet *, bool>::iterator LineRunner;
+  map <class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
+  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) );
+  }
+  if (!ReferencePoint->lines.empty()) {
+    for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++) {
+      LineRunner = TouchedLine.find(runner->second);
+      if (LineRunner == TouchedLine.end()) {
+        *out << Verbose(2) << "ERROR: I could not find " << *runner->second << " in the touched list." << endl;
+      } else if (!LineRunner->second) {
+        LineRunner->second = true;
+        connectedPath = new list<TesselPoint*>;
+        triangle = NULL;
+        CurrentLine = runner->second;
+        StartLine = CurrentLine;
+        CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+        *out << Verbose(3)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;
+        do {
+          // push current one
+          *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
+          connectedPath->push_back(CurrentPoint->node);
+
+          // find next triangle
+          for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
+            *out << Verbose(3) << "INFO: Inspecting triangle " << *Runner->second << "." << endl;
+            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;
+                  *out << Verbose(3) << "INFO: Connecting triangle is " << *triangle << "." << endl;
+                  break;
+                } else {
+                  *out << Verbose(3) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl;
+                  triangle = NULL;
+                }
+              } else {
+                *out << Verbose(2) << "ERROR: I could not find " << *triangle << " in the touched list." << endl;
+                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];
+              *out << Verbose(3) << "INFO: Connecting line is " << *CurrentLine << "." << endl;
+              break;
+            }
+          }
+          LineRunner = TouchedLine.find(CurrentLine);
+          if (LineRunner == TouchedLine.end())
+            *out << Verbose(2) << "ERROR: I could not find " << *CurrentLine << " in the touched list." << endl;
+          else
+            LineRunner->second = true;
+          // find next point
+          CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+
+        } while (CurrentLine != StartLine);
+        // last point is missing, as it's on start line
+        *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
+        if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
+          connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+
+        ListOfPaths->push_back(connectedPath);
+      } else {
+        *out << Verbose(3) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl;
+      }
+    }
+  } else {
+    *out << Verbose(1) << "ERROR: There are no lines attached to " << *ReferencePoint << "." << endl;
+  }
+
+  return ListOfPaths;
+}
+
+/** Gets all closed paths on the circle of points connected to the provided point by triangulation lines, if this very point is removed.
+ * From GetPathsOfConnectedPoints() extracts all single loops of intracrossing paths in the list of closed paths.
+ * @param *out output stream for debugging
+ * @param *Point of which get all connected points
+ * @return list of the closed paths
+ */
+list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(out, Point);
+  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPoint*> *newPath = NULL;
+  int count = 0;
+
+
+  list<TesselPoint*>::iterator CircleRunner;
+  list<TesselPoint*>::iterator CircleStart;
+
+  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+    connectedPath = *ListRunner;
+
+    *out << Verbose(2) << "INFO: Current path is " << connectedPath << "." << endl;
+
+    // go through list, look for reappearance of starting Point and count
+    CircleStart = connectedPath->begin();
+
+    // go through list, look for reappearance of starting Point and create list
+    list<TesselPoint*>::iterator Marker = CircleStart;
+    for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
+      if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
+        // we have a closed circle from Marker to new Marker
+        *out << Verbose(3) << count+1 << ". closed path consists of: ";
+        newPath = new list<TesselPoint*>;
+        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        for (; CircleSprinter != CircleRunner; CircleSprinter++) {
+          newPath->push_back(*CircleSprinter);
+          *out << (**CircleSprinter) << " <-> ";
+        }
+        *out << ".." << endl;
+        count++;
+        Marker = CircleRunner;
+
+        // add to list
+        ListofClosedPaths->push_back(newPath);
+      }
+    }
+  }
+  *out << Verbose(3) << "INFO: " << count << " closed additional path(s) have been created." << endl;
+
+  // delete list of paths
+  while (!ListofPaths->empty()) {
+    connectedPath = *(ListofPaths->begin());
+    ListofPaths->remove(connectedPath);
+    delete(connectedPath);
+  }
+  delete(ListofPaths);
+
+  // exit
+  return ListofClosedPaths;
+};
+
+
+/** Gets all belonging triangles for a given BoundaryPointSet.
+ * \param *out output stream for debugging
+ * \param *Point BoundaryPoint
+ * \return pointer to allocated list of triangles
+ */
+set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(ofstream *out, class BoundaryPointSet *Point)
+{
+  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+
+  if (Point == NULL) {
+    *out << Verbose(1) << "ERROR: Point given is NULL." << endl;
+  } else {
+    // go through its lines and insert all triangles
+    for (LineMap::iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
+      for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
+      connectedTriangles->insert(TriangleRunner->second);
+    }
+  }
+
+  return connectedTriangles;
+};
+
+
 /** Removes a boundary point from the envelope while keeping it closed.
- * We create new triangles and remove the old ones connected to the point.
+ * 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
@@ -2625 +2981 @@
   class BoundaryLineSet *line = NULL;
   class BoundaryTriangleSet *triangle = NULL;
-  Vector OldPoint, TetraederVector[3];
+  Vector OldPoint, NormalVector;
   double volume = 0;
-  int *numbers = NULL;
   int count = 0;
-  int i;
 
   if (point == NULL) {
@@ -2645 +2999 @@
     return 0.;
   }
-  list<TesselPoint*> *CircleofPoints = GetCircleOfConnectedPoints(out, point->node);
-
-  // remove all triangles
+
+  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(out, point->node);
+  list<TesselPoint*> *connectedPath = NULL;
+
+  // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  numbers = new int[count];
-  class BoundaryTriangleSet **Candidates = new BoundaryTriangleSet*[count];
-  i=0;
-  for (LineMap::iterator LineRunner = point->lines.begin(); (point != NULL) && (LineRunner != point->lines.end()); LineRunner++) {
+  map<class BoundaryTriangleSet *, int> Candidates;
+  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates[i] = triangle;
-      numbers[i++] = triangle->Nr;
-    }
-  }
-  for (int j=0;j<i;j++) {
-    RemoveTesselationTriangle(Candidates[j]);
-  }
-  delete[](Candidates);
-  *out << Verbose(1) << i << " triangles were removed." << endl;
-
-  // re-create all triangles by going through connected points list
-  list<TesselPoint*>::iterator CircleRunner = CircleofPoints->begin();
-  list<TesselPoint*>::iterator OtherCircleRunner = CircleofPoints->begin();
-  class TesselPoint *CentralNode = *CircleRunner;
-  // advance two with CircleRunner and one with OtherCircleRunner
-  CircleRunner++;
-  CircleRunner++;
-  OtherCircleRunner++;
-  i=0;
-  cout << Verbose(2) << "INFO: CentralNode is " << *CentralNode << "." << endl;
-  for (; (OtherCircleRunner != CircleofPoints->end()) && (CircleRunner != CircleofPoints->end()); (CircleRunner++), (OtherCircleRunner++)) {
-    cout << Verbose(3) << "INFO: CircleRunner's node is " << **CircleRunner << "." << endl;
-    cout << Verbose(3) << "INFO: OtherCircleRunner's node is " << **OtherCircleRunner << "." << endl;
-    *out << Verbose(4) << "Adding new triangle points."<< endl;
-    AddTesselationPoint(CentralNode, 0);
-    AddTesselationPoint(*OtherCircleRunner, 1);
-    AddTesselationPoint(*CircleRunner, 2);
-    *out << Verbose(4) << "Adding new triangle lines."<< endl;
-    AddTesselationLine(TPS[0], TPS[1], 0);
-    AddTesselationLine(TPS[0], TPS[2], 1);
-    AddTesselationLine(TPS[1], TPS[2], 2);
-    BTS = new class BoundaryTriangleSet(BLS, numbers[i]);
-    TrianglesOnBoundary.insert(TrianglePair(numbers[i], BTS));
-    *out << Verbose(4) << "Created triangle " << *BTS << "." << endl;
-    // calculate volume summand as a general tetraeder
-    for (int j=0;j<3;j++) {
-      TetraederVector[j].CopyVector(TPS[j]->node->node);
-      TetraederVector[j].SubtractVector(&OldPoint);
-    }
-    OldPoint.CopyVector(&TetraederVector[0]);
-    OldPoint.VectorProduct(&TetraederVector[1]);
-    volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
-    // advance number
-    i++;
-    if (i >= count)
-      *out << Verbose(2) << "WARNING: Maximum of numbers reached!" << endl;
-  }
-  *out << Verbose(1) << i << " triangles were created." << endl;
-
-  delete[](numbers);
-
-  return volume;
-};
-
-/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
- * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
- * make it bigger (i.e. closing one (the baseline) and opening two new ones).
- * \param TPS[3] nodes of the triangle
- * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
- */
-bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3])
-{
-  bool result = false;
-  int counter = 0;
-
-  // check all three points
-  for (int i=0;i<3;i++)
-    for (int j=i+1; j<3; j++) {
-      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
-        LineMap::iterator FindLine;
-        pair<LineMap::iterator,LineMap::iterator> FindPair;
-        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
-        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-          // If there is a line with less than two attached triangles, we don't need a new line.
-          if (FindLine->second->triangles.size() < 2) {
-            counter++;
-            break;  // increase counter only once per edge
-          }
+      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+    }
+  }
+
+  // remove all triangles
+  count=0;
+  NormalVector.Zero();
+  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    *out << Verbose(3) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
+    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->first);
+    count++;
+  }
+  *out << Verbose(1) << count << " triangles were removed." << endl;
+
+  list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
+  map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
+  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  double angle;
+  double smallestangle;
+  Vector Point, Reference, OrthogonalVector;
+  if (count > 2) {  // less than three triangles, then nothing will be created
+    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
+      list<class BoundaryLineSet *> NewLines;
+      for (;!connectedPath->empty();) {
+        // search middle node with widest angle to next neighbours
+        EndNode = connectedPath->end();
+        smallestangle = 0.;
+        for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
+          cout << Verbose(3) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+          // construct vectors to next and previous neighbour
+          StartNode = MiddleNode;
+          if (StartNode == connectedPath->begin())
+            StartNode = connectedPath->end();
+          StartNode--;
+          //cout << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
+          Point.CopyVector((*StartNode)->node);
+          Point.SubtractVector((*MiddleNode)->node);
+          StartNode = MiddleNode;
+          StartNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          //cout << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
+          Reference.CopyVector((*StartNode)->node);
+          Reference.SubtractVector((*MiddleNode)->node);
+          OrthogonalVector.CopyVector((*MiddleNode)->node);
+          OrthogonalVector.SubtractVector(&OldPoint);
+          OrthogonalVector.MakeNormalVector(&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;
+            }
         }
-      } else { // no line
-        cout << Verbose(1) << "The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
-        result = true;
+        MiddleNode = EndNode;
+        if (MiddleNode == connectedPath->end()) {
+          cout << Verbose(1) << "CRITICAL: Could not find a smallest angle!" << endl;
+          exit(255);
+        }
+        StartNode = MiddleNode;
+        if (StartNode == connectedPath->begin())
+          StartNode = connectedPath->end();
+        StartNode--;
+        EndNode++;
+        if (EndNode == connectedPath->end())
+          EndNode = connectedPath->begin();
+        cout << Verbose(4) << "INFO: StartNode is " << **StartNode << "." << endl;
+        cout << Verbose(4) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+        cout << Verbose(4) << "INFO: EndNode is " << **EndNode << "." << endl;
+        *out << Verbose(3) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl;
+        TriangleCandidates[0] = *StartNode;
+        TriangleCandidates[1] = *MiddleNode;
+        TriangleCandidates[2] = *EndNode;
+        triangle = GetPresentTriangle(out, TriangleCandidates);
+        if (triangle != NULL) {
+          cout << Verbose(1) << "WARNING: New triangle already present, skipping!" << endl;
+          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;
+        }
+        *out << Verbose(5) << "Adding new triangle points."<< endl;
+        AddTesselationPoint(*StartNode, 0);
+        AddTesselationPoint(*MiddleNode, 1);
+        AddTesselationPoint(*EndNode, 2);
+        *out << Verbose(5) << "Adding new triangle lines."<< endl;
+        AddTesselationLine(TPS[0], TPS[1], 0);
+        AddTesselationLine(TPS[0], TPS[2], 1);
+        NewLines.push_back(BLS[1]);
+        AddTesselationLine(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->node, TPS[1]->node->node, TPS[2]->node->node, &OldPoint);
+        // advance number
+        count++;
+
+        // prepare nodes for next triangle
+        StartNode = EndNode;
+        cout << Verbose(4) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl;
+        connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
+        if (connectedPath->size() == 2) { // we are done
+          connectedPath->remove(*StartNode); // remove the start node
+          connectedPath->remove(*EndNode); // remove the end node
+          break;
+        } else if (connectedPath->size() < 2) { // something's gone wrong!
+          cout << Verbose(1) << "CRITICAL: There are only two endpoints left!" << endl;
+          exit(255);
+        } else {
+          MiddleNode = StartNode;
+          MiddleNode++;
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          EndNode = MiddleNode;
+          EndNode++;
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+        }
       }
-    }
-  if (counter > 1) {
-    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
-    result = true;
-  }
-  return result;
-};
-
-
-/** Sort function for the candidate list.
- */
-bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
-{
-  Vector BaseLineVector, OrthogonalVector, helper;
-  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
-    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
-    //return false;
-    exit(1);
-  }
-  // create baseline vector
-  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
-  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
-  BaseLineVector.Normalize();
-
-  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
-  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
-  helper.SubtractVector(candidate1->point->node);
-  OrthogonalVector.CopyVector(&helper);
-  helper.VectorProduct(&BaseLineVector);
-  OrthogonalVector.SubtractVector(&helper);
-  OrthogonalVector.Normalize();
-
-  // calculate both angles and correct with in-plane vector
-  helper.CopyVector(candidate1->point->node);
-  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
-  double phi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    phi = 2.*M_PI - phi;
-  }
-  helper.CopyVector(candidate2->point->node);
-  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
-  double psi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    psi = 2.*M_PI - psi;
-  }
-
-  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
-  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
-
-  // return comparison
-  return phi < psi;
-};
-
-/**
- * Finds the point which is second closest to the provided one.
- *
- * @param Point to which to find the second closest other point
- * @param linked cell structure
- *
- * @return point which is second closest to the provided one
- */
-TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
-{
-  LinkedNodes *List = NULL;
-  TesselPoint* closestPoint = NULL;
-  TesselPoint* secondClosestPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //cout << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
-            if (currentNorm < distance) {
-              // remember second point
-              secondDistance = distance;
-              secondClosestPoint = closestPoint;
-              // mark down new closest point
-              distance = currentNorm;
-              closestPoint = (*Runner);
-              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
+      // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
+      if (NewLines.size() > 1) {
+        list<class BoundaryLineSet *>::iterator Candidate;
+        class BoundaryLineSet *OtherBase = NULL;
+        double tmp, maxgain;
+        do {
+          maxgain = 0;
+          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+            tmp = PickFarthestofTwoBaselines(out, *Runner);
+            if (maxgain < tmp) {
+              maxgain = tmp;
+              Candidate = Runner;
             }
           }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
-        }
+          if (maxgain != 0) {
+            volume += maxgain;
+            cout << Verbose(3) << "Flipping baseline with highest volume" << **Candidate << "." << endl;
+            OtherBase = FlipBaseline(out, *Candidate);
+            NewLines.erase(Candidate);
+            NewLines.push_back(OtherBase);
+          }
+        } while (maxgain != 0.);
       }
 
-  return secondClosestPoint;
-};
-
-/**
- * Finds the point which is closest to the provided one.
- *
- * @param Point to which to find the closest other point
- * @param SecondPoint the second closest other point on return, NULL if none found
- * @param linked cell structure
- *
- * @return point which is closest to the provided one, NULL if none found
- */
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
-{
-  LinkedNodes *List = NULL;
-  TesselPoint* closestPoint = NULL;
-  SecondPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //cout << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        cout << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
-            if (currentNorm < distance) {
-              secondDistance = distance;
-              SecondPoint = closestPoint;
-              distance = currentNorm;
-              closestPoint = (*Runner);
-              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
-            } else if (currentNorm < secondDistance) {
-              secondDistance = currentNorm;
-              SecondPoint = (*Runner);
-              cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
-            }
-          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-
-  return closestPoint;
-};
+      ListOfClosedPaths->remove(connectedPath);
+      delete(connectedPath);
+    }
+    *out << Verbose(1) << count << " triangles were created." << endl;
+  } else {
+    while (!ListOfClosedPaths->empty()) {
+      ListRunner = ListOfClosedPaths->begin();
+      connectedPath = *ListRunner;
+      ListOfClosedPaths->remove(connectedPath);
+      delete(connectedPath);
+    }
+    *out << Verbose(1) << "No need to create any triangles." << endl;
+  }
+  delete(ListOfClosedPaths);
+
+  *out << Verbose(1) << "Removed volume is " << volume << "." << endl;
+
+  return volume;
+};
+
+
 
 /**
@@ -2939 +3228 @@
               FindTriangle = FindLine->second->triangles.begin();
               for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if ((
-                  (FindTriangle->second->endpoints[0] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[1] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[2] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
-                  )
-                ) {
+                if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
                   result->push_back(FindTriangle->second);
                 }
@@ -2970 +3246 @@
 
 /**
+ * Finds all degenerated lines within the tesselation structure.
+ *
+ * @return map of keys of degenerated line pairs, each line occurs twice
+ *         in the list, once as key and once as value
+ */
+map<int, int> * Tesselation::FindAllDegeneratedLines()
+{
+  map<int, class BoundaryLineSet *> AllLines;
+  map<int, int> * DegeneratedLines = new map<int, int>;
+
+  // sanity check
+  if (LinesOnBoundary.empty()) {
+    cout << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
+    return DegeneratedLines;
+  }
+
+  LineMap::iterator LineRunner1;
+  pair<LineMap::iterator, bool> tester;
+  for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
+    tester = AllLines.insert( pair<int,BoundaryLineSet *> (LineRunner1->second->endpoints[0]->Nr, LineRunner1->second) );
+    if ((!tester.second) && (tester.first->second->endpoints[1]->Nr == LineRunner1->second->endpoints[1]->Nr)) { // found degenerated line
+      DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, tester.first->second->Nr) );
+      DegeneratedLines->insert ( pair<int, int> (tester.first->second->Nr, LineRunner1->second->Nr) );
+    }
+  }
+
+  AllLines.clear();
+
+  cout << Verbose(1) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
+  map<int,int>::iterator it;
+  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++)
+      cout << Verbose(2) << (*it).first << " => " << (*it).second << endl;
+
+  return DegeneratedLines;
+}
+
+/**
  * Finds all degenerated triangles within the tesselation structure.
  *
@@ -2975 +3288 @@
  *         in the list, once as key and once as value
  */
-map<int, int> Tesselation::FindAllDegeneratedTriangles()
-{
-  map<int, int> DegeneratedTriangles;
-
-  // sanity check
-  if (LinesOnBoundary.empty()) {
-    cout << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
-    return DegeneratedTriangles;
-  }
-
-  LineMap::iterator LineRunner1, LineRunner2;
-
-  for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
-    for (LineRunner2 = LinesOnBoundary.begin(); LineRunner2 != LinesOnBoundary.end(); ++LineRunner2) {
-      if ((LineRunner1->second != LineRunner2->second)
-        && (LineRunner1->second->endpoints[0] == LineRunner2->second->endpoints[0])
-        && (LineRunner1->second->endpoints[1] == LineRunner2->second->endpoints[1])
-      ) {
-        TriangleMap::iterator TriangleRunner1 = LineRunner1->second->triangles.begin(),
-          TriangleRunner2 = LineRunner2->second->triangles.begin();
-
-        for (; TriangleRunner1 != LineRunner1->second->triangles.end(); ++TriangleRunner1) {
-          for (; TriangleRunner2 != LineRunner2->second->triangles.end(); ++TriangleRunner2) {
-            if ((TriangleRunner1->second != TriangleRunner2->second)
-              && (TriangleRunner1->second->endpoints[0] == TriangleRunner2->second->endpoints[0])
-              && (TriangleRunner1->second->endpoints[1] == TriangleRunner2->second->endpoints[1])
-              && (TriangleRunner1->second->endpoints[2] == TriangleRunner2->second->endpoints[2])
-            ) {
-              DegeneratedTriangles[TriangleRunner1->second->Nr] = TriangleRunner2->second->Nr;
-              DegeneratedTriangles[TriangleRunner2->second->Nr] = TriangleRunner1->second->Nr;
-            }
-          }
+map<int, int> * Tesselation::FindAllDegeneratedTriangles()
+{
+  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
+  map<int, int> * DegeneratedTriangles = new map<int, int>;
+
+  TriangleMap::iterator TriangleRunner1, TriangleRunner2;
+  LineMap::iterator Liner;
+  class BoundaryLineSet *line1 = NULL, *line2 = NULL;
+
+  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+    // run over both lines' triangles
+    Liner = LinesOnBoundary.find(LineRunner->first);
+    if (Liner != LinesOnBoundary.end())
+      line1 = Liner->second;
+    Liner = LinesOnBoundary.find(LineRunner->second);
+    if (Liner != LinesOnBoundary.end())
+      line2 = Liner->second;
+    for (TriangleRunner1 = line1->triangles.begin(); TriangleRunner1 != line1->triangles.end(); ++TriangleRunner1) {
+      for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) {
+        if ((TriangleRunner1->second != TriangleRunner2->second)
+          && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr) );
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr) );
         }
       }
     }
   }
-
-  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles.size() << " triangles." << endl;
+  delete(DegeneratedLines);
+
+  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
   map<int,int>::iterator it;
-  for (it = DegeneratedTriangles.begin(); it != DegeneratedTriangles.end(); it++)
+  for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       cout << Verbose(2) << (*it).first << " => " << (*it).second << endl;
 
@@ -3026 +3331 @@
 void Tesselation::RemoveDegeneratedTriangles()
 {
-  map<int, int> DegeneratedTriangles = FindAllDegeneratedTriangles();
-
-  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles.begin();
-    TriangleKeyRunner != DegeneratedTriangles.end(); ++TriangleKeyRunner
+  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  TriangleMap::iterator finder;
+  BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
+  int count  = 0;
+
+  cout << Verbose(1) << "Begin of RemoveDegeneratedTriangles" << endl;
+
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+    TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
-    BoundaryTriangleSet *triangle = TrianglesOnBoundary.find(TriangleKeyRunner->first)->second,
-      *partnerTriangle = TrianglesOnBoundary.find(TriangleKeyRunner->second)->second;
+    finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
+    if (finder != TrianglesOnBoundary.end())
+      triangle = finder->second;
+    else
+      break;
+    finder = TrianglesOnBoundary.find(TriangleKeyRunner->second);
+    if (finder != TrianglesOnBoundary.end())
+      partnerTriangle = finder->second;
+    else
+      break;
 
     bool trianglesShareLine = false;
@@ -3044 +3362 @@
       && (triangle->endpoints[0]->LinesCount > 2)
     ) {
+      // check whether we have to fix lines
+      BoundaryTriangleSet *Othertriangle = NULL;
+      BoundaryTriangleSet *OtherpartnerTriangle = NULL;
+      TriangleMap::iterator TriangleRunner;
+      for (int i = 0; i < 3; ++i)
+        for (int j = 0; j < 3; ++j)
+          if (triangle->lines[i] != partnerTriangle->lines[j]) {
+            // get the other two triangles
+            for (TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); ++TriangleRunner)
+              if (TriangleRunner->second != triangle) {
+                Othertriangle = TriangleRunner->second;
+              }
+            for (TriangleRunner = partnerTriangle->lines[i]->triangles.begin(); TriangleRunner != partnerTriangle->lines[i]->triangles.end(); ++TriangleRunner)
+              if (TriangleRunner->second != partnerTriangle) {
+                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);
+            triangle->lines[i]->triangles.insert( TrianglePair( partnerTriangle->Nr, partnerTriangle) );
+            for (int k=0;k<3;k++)
+              if (triangle->lines[i] == Othertriangle->lines[k]) {
+                Othertriangle->lines[k] = partnerTriangle->lines[j];
+                break;
+              }
+            // the line of partnerTriangle receives the non-degenerated ones
+            partnerTriangle->lines[j]->triangles.erase( partnerTriangle->Nr);
+            partnerTriangle->lines[j]->triangles.insert( TrianglePair( Othertriangle->Nr, Othertriangle) );
+            partnerTriangle->lines[j] = triangle->lines[i];
+          }
+
+      // erase the pair
+      count += (int) DegeneratedTriangles->erase(triangle->Nr);
       cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;
       RemoveTesselationTriangle(triangle);
+      count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
       cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;
       RemoveTesselationTriangle(partnerTriangle);
-      DegeneratedTriangles.erase(DegeneratedTriangles.find(partnerTriangle->Nr));
     } else {
       cout << Verbose(1) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle
@@ -3055 +3406 @@
     }
   }
+  delete(DegeneratedTriangles);
+
+  cout << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl;
+  cout << Verbose(1) << "End of RemoveDegeneratedTriangles" << endl;
 }
 
-/** Gets the angle between a point and a reference relative to the provided center.
- * We have two shanks point and reference between which the angle is calculated
- * and by scalar product with OrthogonalVector we decide the interval.
- * @param point to calculate the angle for
- * @param reference to which to calculate the angle
- * @param OrthogonalVector points in direction of [pi,2pi] interval
- *
- * @return angle between point and reference
- */
-double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
-{
-  if (reference.IsZero())
-    return M_PI;
-
-  // calculate both angles and correct with in-plane vector
-  if (point.IsZero())
-    return M_PI;
-  double phi = point.Angle(&reference);
-  if (OrthogonalVector.ScalarProduct(&point) > 0) {
-    phi = 2.*M_PI - phi;
-  }
-
-  cout << Verbose(3) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
-
-  return phi;
-}
-
+/** Adds an outside Tesselpoint to the envelope via (two) degenerated triangles.
+ * We look for the closest point on the boundary, we look through its connected boundary lines and
+ * seek the one with the minimum angle between its center point and the new point and this base line.
+ * We open up the line by adding a degenerated triangle, whose other side closes the base line again.
+ * \param *out output stream for debugging
+ * \param *point point to add
+ * \param *LC Linked Cell structure to find nearest point
+ */
+void Tesselation::AddBoundaryPointByDegeneratedTriangle(ofstream *out, class TesselPoint *point, LinkedCell *LC)
+{
+  *out << Verbose(2) << "Begin of AddBoundaryPointByDegeneratedTriangle" << endl;
+
+  // find nearest boundary point
+  class TesselPoint *BackupPoint = NULL;
+  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
+  class BoundaryPointSet *NearestBoundaryPoint = NULL;
+  PointMap::iterator PointRunner;
+
+  if (NearestPoint == point)
+    NearestPoint = BackupPoint;
+  PointRunner = PointsOnBoundary.find(NearestPoint->nr);
+  if (PointRunner != PointsOnBoundary.end()) {
+    NearestBoundaryPoint = PointRunner->second;
+  } else {
+    *out << Verbose(1) << "ERROR: I cannot find the boundary point." << endl;
+    return;
+  }
+  *out << Verbose(2) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl;
+
+  // go through its lines and find the best one to split
+  Vector CenterToPoint;
+  Vector BaseLine;
+  double angle, BestAngle = 0.;
+  class BoundaryLineSet *BestLine = NULL;
+  for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
+    BaseLine.CopyVector(Runner->second->endpoints[0]->node->node);
+    BaseLine.SubtractVector(Runner->second->endpoints[1]->node->node);
+    CenterToPoint.CopyVector(Runner->second->endpoints[0]->node->node);
+    CenterToPoint.AddVector(Runner->second->endpoints[1]->node->node);
+    CenterToPoint.Scale(0.5);
+    CenterToPoint.SubtractVector(point->node);
+    angle = CenterToPoint.Angle(&BaseLine);
+    if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
+      BestAngle = angle;
+      BestLine = Runner->second;
+    }
+  }
+
+  // remove one triangle from the chosen line
+  class BoundaryTriangleSet *TempTriangle = (BestLine->triangles.begin())->second;
+  BestLine->triangles.erase(TempTriangle->Nr);
+  int nr = -1;
+  for (int i=0;i<3; i++) {
+    if (TempTriangle->lines[i] == BestLine) {
+      nr = i;
+      break;
+    }
+  }
+
+  // create new triangle to connect point (connects automatically with the missing spot of the chosen line)
+  *out << Verbose(5) << "Adding new triangle points."<< endl;
+  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
+  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
+  AddTesselationPoint(point, 2);
+  *out << Verbose(5) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  BTS->GetNormalVector(TempTriangle->NormalVector);
+  BTS->NormalVector.Scale(-1.);
+  *out << Verbose(3) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl;
+  AddTesselationTriangle();
+
+  // create other side of this triangle and close both new sides of the first created triangle
+  *out << Verbose(5) << "Adding new triangle points."<< endl;
+  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
+  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
+  AddTesselationPoint(point, 2);
+  *out << Verbose(5) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  BTS->GetNormalVector(TempTriangle->NormalVector);
+  *out << Verbose(3) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;
+  AddTesselationTriangle();
+
+  // add removed triangle to the last open line of the second triangle
+  for (int i=0;i<3;i++) { // look for the same line as BestLine (only it's its degenerated companion)
+    if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
+      if (BestLine == BTS->lines[i]){
+        *out << Verbose(1) << "CRITICAL: BestLine is same as found line, something's wrong here!" << endl;
+        exit(255);
+      }
+      BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) );
+      TempTriangle->lines[nr] = BTS->lines[i];
+      break;
+    }
+  }
+
+  // exit
+  *out << Verbose(2) << "End of AddBoundaryPointByDegeneratedTriangle" << endl;
+};
+
+/** Writes the envelope to file.
+ * \param *out otuput stream for debugging
+ * \param *filename basename of output file
+ * \param *cloud PointCloud structure with all nodes
+ */
+void Tesselation::Output(ofstream *out, const char *filename, PointCloud *cloud)
+{
+  ofstream *tempstream = NULL;
+  string NameofTempFile;
+  char NumberName[255];
+
+  if (LastTriangle != NULL) {
+    sprintf(NumberName, "-%04d-%s_%s_%s", (int)TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name);
+    if (DoTecplotOutput) {
+      string NameofTempFile(filename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+      NameofTempFile.erase(npos, 1);
+      NameofTempFile.append(TecplotSuffix);
+      *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      WriteTecplotFile(out, tempstream, this, cloud, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+
+    if (DoRaster3DOutput) {
+      string NameofTempFile(filename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+      NameofTempFile.erase(npos, 1);
+      NameofTempFile.append(Raster3DSuffix);
+      *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      WriteRaster3dFile(out, tempstream, this, cloud);
+      IncludeSphereinRaster3D(out, tempstream, this, cloud);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+  }
+  if (DoTecplotOutput || DoRaster3DOutput)
+    TriangleFilesWritten++;
+};