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

File:

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

src/tesselation.cpp

@@ -6 +6 @@
  */
 
-#include <fstream>
-
-#include "linkedcell.hpp"
 #include "tesselation.hpp"
-#include "tesselationhelpers.hpp"
-#include "vector.hpp"
-#include "verbose.hpp"
-
-class molecule;
+#include "memoryallocator.hpp"
 
 // ======================================== Points on Boundary =================================
@@ -44 +37 @@
 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;
@@ -74 +67 @@
 ostream & operator <<(ostream &ost, BoundaryPointSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.node->Name << " at " << *a.node->node << "]";
+  ost << "[" << a.Nr << "|" << a.node->Name << "]";
   return ost;
 }
@@ -132 +125 @@
         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]);
@@ -176 +168 @@
 
 /** Checks whether the adjacent triangles of a baseline are convex or not.
- * We sum the two angles of each height vector with respect to the center of the baseline.
+ * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
  * If greater/equal M_PI than we are convex.
  * \param *out output stream for debugging
@@ -186 +178 @@
   // get the two triangles
   if (triangles.size() != 2) {
-    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;
+    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
     return true;
   }
@@ -209 +201 @@
     NormalCheck.Scale(sign);
     sign = -sign;
-    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);
-    }
+    BaseLineNormal.SubtractVector(&runner->second->NormalVector);   // we subtract as BaseLineNormal has to point inward in direction of [pi,2pi]
     node = runner->second->GetThirdEndpoint(this);
     if (node != NULL) {
@@ -224 +211 @@
       i++;
     } else {
-      //*out << Verbose(2) << "ERROR: I cannot find third node in triangle, something's wrong." << endl;
+      //*out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
       return true;
     }
@@ -230 +217 @@
   //*out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
   if (NormalCheck.NormSquared() < MYEPSILON) {
-    *out << Verbose(3) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
+    *out << Verbose(2) << "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(3) << "ACCEPT: Angle is greater than pi: convex." << endl;
+    *out << Verbose(2) << "ACCEPT: Angle is greater than pi: convex." << endl;
     return true;
   } else {
-    *out << Verbose(3) << "REJECT: Angle is less than pi: concave." << endl;
+    *out << Verbose(2) << "REJECT: Angle is less than pi: concave." << endl;
     return false;
   }
@@ -276 +262 @@
 ostream & operator <<(ostream &ost, BoundaryLineSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
   return ost;
 };
@@ -346 +332 @@
   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;
@@ -356 +341 @@
     }
   }
-  //cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+  cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
 };
 
@@ -376 +361 @@
  * 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: 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.
+ * 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.
  * \param *out output stream for debugging
  * \param *MolCenter offset vector of line
@@ -410 +395 @@
     exit(255);
   }
-  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+  CrossPoint.SubtractVector(endpoints[i%3]->node->node);
 
   // 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!
@@ -441 +426 @@
   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;
@@ -478 +451 @@
 };
 
-/** 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
@@ -533 +485 @@
 ostream &operator <<(ostream &ost, BoundaryTriangleSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
-      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
   return ost;
 };
@@ -553 +505 @@
 TesselPoint::~TesselPoint()
 {
+  Free(&Name);
 };
 
@@ -559 +512 @@
 ostream & operator << (ostream &ost, const TesselPoint &a)
 {
-  ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
+  ost << "[" << (a.Name) << "|" << &a << "]";
   return ost;
 };
@@ -616 +569 @@
   TrianglesOnBoundaryCount = 0;
   InternalPointer = PointsOnBoundary.begin();
-  LastTriangle = NULL;
-  TriangleFilesWritten = 0;
 }
 ;
@@ -634 +585 @@
       cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
   }
-  cout << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
 }
 ;
@@ -1102 +1052 @@
   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
-    if (AddFlag) {
-      delete(BoundaryPoints);
-      BoundaryPoints = new LinkedCell(this, 5.);
-      AddFlag = false;
-    }
+    LinkedCell BoundaryPoints(this, 5.);
     Walker = cloud->GetPoint();
     *out << Verbose(2) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    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;
+    triangles = FindClosestTrianglesToPoint(out, Walker->node, &BoundaryPoints);
+    if (triangles == NULL) {
+      *out << Verbose(1) << "No triangles found, probably a tesselation point itself." << endl;
       cloud->GoToNext();
       continue;
     } else {
+      BTS = triangles->front();
     }
     *out << Verbose(2) << "Closest triangle is " << *BTS << "." << endl;
@@ -1147 +1090 @@
         // add Walker to boundary points
         *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
-        AddFlag = true;
         if (AddBoundaryPoint(Walker,0))
           NewPoint = BPS[0];
@@ -1238 +1180 @@
   } else {
     delete TPS[n];
-    cout << Verbose(4) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
+    cout << Verbose(3) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
     TPS[n] = (InsertUnique.first)->second;
   }
@@ -1255 +1197 @@
 
   if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine = a->lines.find(b->node->nr);
+    LineMap::iterator FindLine;
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-    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++) {
+
+    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(4) << "Using existing line " << *FindLine->second << endl;
+        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
 
         BPS[0] = FindLine->second->endpoints[0];
@@ -1291 +1232 @@
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
 {
-  cout << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
+  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
@@ -1301 +1242 @@
 };
 
-/** Function adds triangle to global list.
- * Furthermore, the triangle receives the next free id and id counter \a TrianglesOnBoundaryCount is increased.
+/** 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.
  */
 void Tesselation::AddTesselationTriangle()
@@ -1311 +1252 @@
   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
@@ -1351 +1272 @@
           cout << Verbose(5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
           RemoveTesselationLine(triangle->lines[i]);
-      } 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
+          triangle->lines[i] = NULL;
+      } else
+        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle." << endl;
     } else
       cerr << "ERROR: This line " << i << " has already been free'd." << endl;
@@ -1383 +1294 @@
   if (line == NULL)
     return;
-  // get other endpoint number for finding copies of same line
+  // get other endpoint number of finding copies of same line
   if (line->endpoints[1] != NULL)
     Numbers[0] = line->endpoints[1]->Nr;
@@ -1410 +1321 @@
         cout << Verbose(5) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         RemoveTesselationPoint(line->endpoints[i]);
-      } 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
+        line->endpoints[i] = NULL;
+      } else
+        cout << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to." << endl;
     } else
       cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
@@ -1483 +1390 @@
           }
           // 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;
         }
       }
@@ -1493 +1400 @@
   *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;
 };
 
@@ -1603 +1461 @@
   CandidateList *OptCandidates = new CandidateList();
   for (int k=0;k<NDIM;k++) {
-    Oben.Zero();
     Oben.x[k] = 1.;
     FirstPoint = MaxPoint[k];
@@ -1612 +1469 @@
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
 
-    FindSecondPointForTesselation(FirstPoint, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(FirstPoint, NULL, 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);
@@ -1637 +1496 @@
 
     // 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);
@@ -1712 +1569 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC)
 {
   cout << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
@@ -1747 +1604 @@
     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
@@ -1834 +1691 @@
         result = false;
       }
-    } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
+    } else if (existentTrianglesCount == 1) { // 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);
@@ -1875 +1732 @@
     // 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);
-    }
-
   }
 
@@ -1898 +1750 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return NULL - convex, otherwise endpoint that makes it concave
+ * \return NULL - concave, otherwise endpoint that makes it concave
  */
 class BoundaryPointSet *Tesselation::IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base)
@@ -1975 +1827 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return volume change due to flipping (0 - then no flipped occured)
- */
-double Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
+ * \return true - line was changed, false - same line as before
+ */
+bool Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
 {
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint[2];
-  double volume;
 
   int m=0;
@@ -2001 +1852 @@
   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
@@ -2018 +1866 @@
     if (Base->triangles.size() < 2) {
       *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
-      return 0.;
+      return false;
     }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
@@ -2028 +1876 @@
     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;
-      // calculate volume summand as a general tetraeder
-      return volume;
+      FlipBaseline(out, Base);
+      return true;
     } else {  // Base higher than OtherBase -> do nothing
       *out << Verbose(2) << "REJECT: Base line is higher: Nothing to do." << endl;
-      return 0.;
-    }
-  }
+      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;
 };
 
@@ -2042 +1930 @@
  * \param *out output stream for debugging
  * \param *Base line to be flipped
- * \return pointer to allocated new baseline - flipping successful, NULL - something went awry
- */
-class BoundaryLineSet * Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
+ * \return true - flipping successful, false - something went awry
+ */
+bool Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
 {
   class BoundaryLineSet *OldLines[4], *NewLine;
@@ -2058 +1946 @@
   if (Base->triangles.size() < 2) {
     *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
-    return NULL;
+    return false;
   }
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
@@ -2094 +1982 @@
   if (i<4) {
     *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-    return NULL;
+    return false;
   }
   for (int j=0;j<4;j++)
     if (OldLines[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-      return NULL;
+      return false;
     }
   for (int j=0;j<2;j++)
     if (OldPoints[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
-      return NULL;
+      return false;
     }
 
@@ -2136 +2024 @@
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
     BTS->GetNormalVector(BaseLineNormal);
-    AddTesselationTriangle(OldTriangleNrs[0]);
+    TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[0], BTS));
     *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
 
@@ -2144 +2032 @@
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
     BTS->GetNormalVector(BaseLineNormal);
-    AddTesselationTriangle(OldTriangleNrs[1]);
+    TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[1], BTS));
     *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 NULL;
+    return false;
   }
 
   *out << Verbose(1) << "End of FlipBaseline" << endl;
-  return NewLine;
+  return true;
 };
 
@@ -2158 +2046 @@
 /** 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
@@ -2164 +2053 @@
  * \param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
+void Tesselation::FindSecondPointForTesselation(TesselPoint* a, TesselPoint* Candidate, 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;
@@ -2303 +2191 @@
   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
@@ -2328 +2216 @@
     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 
       // check SearchDirection
@@ -2505 +2392 @@
   TesselPoint *trianglePoints[3];
   TesselPoint *SecondPoint = NULL;
-  list<BoundaryTriangleSet*> *triangles = NULL;
 
   if (LinesOnBoundary.empty()) {
@@ -2516 +2402 @@
   // check whether closest point is "too close" :), then it's inside
   if (trianglePoints[0] == NULL) {
-    *out << Verbose(2) << "Is the only point, no one else is closeby." << endl;
+    *out << Verbose(1) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
   if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
-    *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);
-  }
+    *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);
   
   if (triangles->empty()) {
-    *out << Verbose(0) << "ERROR: There is no nearest triangle. Please check the tesselation structure.";
-    delete(triangles);
+    *out << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
     return NULL;
   } else
@@ -2581 +2443 @@
   class BoundaryTriangleSet *result = NULL;
   list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(out, x, LC);
-  Vector Center;
 
   if (triangles == NULL)
     return NULL;
 
-  if (triangles->size() == 1) { // there is no degenerate case
+  if (x->ScalarProduct(&triangles->front()->NormalVector) < 0)
+    result = triangles->back();
+  else
     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;
@@ -2616 +2466 @@
 {
   class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, &Point, LC);
-  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;
+  if (result == NULL)
+    return true;
+  if (Point.ScalarProduct(&result->NormalVector) < 0)
+    return true;
+  else
     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;
-  } else {
-    *out << Verbose(1) << Point << " is NOT an inner point." << endl;
-    return false;
-  }
 }
 
@@ -2645 +2483 @@
 bool Tesselation::IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC)
 {
-  return IsInnerPoint(out, *(Point->node), 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;
 }
 
@@ -2652 +2496 @@
  * @param *Point of which get all connected points
  *
- * @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*>;
+ * @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*>;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -2666 +2510 @@
     ReferencePoint = PointRunner->second;
   } else {
-    *out << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    *out << Verbose(2) << "getCircleOfConnectedPoints() 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;
@@ -2677 +2521 @@
   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(5) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
-     connectedPoints->insert(current);
-   }
-
-   findLines++;
+    takePoint = false;
+
+    if (findLines->second->endpoints[0]->Nr == Point->nr) {
+      takePoint = true;
+      current = findLines->second->endpoints[1]->node;
+    } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
+      takePoint = true;
+      current = findLines->second->endpoints[0]->node;
+    }
+    
+    if (takePoint) {
+      *out << Verbose(3) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is taken into the circle." << endl;
+      connectedPoints->push_back(current);
+    }
+
+    findLines++;
   }
 
@@ -2699 +2543 @@
     return NULL;
   }
-
   return connectedPoints;
-};
-
-
-/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
+}
+
+/** Gets the two neighbouring points with respect to a reference line to the provided 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
@@ -2711 +2553 @@
  *
  * @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 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)
+ * @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)
 {
   map<double, TesselPoint*> anglesOfPoints;
-  set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(out, Point);
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
-  Vector center;
-  Vector PlaneNormal;
-  Vector AngleZero;
-  Vector OrthogonalVector;
-  Vector helper;
+  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;
+  }
 
   // calculate central point
-  for (set<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
+  for (list<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //*out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -2741 +2586 @@
 
   // construct one orthogonal vector
-  if (Reference != NULL)
-    AngleZero.CopyVector(Reference);
-  else
-    AngleZero.CopyVector((*connectedPoints->begin())->node);
+  AngleZero.CopyVector(Reference);
   AngleZero.SubtractVector(Point->node);
   AngleZero.ProjectOntoPlane(&PlaneNormal);
@@ -2752 +2594 @@
 
   // go through all connected points and calculate angle
-  for (set<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
+  for (list<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(3) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
+    *out << Verbose(2) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
     anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
   }
 
-  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
-    connectedCircle->push_back(AngleRunner->second);
-  }
-
-  delete(connectedPoints);
-  return connectedCircle;
+  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;
 }
 
-/** 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 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).
+ * We create new triangles and remove the old ones connected to the point.
  * \param *out output stream for debugging
  * \param *point point to be removed
@@ -2983 +2625 @@
   class BoundaryLineSet *line = NULL;
   class BoundaryTriangleSet *triangle = NULL;
-  Vector OldPoint, NormalVector;
+  Vector OldPoint, TetraederVector[3];
   double volume = 0;
+  int *numbers = NULL;
   int count = 0;
+  int i;
 
   if (point == NULL) {
@@ -3001 +2645 @@
     return 0.;
   }
-
-  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(out, point->node);
-  list<TesselPoint*> *connectedPath = NULL;
-
-  // gather all triangles
+  list<TesselPoint*> *CircleofPoints = GetCircleOfConnectedPoints(out, point->node);
+
+  // remove all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  map<class BoundaryTriangleSet *, int> Candidates;
-  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+  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++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
-    }
-  }
-
-  // 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;
+      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
+          }
+        }
+      } 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;
+      }
+    }
+  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;
             }
-        }
-        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();
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << 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;
+
+  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;
             }
           }
-          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.);
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
       }
 
-      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;
-};
-
-
+  return closestPoint;
+};
 
 /**
@@ -3230 +2939 @@
               FindTriangle = FindLine->second->triangles.begin();
               for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+                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])
+                  )
+                ) {
                   result->push_back(FindTriangle->second);
                 }
@@ -3248 +2970 @@
 
 /**
- * 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.
  *
@@ -3290 +2975 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedTriangles()
-{
-  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
-  map<int, int> * DegeneratedTriangles = new map<int, int>;
-
-  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) );
+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;
+            }
+          }
         }
       }
     }
   }
-  delete(DegeneratedLines);
-
-  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
+
+  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;
 
@@ -3333 +3026 @@
 void Tesselation::RemoveDegeneratedTriangles()
 {
-  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
+  map<int, int> DegeneratedTriangles = FindAllDegeneratedTriangles();
+
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles.begin();
+    TriangleKeyRunner != DegeneratedTriangles.end(); ++TriangleKeyRunner
   ) {
-    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;
+    BoundaryTriangleSet *triangle = TrianglesOnBoundary.find(TriangleKeyRunner->first)->second,
+      *partnerTriangle = TrianglesOnBoundary.find(TriangleKeyRunner->second)->second;
 
     bool trianglesShareLine = false;
@@ -3364 +3044 @@
       && (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
@@ -3408 +3055 @@
     }
   }
-  delete(DegeneratedTriangles);
-
-  cout << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl;
-  cout << Verbose(1) << "End of RemoveDegeneratedTriangles" << endl;
 }
 
-/** 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++;
-};
+/** 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;
+}
+