Changes in src/tesselation.cpp [244a84:125b3c]

File:

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

src/tesselation.cpp

@@ -35 +35 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
   LinesCount(0),
   node(Walker),
@@ -61 +61 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -105 +105 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
 {
         Info FunctionInfo(__func__);
@@ -115 +115 @@
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
-  // set skipped to false
-  skipped = false;
-  // clear triangles list
-  Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
-};
-
-/** Constructor of BoundaryLineSet with two endpoints.
- * Adds line automatically to each endpoints' LineMap
- * \param *Point1 first boundary point
- * \param *Point2 second boundary point
- * \param number number of the list
- */
-BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number)
-{
-  Info FunctionInfo(__func__);
-  // set number
-  Nr = number;
-  // set endpoints in ascending order
-  SetEndpointsOrdered(endpoints, Point1, Point2);
-  // add this line to the hash maps of both endpoints
-  Point1->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
-  Point2->AddLine(this); //
   // set skipped to false
   skipped = false;
@@ -193 +171 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
         Info FunctionInfo(__func__);
@@ -204 +182 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -219 +197 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion() const
+bool BoundaryLineSet::CheckConvexityCriterion()
 {
         Info FunctionInfo(__func__);
@@ -243 +221 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
     //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     NormalCheck.AddVector(&runner->second->NormalVector);
@@ -286 +264 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -299 +277 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -339 +317 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
   Nr(number)
 {
@@ -398 +376 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
 {
         Info FunctionInfo(__func__);
@@ -410 +388 @@
 };
 
-/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
+/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * 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
@@ -422 +400 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
-{
-  Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
+{
+        Info FunctionInfo(__func__);
   Vector CrossPoint;
   Vector helper;
@@ -433 +411 @@
   }
 
-  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
-  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
-  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
-
-  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
-    Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl;
-    return true;
-  }   else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
-    Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl;
-    return true;
-  }   else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {
-    Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl;
-    return true;
-  }
   // Calculate cross point between one baseline and the line from the third endpoint to intersection
   int i=0;
@@ -453 +417 @@
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
-      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-      const double s = CrossPoint.ScalarProduct(&helper)/helper.NormSquared();
-      Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl;
-      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
-        Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl;
-        i=4;
-        break;
-      }
+    } else 
       i++;
-    } else 
+    if (i>2)
       break;
-  } while (i<3);
+  } while (CrossPoint.NormSquared() < MYEPSILON);
   if (i==3) {
-    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
+    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
+  }
+  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+
+  // check whether intersection is inside or not by comparing length of intersection and length of cross point
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
     return true;
-  } else {
-    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
+  } else { // outside!
+    Intersection->Zero();
     return false;
   }
-};
-
-/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
- * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * Thus 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.
- * \param *x point
- * \param *ClosestPoint desired closest point inside triangle to \a *x, is absolute vector
- * \return Distance squared between \a *x and closest point inside triangle
- */
-double BoundaryTriangleSet::GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const
-{
-  Info FunctionInfo(__func__);
-  Vector Direction;
-
-  // 1. get intersection with plane
-  Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
-  GetCenter(&Direction);
-  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
-    ClosestPoint->CopyVector(x);
-  }
-
-  // 2. Calculate in plane part of line (x, intersection)
-  Vector InPlane;
-  InPlane.CopyVector(x);
-  InPlane.SubtractVector(ClosestPoint);  // points from plane intersection to straight-down point
-  InPlane.ProjectOntoPlane(&NormalVector);
-  InPlane.AddVector(ClosestPoint);
-
-  Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl;
-  Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl;
-  Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl;
-
-  // Calculate cross point between one baseline and the desired point such that distance is shortest
-  double ShortestDistance = -1.;
-  bool InsideFlag = false;
-  Vector CrossDirection[3];
-  Vector CrossPoint[3];
-  Vector helper;
-  for (int i=0;i<3;i++) {
-    // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
-    Direction.CopyVector(endpoints[(i+1)%3]->node->node);
-    Direction.SubtractVector(endpoints[i%3]->node->node);
-    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
-    CrossDirection[i].CopyVector(&CrossPoint[i]);
-    CrossDirection[i].SubtractVector(&InPlane);
-    CrossPoint[i].SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-    const double s = CrossPoint[i].ScalarProduct(&Direction)/Direction.NormSquared();
-    Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl;
-    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
-      CrossPoint[i].AddVector(endpoints[i%3]->node->node);  // make cross point absolute again
-      Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i%3]->node->node << " and " << *endpoints[(i+1)%3]->node->node << "." << endl;
-      const double distance = CrossPoint[i].DistanceSquared(x);
-      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
-        ShortestDistance = distance;
-        ClosestPoint->CopyVector(&CrossPoint[i]);
-      }
-    } else
-      CrossPoint[i].Zero();
-  }
-  InsideFlag = true;
-  for (int i=0;i<3;i++) {
-    const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+1)%3]);
-    const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+2)%3]);;
-    if ((sign > -MYEPSILON) && (othersign > -MYEPSILON))  // have different sign
-      InsideFlag = false;
-  }
-  if (InsideFlag) {
-    ClosestPoint->CopyVector(&InPlane);
-    ShortestDistance = InPlane.DistanceSquared(x);
-  } else {  // also check endnodes
-    for (int i=0;i<3;i++) {
-      const double distance = x->DistanceSquared(endpoints[i]->node->node);
-      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
-        ShortestDistance = distance;
-        ClosestPoint->CopyVector(endpoints[i]->node->node);
-      }
-    }
-  }
-  Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl;
-  return ShortestDistance;
 };
 
@@ -562 +440 @@
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -575 +453 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
 {
         Info FunctionInfo(__func__);
@@ -588 +466 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
 {
         Info FunctionInfo(__func__);
@@ -601 +479 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
-{
-        Info FunctionInfo(__func__);
-        Log() << Verbose(1) << "INFO: Checking " << Points[0] << ","  << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl;
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
+{
+        Info FunctionInfo(__func__);
   return (((endpoints[0] == Points[0])
             || (endpoints[0] == Points[1])
@@ -624 +501 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
 {
         Info FunctionInfo(__func__);
@@ -646 +523 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
 {
         Info FunctionInfo(__func__);
@@ -663 +540 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector * const center) const
+void BoundaryTriangleSet::GetCenter(Vector *center)
 {
         Info FunctionInfo(__func__);
@@ -670 +547 @@
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
-  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -939 +815 @@
 TesselPoint::TesselPoint()
 {
-  //Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
   node = NULL;
   nr = -1;
@@ -949 +825 @@
 TesselPoint::~TesselPoint()
 {
-  //Info FunctionInfo(__func__);
+  Info FunctionInfo(__func__);
 };
 
@@ -976 +852 @@
 PointCloud::PointCloud()
 {
-        //Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
 };
 
@@ -983 +859 @@
 PointCloud::~PointCloud()
 {
-        //Info FunctionInfo(__func__);
+        Info FunctionInfo(__func__);
 };
 
@@ -1174 +1050 @@
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
  */
-void Tesselation::GuessStartingTriangle()
+void
+Tesselation::GuessStartingTriangle()
 {
         Info FunctionInfo(__func__);
@@ -1545 +1422 @@
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  TriangleList *triangles = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
@@ -1560 +1437 @@
     Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
@@ -2461 +2338 @@
 
   // fill the set of neighbours
-  TesselPointSet SetOfNeighbours;
+  Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);
+  Center.SubtractVector(TurningPoint->node);
+  set<TesselPoint*> SetOfNeighbours;
   SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
-  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
 
   // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
-  Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl;
-  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
-    Log() << Verbose(0) << **TesselRunner << endl;
   TesselPointList::iterator Runner = connectedClosestPoints->begin();
   TesselPointList::iterator Sprinter = Runner;
@@ -2480 +2356 @@
     AddTesselationPoint((*Sprinter), 2);
 
+
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -2496 +2373 @@
     Runner = Sprinter;
     Sprinter++;
-    Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl;
-    if (Sprinter != connectedClosestPoints->end())
-      Log() << Verbose(0) << " There are still more triangles to add." << endl;
   }
   delete(connectedClosestPoints);
@@ -3083 +2957 @@
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  PointMap OrderMap;
-  PointTestPair OrderTest;
+  map<int, class BoundaryPointSet *> OrderMap;
+  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
@@ -3104 +2978 @@
 };
 
-/** Finds the boundary points that are closest to a given Vector \a *x.
- * \param *out output stream for debugging
- * \param *x Vector to look from
- * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
- */
-DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
-{
-  Info FunctionInfo(__func__);
-  PointMap::const_iterator FindPoint;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  if (LinesOnBoundary.empty()) {
-    eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl;
-    return NULL;
-  }
-
-  // gather all points close to the desired one
-  LC->SetIndexToVector(x); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  DistanceToPointMap * points = new DistanceToPointMap;
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //Log() << Verbose(1) << 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]++) {
-        const LinkedNodes *List = LC->GetCurrentCell();
-        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            FindPoint = PointsOnBoundary.find((*Runner)->nr);
-            if (FindPoint != PointsOnBoundary.end()) {
-              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
-              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
-            }
-          }
-        } else {
-          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-
-  // check whether we found some points
-  if (points->empty()) {
-    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
-    delete(points);
-    return NULL;
-  }
-  return points;
-};
-
-/** Finds the boundary line that is closest to a given Vector \a *x.
- * \param *out output stream for debugging
- * \param *x Vector to look from
- * \return closest BoundaryLineSet or NULL in degenerate case.
- */
-BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
-{
-  Info FunctionInfo(__func__);
-
-  // get closest points
-  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
-  if (points == NULL) {
-    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
-    return NULL;
-  }
-
-  // for each point, check its lines, remember closest
-  Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl;
-  BoundaryLineSet *ClosestLine = NULL;
-  double MinDistance = -1.;
-  Vector helper;
-  Vector Center;
-  Vector BaseLine;
-  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
-    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
-      // calculate closest point on line to desired point
-      helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-      helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
-      helper.Scale(0.5);
-      Center.CopyVector(x);
-      Center.SubtractVector(&helper);
-      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-      Center.ProjectOntoPlane(&BaseLine);
-      const double distance = Center.NormSquared();
-      if ((ClosestLine == NULL) || (distance < MinDistance)) {
-        // additionally calculate intersection on line (whether it's on the line section or not)
-        helper.CopyVector(x);
-        helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
-        helper.SubtractVector(&Center);
-        const double lengthA = helper.ScalarProduct(&BaseLine);
-        helper.CopyVector(x);
-        helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-        helper.SubtractVector(&Center);
-        const double lengthB = helper.ScalarProduct(&BaseLine);
-        if (lengthB*lengthA < 0) {  // if have different sign
-          ClosestLine = LineRunner->second;
-          MinDistance = distance;
-          Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
-        } else {
-          Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
-        }
-      } else {
-        Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
-      }
-    }
-  }
-  delete(points);
-  // check whether closest line is "too close" :), then it's inside
-  if (ClosestLine == NULL) {
-    Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
-    return NULL;
-  }
-  return ClosestLine;
-};
-
-
 /** Finds the triangle that is closest to a given Vector \a *x.
  * \param *out output stream for debugging
  * \param *x Vector to look from
- * \return BoundaryTriangleSet of nearest triangle or NULL.
- */
-TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
-{
-        Info FunctionInfo(__func__);
-
-        // get closest points
-        DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
-  if (points == NULL) {
-    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+ * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
+ */
+list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const
+{
+        Info FunctionInfo(__func__);
+  TesselPoint *trianglePoints[3];
+  TesselPoint *SecondPoint = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
+
+  if (LinesOnBoundary.empty()) {
+    eLog() << Verbose(1) << "Error: There is no tesselation structure to compare the point with, please create one first.";
     return NULL;
   }
-
-  // for each point, check its lines, remember closest
-  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
-  LineSet ClosestLines;
-  double MinDistance = 1e+16;
-  Vector BaseLineIntersection;
-  Vector Center;
-  Vector BaseLine;
-  Vector BaseLineCenter;
-  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
-    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
-
-      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-      const double lengthBase = BaseLine.NormSquared();
-
-      BaseLineIntersection.CopyVector(x);
-      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
-      const double lengthEndA = BaseLineIntersection.NormSquared();
-
-      BaseLineIntersection.CopyVector(x);
-      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-      const double lengthEndB = BaseLineIntersection.NormSquared();
-
-      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
-        const double lengthEnd = Min(lengthEndA, lengthEndB);
-        if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
-          ClosestLines.clear();
-          ClosestLines.insert(LineRunner->second);
-          MinDistance = lengthEnd;
-          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl;
-        } else if  (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
-          ClosestLines.insert(LineRunner->second);
-          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl;
-        } else { // line is worse
-          Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
-        }
-      } else { // intersection is closer, calculate
-        // calculate closest point on line to desired point
-        BaseLineIntersection.CopyVector(x);
-        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-        Center.CopyVector(&BaseLineIntersection);
-        Center.ProjectOntoPlane(&BaseLine);
-        BaseLineIntersection.SubtractVector(&Center);
-        const double distance = BaseLineIntersection.NormSquared();
-        if (Center.NormSquared() > BaseLine.NormSquared()) {
-          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
-        }
-        if ((ClosestLines.empty()) || (distance < MinDistance)) {
-          ClosestLines.insert(LineRunner->second);
-          MinDistance = distance;
-          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl;
-        } else {
-          Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
-        }
-      }
-    }
-  }
-  delete(points);
-
-  // check whether closest line is "too close" :), then it's inside
-  if (ClosestLines.empty()) {
+  Log() << Verbose(1) << "Finding closest Tesselpoint to " << *x << " ... " << endl;
+  trianglePoints[0] = FindClosestPoint(x, SecondPoint, LC);
+  
+  // check whether closest point is "too close" :), then it's inside
+  if (trianglePoints[0] == NULL) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
-  TriangleList * candidates = new TriangleList;
-  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
-    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
-    candidates->push_back(Runner->second);
-  }
-  return candidates;
+  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
+    Log() << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
+    PointMap::const_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 {
+        eLog() << Verbose(1) << "I cannot find a boundary point to the tessel point " << *trianglePoints[0] << "." << endl;
+        return NULL;
+      }
+    }
+  } else {
+    set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(trianglePoints[0]);
+    TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(connectedPoints, trianglePoints[0], x);
+    delete(connectedPoints);
+    if (connectedClosestPoints != NULL) {
+      trianglePoints[1] = connectedClosestPoints->front();
+      trianglePoints[2] = connectedClosestPoints->back();
+      for (int i=0;i<3;i++) {
+        if (trianglePoints[i] == NULL) {
+          eLog() << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+        }
+        //Log() << Verbose(1) << "List of triangle points:" << endl;
+        //Log() << Verbose(2) << *trianglePoints[i] << endl;
+      }
+
+      triangles = FindTriangles(trianglePoints);
+      Log() << Verbose(1) << "List of possible triangles:" << endl;
+      for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+        Log() << Verbose(2) << **Runner << endl;
+
+      delete(connectedClosestPoints);
+    } else {
+      triangles = NULL;
+      eLog() << Verbose(2) << "There is no circle of connected points!" << endl;
+    }
+  }
+  
+  if ((triangles == NULL) || (triangles->empty())) {
+    eLog() << Verbose(1) << "There is no nearest triangle. Please check the tesselation structure.";
+    delete(triangles);
+    return NULL;
+  } else
+    return triangles;
 };
 
@@ -3317 +3067 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
 {
         Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
-  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
-  TriangleList candidates;
+  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
   Vector Center;
-  Vector helper;
-
-  if ((triangles == NULL) || (triangles->empty()))
+
+  if (triangles == NULL)
     return NULL;
 
-  // go through all and pick the one with the best alignment to x
-  double MinAlignment = 2.*M_PI;
-  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
-    (*Runner)->GetCenter(&Center);
-    helper.CopyVector(x);
-    helper.SubtractVector(&Center);
-    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
-    if (Alignment < MinAlignment) {
-      result = *Runner;
-      MinAlignment = Alignment;
-      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
-    } else {
-      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
+  if (triangles->size() == 1) { // there is no degenerate case
+    result = triangles->front();
+    Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
+  } else {
+    result = triangles->front();
+    result->GetCenter(&Center);
+    Center.SubtractVector(x);
+    Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
+    if (Center.ScalarProduct(&result->NormalVector) < 0) {
+      result = triangles->back();
+      Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
+      if (Center.ScalarProduct(&result->NormalVector) < 0) {
+        eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
+      }
     }
   }
   delete(triangles);
-
   return result;
 };
 
-/** Checks whether the provided Vector is within the Tesselation structure.
- * Basically calls Tesselation::GetDistanceToSurface() and checks the sign of the return value.
- * @param point of which to check the position
- * @param *LC LinkedCell structure
- *
- * @return true if the point is inside the Tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
-{
-  return (GetDistanceSquaredToSurface(Point, LC) < MYEPSILON);
-}
-
-/** Returns the distance to the surface given by the tesselation.
- * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
- * towards or away from the given \a &Point. Additionally, we check whether it's normal to the normal vector, i.e. on the
- * closest triangle's plane. Then, we have to check whether \a Point is inside the triangle or not to determine whether it's
- * an inside or outside point. This is done by calling BoundaryTriangleSet::GetIntersectionInsideTriangle().
- * In the end we additionally find the point on the triangle who was smallest distance to \a Point:
- *  -# Separate distance from point to center in vector in NormalDirection and on the triangle plane.
- *  -# Check whether vector on triangle plane points inside the triangle or crosses triangle bounds.
- *  -# If inside, take it to calculate closest distance
- *  -# If not, take intersection with BoundaryLine as distance
- *
- * @note distance is squared despite it still contains a sign to determine in-/outside!
+/** Checks whether the provided Vector is within the tesselation structure.
  *
  * @param point of which to check the position
  * @param *LC LinkedCell structure
  *
- * @return >0 if outside, ==0 if on surface, <0 if inside
- */
-double Tesselation::GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const
-{
-  Info FunctionInfo(__func__);
+ * @return true if the point is inside the tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
+{
+        Info FunctionInfo(__func__);
+  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
   Vector Center;
-  Vector helper;
-  Vector DistanceToCenter;
-  Vector Intersection;
-  double distance = 0.;
-
-  if (triangle == NULL) {// is boundary point or only point in point cloud?
-    Log() << Verbose(1) << "No triangle given!" << endl;
-    return -1.;
+
+  if (result == NULL) {// is boundary point or only point in point cloud?
+    Log() << Verbose(1) << Point << " is the only point in vicinity." << endl;
+    return false;
+  }
+
+  result->GetCenter(&Center);
+  Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
+  Center.SubtractVector(&Point);
+  Log() << Verbose(2) << "INFO: Vector from center to point to test is " << Center << "." << endl;
+  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
+    Log() << Verbose(1) << Point << " is an inner point." << endl;
+    return true;
   } else {
-    Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl;
-  }
-
-  triangle->GetCenter(&Center);
-  Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
-  DistanceToCenter.CopyVector(&Center);
-  DistanceToCenter.SubtractVector(&Point);
-  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl;
-
-  // check whether we are on boundary
-  if (fabs(DistanceToCenter.ScalarProduct(&triangle->NormalVector)) < MYEPSILON) {
-    // calculate whether inside of triangle
-    DistanceToCenter.CopyVector(&Point);
-    Center.CopyVector(&Point);
-    Center.SubtractVector(&triangle->NormalVector); // points towards MolCenter
-    DistanceToCenter.AddVector(&triangle->NormalVector); // points outside
-    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
-    if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
-      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
-      return 0.;
-    } else {
-      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
-      return false;
-    }
-  } else {
-    // calculate smallest distance
-    distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection);
-    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
-
-    // then check direction to boundary
-    if (DistanceToCenter.ScalarProduct(&triangle->NormalVector) > MYEPSILON) {
-      Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
-      return -distance;
-    } else {
-      Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl;
-      return +distance;
-    }
-  }
-};
-
-/** Calculates distance to a tesselated surface.
- * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
- * \param &Point point to calculate distance from
- * \param *LC needed for finding closest points fast
- * \return distance squared to closest point on surface
- */
-double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
-{
-  BoundaryTriangleSet *triangle = FindClosestTriangleToVector(&Point, LC);
-  const double distance = GetDistanceSquaredToTriangle(Point, triangle);
-  return Min(distance, LC->RADIUS);
-};
+    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
+    return false;
+  }
+}
+
+/** Checks whether the provided TesselPoint is within the tesselation structure.
+ *
+ * @param *Point of which to check the position
+ * @param *LC Linked Cell structure
+ *
+ * @return true if the point is inside the tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const
+{
+        Info FunctionInfo(__func__);
+  return IsInnerPoint(*(Point->node), LC);
+}
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -3451 +3147 @@
  * @return set of the all points linked to the provided one
  */
-TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
-{
-        Info FunctionInfo(__func__);
-        TesselPointSet *connectedPoints = new TesselPointSet;
+set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -3495 +3191 @@
   }
 
-  if (connectedPoints->empty()) { // if have not found any points
+  if (connectedPoints->size() == 0) { // if have not found any points
     eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl;
     return NULL;
@@ -3516 +3212 @@
  * @return list of the all points linked to the provided one
  */
-TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  TesselPointList *connectedCircle = new TesselPointList;
+  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  Vector center;
   Vector PlaneNormal;
   Vector AngleZero;
   Vector OrthogonalVector;
   Vector helper;
-  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
-  TriangleList *triangles = NULL;
 
   if (SetOfNeighbours == NULL) {
@@ -3535 +3230 @@
 
   // calculate central point
-  triangles = FindTriangles(TrianglePoints);
-  if ((triangles != NULL) && (!triangles->empty())) {
-    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-      PlaneNormal.AddVector(&(*Runner)->NormalVector);
-  } else {
-    eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl;
-    performCriticalExit();
-  }
-  PlaneNormal.Scale(1.0/triangles->size());
-  Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl;
+  for (set<TesselPoint*>::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+    center.AddVector((*TesselRunner)->node);
+  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
+  //  << "; scale factor " << 1.0/connectedPoints.size();
+  center.Scale(1.0/SetOfNeighbours->size());
+  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
+
+  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
+  PlaneNormal.CopyVector(Point->node);
+  PlaneNormal.SubtractVector(&center);
   PlaneNormal.Normalize();
+  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
 
   // construct one orthogonal vector
@@ -3571 +3267 @@
 
   // go through all connected points and calculate angle
-  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
@@ -3587 +3283 @@
 }
 
-/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
- * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
- * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
- * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
- * triangle we are looking for.
- *
- * @param *SetOfNeighbours all points for which the angle should be calculated
- * @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
- */
-TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
-{
-  Info FunctionInfo(__func__);
-  map<double, TesselPoint*> anglesOfPoints;
-  TesselPointList *connectedCircle = new TesselPointList;
-  Vector center;
-  Vector PlaneNormal;
-  Vector AngleZero;
-  Vector OrthogonalVector;
-  Vector helper;
-
-  if (SetOfNeighbours == NULL) {
-    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
-    delete(connectedCircle);
-    return NULL;
-  }
-
-  // check whether there's something to do
-  if (SetOfNeighbours->size() < 3) {
-    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
-      connectedCircle->push_back(*TesselRunner);
-    return connectedCircle;
-  }
-
-  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
-  // calculate central point
-
-  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
-  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
-  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
-  TesselB++;
-  TesselC++;
-  TesselC++;
-  int counter = 0;
-  while (TesselC != SetOfNeighbours->end()) {
-    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
-    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
-    counter++;
-    TesselA++;
-    TesselB++;
-    TesselC++;
-    PlaneNormal.AddVector(&helper);
-  }
-  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
-  //  << "; scale factor " << counter;
-  PlaneNormal.Scale(1.0/(double)counter);
-//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
-//
-//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
-//  PlaneNormal.CopyVector(Point->node);
-//  PlaneNormal.SubtractVector(&center);
-//  PlaneNormal.Normalize();
-  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
-
-  // construct one orthogonal vector
-  if (Reference != NULL) {
-    AngleZero.CopyVector(Reference);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
-  }
-  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
-    Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;
-    AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
-    if (AngleZero.NormSquared() < MYEPSILON) {
-      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
-      performCriticalExit();
-    }
-  }
-  Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
-  if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
-  else
-    OrthogonalVector.MakeNormalVector(&PlaneNormal);
-  Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
-
-  // go through all connected points and calculate angle
-  pair <map<double, TesselPoint*>::iterator, bool > InserterTest;
-  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper.CopyVector((*listRunner)->node);
-    helper.SubtractVector(Point->node);
-    helper.ProjectOntoPlane(&PlaneNormal);
-    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
-    if (angle > M_PI) // the correction is of no use here (and not desired)
-      angle = 2.*M_PI - angle;
-    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
-    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
-    if (!InserterTest.second) {
-      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
-      performCriticalExit();
-    }
-  }
-
-  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
-    connectedCircle->push_back(AngleRunner->second);
-  }
-
-  return connectedCircle;
-}
-
 /** Gets all points connected to the provided point by triangulation lines, ordered such that we walk along a closed path.
  *
@@ -3705 +3289 @@
  * @return list of the all points linked to the provided one
  */
-ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
         Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
-  TesselPointList *connectedPath = NULL;
+  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
@@ -3747 +3331 @@
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new TesselPointList;
+        connectedPath = new list<TesselPoint*>;
         triangle = NULL;
         CurrentLine = runner->second;
@@ -3821 +3405 @@
  * @return list of the closed paths
  */
-ListOfTesselPointList * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
-        Info FunctionInfo(__func__);
-  list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
-  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
-  TesselPointList *connectedPath = NULL;
-  TesselPointList *newPath = NULL;
+list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
+        Info FunctionInfo(__func__);
+  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point);
+  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPoint*> *newPath = NULL;
   int count = 0;
 
 
-  TesselPointList::iterator CircleRunner;
-  TesselPointList::iterator CircleStart;
-
-  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  list<TesselPoint*>::iterator CircleRunner;
+  list<TesselPoint*>::iterator CircleStart;
+
+  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
@@ -3843 +3427 @@
 
     // go through list, look for reappearance of starting Point and create list
-    TesselPointList::iterator Marker = CircleStart;
+    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
         Log() << Verbose(1) << count+1 << ". closed path consists of: ";
-        newPath = new TesselPointList;
-        TesselPointList::iterator CircleSprinter = Marker;
+        newPath = new list<TesselPoint*>;
+        list<TesselPoint*>::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
@@ -3883 +3467 @@
  * \return pointer to allocated list of triangles
  */
-TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
-{
-        Info FunctionInfo(__func__);
-        TriangleSet *connectedTriangles = new TriangleSet;
+set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+        Info FunctionInfo(__func__);
+  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
 
   if (Point == NULL) {
@@ -3935 +3519 @@
   }
 
-  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  TesselPointList *connectedPath = NULL;
+  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  list<TesselPoint*> *connectedPath = NULL;
 
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  TriangleMap Candidates;
+  map<class BoundaryTriangleSet *, int> Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
+      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
     }
   }
@@ -3953 +3537 @@
   count=0;
   NormalVector.Zero();
-  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl;
-    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
-    RemoveTesselationTriangle(Runner->second);
+  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
+    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->first);
     count++;
   }
   Log() << Verbose(1) << count << " triangles were removed." << endl;
 
-  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
-  list<TesselPointList *>::iterator ListRunner = ListAdvance;
-  TriangleMap::iterator NumberRunner = Candidates.begin();
-  TesselPointList::iterator StartNode, MiddleNode, EndNode;
+  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;
@@ -3978 +3562 @@
 
       // re-create all triangles by going through connected points list
-      LineList NewLines;
+      list<class BoundaryLineSet *> NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
@@ -4084 +3668 @@
       // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
       if (NewLines.size() > 1) {
-        LineList::iterator Candidate;
+        list<class BoundaryLineSet *>::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
-          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -4131 +3715 @@
  * Finds triangles belonging to the three provided points.
  *
- * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
+ * @param *Points[3] list, is expected to contain three points
  *
  * @return triangles which belong to the provided points, will be empty if there are none,
  *         will usually be one, in case of degeneration, there will be two
  */
-TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
-{
-        Info FunctionInfo(__func__);
-        TriangleList *result = new TriangleList;
+list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+        Info FunctionInfo(__func__);
+  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
   class BoundaryPointSet *TrianglePoints[3];
-  size_t NoOfWildcards = 0;
 
   for (int i = 0; i < 3; i++) {
-    if (Points[i] == NULL) {
-      NoOfWildcards++;
+    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      TrianglePoints[i] = FindPoint->second;
+    } else {
       TrianglePoints[i] = NULL;
-    } else {
-      PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
-      if (FindPoint != PointsOnBoundary.end()) {
-        TrianglePoints[i] = FindPoint->second;
-      } else {
-        TrianglePoints[i] = NULL;
-      }
-    }
-  }
-
-  switch (NoOfWildcards) {
-    case 0: // checks lines between the points in the Points for their adjacent triangles
-      for (int i = 0; i < 3; i++) {
-        if (TrianglePoints[i] != NULL) {
-          for (int j = i+1; j < 3; j++) {
-            if (TrianglePoints[j] != NULL) {
-              for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
-                  (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
-                  FindLine++) {
-                for (FindTriangle = FindLine->second->triangles.begin();
-                    FindTriangle != FindLine->second->triangles.end();
-                    FindTriangle++) {
-                  if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
-                    result->push_back(FindTriangle->second);
-                  }
-                }
+    }
+  }
+
+  // checks lines between the points in the Points for their adjacent triangles
+  for (int i = 0; i < 3; i++) {
+    if (TrianglePoints[i] != NULL) {
+      for (int j = i+1; j < 3; j++) {
+        if (TrianglePoints[j] != NULL) {
+          for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
+              (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
+              FindLine++) {
+            for (FindTriangle = FindLine->second->triangles.begin();
+                FindTriangle != FindLine->second->triangles.end();
+                FindTriangle++) {
+              if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+                result->push_back(FindTriangle->second);
               }
-              // Is it sufficient to consider one of the triangle lines for this.
-              return result;
             }
           }
+          // Is it sufficient to consider one of the triangle lines for this.
+          return result;
         }
       }
-      break;
-    case 1: // copy all triangles of the respective line
-    {
-      int i=0;
-      for (; i < 3; i++)
-        if (TrianglePoints[i] == NULL)
-          break;
-      for (FindLine = TrianglePoints[(i+1)%3]->lines.find(TrianglePoints[(i+2)%3]->node->nr); // is a multimap!
-          (FindLine != TrianglePoints[(i+1)%3]->lines.end()) && (FindLine->first == TrianglePoints[(i+2)%3]->node->nr);
-          FindLine++) {
-        for (FindTriangle = FindLine->second->triangles.begin();
-            FindTriangle != FindLine->second->triangles.end();
-            FindTriangle++) {
-          if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
-            result->push_back(FindTriangle->second);
-          }
-        }
-      }
-      break;
-    }
-    case 2: // copy all triangles of the respective point
-    {
-      int i=0;
-      for (; i < 3; i++)
-        if (TrianglePoints[i] != NULL)
-          break;
-      for (LineMap::const_iterator line = TrianglePoints[i]->lines.begin(); line != TrianglePoints[i]->lines.end(); line++)
-        for (TriangleMap::const_iterator triangle = line->second->triangles.begin(); triangle != line->second->triangles.end(); triangle++)
-          result->push_back(triangle->second);
-      result->sort();
-      result->unique();
-      break;
-    }
-    case 3: // copy all triangles
-    {
-      for (TriangleMap::const_iterator triangle = TrianglesOnBoundary.begin(); triangle != TrianglesOnBoundary.end(); triangle++)
-        result->push_back(triangle->second);
-      break;
-    }
-    default:
-      eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl;
-      performCriticalExit();
-      break;
+    }
   }
 
@@ -4267 +3800 @@
  *         in the list, once as key and once as value
  */
-IndexToIndex * Tesselation::FindAllDegeneratedLines()
+map<int, int> * Tesselation::FindAllDegeneratedLines()
 {
         Info FunctionInfo(__func__);
         UniqueLines AllLines;
-  IndexToIndex * DegeneratedLines = new IndexToIndex;
+  map<int, int> * DegeneratedLines = new map<int, int>;
 
   // sanity check
@@ -4292 +3825 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
-  IndexToIndex::iterator it;
+  map<int,int>::iterator it;
   for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
     const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
@@ -4311 +3844 @@
  *         in the list, once as key and once as value
  */
-IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
-{
-        Info FunctionInfo(__func__);
-  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
-  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
+map<int, int> * Tesselation::FindAllDegeneratedTriangles()
+{
+        Info FunctionInfo(__func__);
+  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
+  map<int, int> * DegeneratedTriangles = new map<int, int>;
 
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
@@ -4321 +3854 @@
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
 
-  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
@@ -4342 +3875 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
-  IndexToIndex::iterator it;
+  map<int,int>::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
@@ -4356 +3889 @@
 {
         Info FunctionInfo(__func__);
-  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
 
-  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
@@ -4449 +3982 @@
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
@@ -4616 +4149 @@
 
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
-  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
+  map <int, int> *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
@@ -4768 +4301 @@
   }
 
-  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
   Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
-  IndexToIndex::iterator it;
+  map<int,int>::iterator it;
   for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;