source: src/tesselation.cpp@ d3dff06

/*
 * tesselation.cpp
 *
 *  Created on: Aug 3, 2009
 *      Author: heber
 */

#include "Helpers/MemDebug.hpp"

#include <fstream>
#include <assert.h>

#include "helpers.hpp"
#include "info.hpp"
#include "linkedcell.hpp"
#include "log.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "triangleintersectionlist.hpp"
#include "vector.hpp"
#include "Line.hpp"
#include "vector_ops.hpp"
#include "verbose.hpp"
#include "Plane.hpp"
#include "Exceptions/LinearDependenceException.hpp"
#include "Helpers/Assert.hpp"

#include "Helpers/Assert.hpp"

class molecule;

// ======================================== Points on Boundary =================================

/** Constructor of BoundaryPointSet.
 */
BoundaryPointSet::BoundaryPointSet() :
  LinesCount(0), value(0.), Nr(-1)
{
  Info FunctionInfo(__func__);
  DoLog(1) && (Log() << Verbose(1) << "Adding noname." << endl);
}
;

/** Constructor of BoundaryPointSet with Tesselpoint.
 * \param *Walker TesselPoint this boundary point represents
 */
BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
  LinesCount(0), node(Walker), value(0.), Nr(Walker->nr)
{
  Info FunctionInfo(__func__);
  DoLog(1) && (Log() << Verbose(1) << "Adding Node " << *Walker << endl);
}
;

/** Destructor of BoundaryPointSet.
 * Sets node to NULL to avoid removing the original, represented TesselPoint.
 * \note When removing point from a class Tesselation, use RemoveTesselationPoint()
 */
BoundaryPointSet::~BoundaryPointSet()
{
  Info FunctionInfo(__func__);
  //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl;
  if (!lines.empty())
    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);
  node = NULL;
}
;

/** Add a line to the LineMap of this point.
 * \param *line line to add
 */
void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
{
  Info FunctionInfo(__func__);
  DoLog(1) && (Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." << endl);
  if (line->endpoints[0] == this) {
    lines.insert(LinePair(line->endpoints[1]->Nr, line));
  } else {
    lines.insert(LinePair(line->endpoints[0]->Nr, line));
  }
  LinesCount++;
}
;

/** output operator for BoundaryPointSet.
 * \param &ost output stream
 * \param &a boundary point
 */
ostream & operator <<(ostream &ost, const BoundaryPointSet &a)
{
  ost << "[" << a.Nr << "|" << a.node->getName() << " at " << *a.node->node << "]";
  return ost;
}
;

// ======================================== Lines on Boundary =================================

/** Constructor of BoundaryLineSet.
 */
BoundaryLineSet::BoundaryLineSet() :
  Nr(-1)
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 2; i++)
    endpoints[i] = NULL;
}
;

/** Constructor of BoundaryLineSet with two endpoints.
 * Adds line automatically to each endpoints' LineMap
 * \param *Point[2] array of two boundary points
 * \param number number of the list
 */
BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
{
  Info FunctionInfo(__func__);
  // set number
  Nr = number;
  // set endpoints in ascending order
  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
  // add this line to the hash maps of both endpoints
  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
  DoLog(0) && (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;
  // clear triangles list
  DoLog(0) && (Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl);
}
;

/** Destructor for BoundaryLineSet.
 * Removes itself from each endpoints' LineMap, calling RemoveTrianglePoint() when point not connected anymore.
 * \note When removing lines from a class Tesselation, use RemoveTesselationLine()
 */
BoundaryLineSet::~BoundaryLineSet()
{
  Info FunctionInfo(__func__);
  int Numbers[2];

  // get other endpoint number of finding copies of same line
  if (endpoints[1] != NULL)
    Numbers[0] = endpoints[1]->Nr;
  else
    Numbers[0] = -1;
  if (endpoints[0] != NULL)
    Numbers[1] = endpoints[0]->Nr;
  else
    Numbers[1] = -1;

  for (int i = 0; i < 2; i++) {
    if (endpoints[i] != NULL) {
      if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
        pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
        for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
          if ((*Runner).second == this) {
            //Log() << Verbose(0) << "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)) {
          //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
        }
      }
      if (endpoints[i]->lines.empty()) {
        //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
        if (endpoints[i] != NULL) {
          delete (endpoints[i]);
          endpoints[i] = NULL;
        }
      }
    }
  }
  if (!triangles.empty())
    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl);
}
;

/** Add triangle to TriangleMap of this boundary line.
 * \param *triangle to add
 */
void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
{
  Info FunctionInfo(__func__);
  DoLog(0) && (Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl);
  triangles.insert(TrianglePair(triangle->Nr, triangle));
}
;

/** Checks whether we have a common endpoint with given \a *line.
 * \param *line other line to test
 * \return true - common endpoint present, false - not connected
 */
bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
{
  Info FunctionInfo(__func__);
  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
    return true;
  else
    return false;
}
;

/** 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.
 * If greater/equal M_PI than we are convex.
 * \param *out output stream for debugging
 * \return true - triangles are convex, false - concave or less than two triangles connected
 */
bool BoundaryLineSet::CheckConvexityCriterion() const
{
  Info FunctionInfo(__func__);
  Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck;
  // get the two triangles
  if (triangles.size() != 2) {
    DoeLog(0) && (eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);
    return true;
  }
  // check normal vectors
  // have a normal vector on the base line pointing outwards
  //Log() << Verbose(0) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
  BaseLineCenter = (1./2.)*((*endpoints[0]->node->node) + (*endpoints[1]->node->node));
  BaseLine = (*endpoints[0]->node->node) - (*endpoints[1]->node->node);

  //Log() << Verbose(0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;

  BaseLineNormal.Zero();
  NormalCheck.Zero();
  double sign = -1.;
  int i = 0;
  class BoundaryPointSet *node = NULL;
  for (TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
    //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
    NormalCheck += runner->second->NormalVector;
    NormalCheck *= sign;
    sign = -sign;
    if (runner->second->NormalVector.NormSquared() > MYEPSILON)
      BaseLineNormal = runner->second->NormalVector;   // yes, copy second on top of first
    else {
      DoeLog(0) && (eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);
    }
    node = runner->second->GetThirdEndpoint(this);
    if (node != NULL) {
      //Log() << Verbose(0) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
      helper[i] = (*node->node->node) - BaseLineCenter;
      helper[i].MakeNormalTo(BaseLine);  // we want to compare the triangle's heights' angles!
      //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
      i++;
    } else {
      DoeLog(1) && (eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);
      return true;
    }
  }
  //Log() << Verbose(0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
  if (NormalCheck.NormSquared() < MYEPSILON) {
    DoLog(0) && (Log() << Verbose(0) << "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) {
    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);
    return true;
  } else {
    DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);
    return false;
  }
}

/** Checks whether point is any of the two endpoints this line contains.
 * \param *point point to test
 * \return true - point is of the line, false - is not
 */
bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 2; i++)
    if (point == endpoints[i])
      return true;
  return false;
}
;

/** Returns other endpoint of the line.
 * \param *point other endpoint
 * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
 */
class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
{
  Info FunctionInfo(__func__);
  if (endpoints[0] == point)
    return endpoints[1];
  else if (endpoints[1] == point)
    return endpoints[0];
  else
    return NULL;
}
;

/** output operator for BoundaryLineSet.
 * \param &ost output stream
 * \param &a boundary line
 */
ostream & operator <<(ostream &ost, const BoundaryLineSet &a)
{
  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->getName() << " at " << *a.endpoints[1]->node->node << "]";
  return ost;
}
;

// ======================================== Triangles on Boundary =================================

/** Constructor for BoundaryTriangleSet.
 */
BoundaryTriangleSet::BoundaryTriangleSet() :
  Nr(-1)
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 3; i++) {
    endpoints[i] = NULL;
    lines[i] = NULL;
  }
}
;

/** Constructor for BoundaryTriangleSet with three lines.
 * \param *line[3] lines that make up the triangle
 * \param number number of triangle
 */
BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
  Nr(number)
{
  Info FunctionInfo(__func__);
  // set number
  // set lines
  for (int i = 0; i < 3; i++) {
    lines[i] = line[i];
    lines[i]->AddTriangle(this);
  }
  // get ascending order of endpoints
  PointMap OrderMap;
  for (int i = 0; i < 3; i++) {
    // for all three lines
    for (int j = 0; j < 2; j++) { // for both endpoints
      OrderMap.insert(pair<int, class BoundaryPointSet *> (line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
      // and we don't care whether insertion fails
    }
  }
  // set endpoints
  int Counter = 0;
  DoLog(0) && (Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl);
  for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) {
    endpoints[Counter] = runner->second;
    DoLog(0) && (Log() << Verbose(0) << " " << *endpoints[Counter] << endl);
    Counter++;
  }
  ASSERT(Counter >= 3,"We have a triangle with only two distinct endpoints!");
};


/** Destructor of BoundaryTriangleSet.
 * Removes itself from each of its lines' LineMap and removes them if necessary.
 * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
 */
BoundaryTriangleSet::~BoundaryTriangleSet()
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 3; i++) {
    if (lines[i] != NULL) {
      if (lines[i]->triangles.erase(Nr)) {
        //Log() << Verbose(0) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
      }
      if (lines[i]->triangles.empty()) {
        //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
        delete (lines[i]);
        lines[i] = NULL;
      }
    }
  }
  //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl;
}
;

/** Calculates the normal vector for this triangle.
 * Is made unique by comparison with \a OtherVector to point in the other direction.
 * \param &OtherVector direction vector to make normal vector unique.
 */
void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
{
  Info FunctionInfo(__func__);
  // get normal vector
  NormalVector = Plane(*(endpoints[0]->node->node),
                       *(endpoints[1]->node->node),
                       *(endpoints[2]->node->node)).getNormal();

  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
  if (NormalVector.ScalarProduct(OtherVector) > 0.)
    NormalVector.Scale(-1.);
  DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl);
}
;

/** 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 *out output stream for debugging
 * \param *MolCenter offset vector of line
 * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
 * \param *Intersection intersection on plane on return
 * \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__);
  Vector CrossPoint;
  Vector helper;

  try {
    Line centerLine = makeLineThrough(*MolCenter, *x);
    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(centerLine);

    DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
    DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
    DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);

    if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
      return true;
    }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
      return true;
    }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
      DoLog(1) && (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;
    do {
      Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
      Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
      CrossPoint = line1.getIntersection(line2);
      helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
      CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
      const double s = CrossPoint.ScalarProduct(helper)/helper.NormSquared();
      DoLog(1) && (Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl);
      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
        DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
        return false;
      }
      i++;
    } while (i < 3);
    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
    return true;
  }
  catch (MathException &excp) {
    Log() << Verbose(1) << excp;
    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
    return false;
  }
}
;

/** Finds the point on the triangle to the point \a *x.
 * We call Vector::GetIntersectionWithPlane() with \a * and the center of the triangle to receive an intersection point.
 * Then we check the in-plane part (the part projected down onto plane). We check whether it crosses one of the
 * boundary lines. If it does, we return this intersection as closest point, otherwise the projected point down.
 * 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
  DoLog(1) && (Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl);
  GetCenter(&Direction);
  try {
    Line l = makeLineThrough(*x, Direction);
    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(l);
  }
  catch (MathException &excp) {
    (*ClosestPoint) = (*x);
  }

  // 2. Calculate in plane part of line (x, intersection)
  Vector InPlane = (*x) - (*ClosestPoint); // points from plane intersection to straight-down point
  InPlane.ProjectOntoPlane(NormalVector);
  InPlane += *ClosestPoint;

  DoLog(2) && (Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl);
  DoLog(2) && (Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl);
  DoLog(2) && (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 = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
    Line l = makeLineThrough(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(l);
    CrossDirection[i] = CrossPoint[i] - InPlane;
    CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
    const double s = CrossPoint[i].ScalarProduct(Direction)/Direction.NormSquared();
    DoLog(2) && (Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl);
    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
          CrossPoint[i] += (*endpoints[i%3]->node->node);  // make cross point absolute again
      DoLog(2) && (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) = 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) = 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) = (*endpoints[i]->node->node);
      }
    }
  }
  DoLog(1) && (Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl);
  return ShortestDistance;
}
;

/** Checks whether lines is any of the three boundary lines this triangle contains.
 * \param *line line to test
 * \return true - line is of the triangle, false - is not
 */
bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 3; i++)
    if (line == lines[i])
      return true;
  return false;
}
;

/** Checks whether point is any of the three endpoints this triangle contains.
 * \param *point point to test
 * \return true - point is of the triangle, false - is not
 */
bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
{
  Info FunctionInfo(__func__);
  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(const TesselPoint * const point) const
{
  Info FunctionInfo(__func__);
  for (int i = 0; i < 3; i++)
    if (point == endpoints[i]->node)
      return true;
  return false;
}
;

/** 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(const BoundaryPointSet * const Points[3]) const
{
  Info FunctionInfo(__func__);
  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking " << Points[0] << "," << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl);
  return (((endpoints[0] == Points[0]) || (endpoints[0] == Points[1]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2])

  ));
}
;

/** 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(const BoundaryTriangleSet * const T) const
{
  Info FunctionInfo(__func__);
  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
 * \return pointer third endpoint or NULL if line does not belong to triangle.
 */
class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
{
  Info FunctionInfo(__func__);
  // sanity check
  if (!ContainsBoundaryLine(line))
    return NULL;
  for (int i = 0; i < 3; i++)
    if (!line->ContainsBoundaryPoint(endpoints[i]))
      return endpoints[i];
  // actually, that' impossible :)
  return NULL;
}
;

/** Calculates the center point of the triangle.
 * Is third of the sum of all endpoints.
 * \param *center central point on return.
 */
void BoundaryTriangleSet::GetCenter(Vector * const center) const
{
  Info FunctionInfo(__func__);
  center->Zero();
  for (int i = 0; i < 3; i++)
    (*center) += (*endpoints[i]->node->node);
  center->Scale(1. / 3.);
  DoLog(1) && (Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl);
}

/**
 * gets the Plane defined by the three triangle Basepoints
 */
Plane BoundaryTriangleSet::getPlane() const{
  ASSERT(endpoints[0] && endpoints[1] && endpoints[2], "Triangle not fully defined");

  return Plane(*endpoints[0]->node->node,
               *endpoints[1]->node->node,
               *endpoints[2]->node->node);
}

Vector BoundaryTriangleSet::getEndpoint(int i) const{
  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");

  return *endpoints[i]->node->node;
}

string BoundaryTriangleSet::getEndpointName(int i) const{
  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");

  return endpoints[i]->node->getName();
}

/** output operator for BoundaryTriangleSet.
 * \param &ost output stream
 * \param &a boundary triangle
 */
ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a)
{
  ost << "[" << a.Nr << "|" << a.getEndpointName(0) << "," << a.getEndpointName(1) << "," << a.getEndpointName(2) << "]";
  //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
  //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
  return ost;
}
;

// ======================================== Polygons on Boundary =================================

/** Constructor for BoundaryPolygonSet.
 */
BoundaryPolygonSet::BoundaryPolygonSet() :
  Nr(-1)
{
  Info FunctionInfo(__func__);
}
;

/** Destructor of BoundaryPolygonSet.
 * Just clears endpoints.
 * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()
 */
BoundaryPolygonSet::~BoundaryPolygonSet()
{
  Info FunctionInfo(__func__);
  endpoints.clear();
  DoLog(1) && (Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl);
}
;

/** Calculates the normal vector for this triangle.
 * Is made unique by comparison with \a OtherVector to point in the other direction.
 * \param &OtherVector direction vector to make normal vector unique.
 * \return allocated vector in normal direction
 */
Vector * BoundaryPolygonSet::GetNormalVector(const Vector &OtherVector) const
{
  Info FunctionInfo(__func__);
  // get normal vector
  Vector TemporaryNormal;
  Vector *TotalNormal = new Vector;
  PointSet::const_iterator Runner[3];
  for (int i = 0; i < 3; i++) {
    Runner[i] = endpoints.begin();
    for (int j = 0; j < i; j++) { // go as much further
      Runner[i]++;
      if (Runner[i] == endpoints.end()) {
        DoeLog(0) && (eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl);
        performCriticalExit();
      }
    }
  }
  TotalNormal->Zero();
  int counter = 0;
  for (; Runner[2] != endpoints.end();) {
    TemporaryNormal = Plane(*((*Runner[0])->node->node),
                            *((*Runner[1])->node->node),
                            *((*Runner[2])->node->node)).getNormal();
    for (int i = 0; i < 3; i++) // increase each of them
      Runner[i]++;
    (*TotalNormal) += TemporaryNormal;
  }
  TotalNormal->Scale(1. / (double) counter);

  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
  if (TotalNormal->ScalarProduct(OtherVector) > 0.)
    TotalNormal->Scale(-1.);
  DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl);

  return TotalNormal;
}
;

/** Calculates the center point of the triangle.
 * Is third of the sum of all endpoints.
 * \param *center central point on return.
 */
void BoundaryPolygonSet::GetCenter(Vector * const center) const
{
  Info FunctionInfo(__func__);
  center->Zero();
  int counter = 0;
  for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
    (*center) += (*(*Runner)->node->node);
    counter++;
  }
  center->Scale(1. / (double) counter);
  DoLog(1) && (Log() << Verbose(1) << "Center is at " << *center << "." << endl);
}

/** Checks whether the polygons contains all three endpoints of the triangle.
 * \param *triangle triangle to test
 * \return true - triangle is contained polygon, false - is not
 */
bool BoundaryPolygonSet::ContainsBoundaryTriangle(const BoundaryTriangleSet * const triangle) const
{
  Info FunctionInfo(__func__);
  return ContainsPresentTupel(triangle->endpoints, 3);
}
;

/** Checks whether the polygons contains both endpoints of the line.
 * \param *line line to test
 * \return true - line is of the triangle, false - is not
 */
bool BoundaryPolygonSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
{
  Info FunctionInfo(__func__);
  return ContainsPresentTupel(line->endpoints, 2);
}
;

/** Checks whether point is any of the three endpoints this triangle contains.
 * \param *point point to test
 * \return true - point is of the triangle, false - is not
 */
bool BoundaryPolygonSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
{
  Info FunctionInfo(__func__);
  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
    DoLog(0) && (Log() << Verbose(0) << "Checking against " << **Runner << endl);
    if (point == (*Runner)) {
      DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
      return true;
    }
  }
  DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
  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 BoundaryPolygonSet::ContainsBoundaryPoint(const TesselPoint * const point) const
{
  Info FunctionInfo(__func__);
  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
    if (point == (*Runner)->node) {
      DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);
      return true;
    }
  DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);
  return false;
}
;

/** Checks whether given array of \a *Points coincide with polygons's endpoints.
 * \param **Points pointer to an array of BoundaryPointSet
 * \param dim dimension of array
 * \return true - set of points is contained in polygon, false - is not
 */
bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const
{
  Info FunctionInfo(__func__);
  int counter = 0;
  DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
  for (int i = 0; i < dim; i++) {
    DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl);
    if (ContainsBoundaryPoint(Points[i])) {
      counter++;
    }
  }

  if (counter == dim)
    return true;
  else
    return false;
}
;

/** Checks whether given PointList coincide with polygons's endpoints.
 * \param &endpoints PointList
 * \return true - set of points is contained in polygon, false - is not
 */
bool BoundaryPolygonSet::ContainsPresentTupel(const PointSet &endpoints) const
{
  Info FunctionInfo(__func__);
  size_t counter = 0;
  DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);
  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
    DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << **Runner << endl);
    if (ContainsBoundaryPoint(*Runner))
      counter++;
  }

  if (counter == endpoints.size())
    return true;
  else
    return false;
}
;

/** Checks whether given set of \a *Points coincide with polygons's endpoints.
 * \param *P pointer to BoundaryPolygonSet
 * \return true - is the very triangle, false - is not
 */
bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const
{
  return ContainsPresentTupel((const PointSet) P->endpoints);
}
;

/** Gathers all the endpoints' triangles in a unique set.
 * \return set of all triangles
 */
TriangleSet * BoundaryPolygonSet::GetAllContainedTrianglesFromEndpoints() const
{
  Info FunctionInfo(__func__);
  pair<TriangleSet::iterator, bool> Tester;
  TriangleSet *triangles = new TriangleSet;

  for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)
    for (LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)
      for (TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) {
        //Log() << Verbose(0) << " Testing triangle " << *(Sprinter->second) << endl;
        if (ContainsBoundaryTriangle(Sprinter->second)) {
          Tester = triangles->insert(Sprinter->second);
          if (Tester.second)
            DoLog(0) && (Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl);
        }
      }

  DoLog(1) && (Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl);
  return triangles;
}
;

/** Fills the endpoints of this polygon from the triangles attached to \a *line.
 * \param *line lines with triangles attached
 * \return true - polygon contains endpoints, false - line was NULL
 */
bool BoundaryPolygonSet::FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line)
{
  Info FunctionInfo(__func__);
  pair<PointSet::iterator, bool> Tester;
  if (line == NULL)
    return false;
  DoLog(1) && (Log() << Verbose(1) << "Filling polygon from line " << *line << endl);
  for (TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {
    for (int i = 0; i < 3; i++) {
      Tester = endpoints.insert((Runner->second)->endpoints[i]);
      if (Tester.second)
        DoLog(1) && (Log() << Verbose(1) << "  Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl);
    }
  }

  return true;
}
;

/** output operator for BoundaryPolygonSet.
 * \param &ost output stream
 * \param &a boundary polygon
 */
ostream &operator <<(ostream &ost, const BoundaryPolygonSet &a)
{
  ost << "[" << a.Nr << "|";
  for (PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) {
    ost << (*Runner)->node->getName();
    Runner++;
    if (Runner != a.endpoints.end())
      ost << ",";
  }
  ost << "]";
  return ost;
}
;

// =========================================================== class TESSELPOINT ===========================================

/** Constructor of class TesselPoint.
 */
TesselPoint::TesselPoint()
{
  //Info FunctionInfo(__func__);
  node = NULL;
  nr = -1;
}
;

/** Destructor for class TesselPoint.
 */
TesselPoint::~TesselPoint()
{
  //Info FunctionInfo(__func__);
}
;

/** Prints LCNode to screen.
 */
ostream & operator <<(ostream &ost, const TesselPoint &a)
{
  ost << "[" << a.getName() << "|" << *a.node << "]";
  return ost;
}
;

/** Prints LCNode to screen.
 */
ostream & TesselPoint::operator <<(ostream &ost)
{
  Info FunctionInfo(__func__);
  ost << "[" << (nr) << "|" << this << "]";
  return ost;
}
;

// =========================================================== class POINTCLOUD ============================================

/** Constructor of class PointCloud.
 */
PointCloud::PointCloud()
{
  //Info FunctionInfo(__func__);
}
;

/** Destructor for class PointCloud.
 */
PointCloud::~PointCloud()
{
  //Info FunctionInfo(__func__);
}
;

// ============================ CandidateForTesselation =============================

/** Constructor of class CandidateForTesselation.
 */
CandidateForTesselation::CandidateForTesselation(BoundaryLineSet* line) :
  BaseLine(line), ThirdPoint(NULL), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
{
  Info FunctionInfo(__func__);
}
;

/** Constructor of class CandidateForTesselation.
 */
CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) :
  BaseLine(line), ThirdPoint(point), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI)
{
        Info FunctionInfo(__func__);
  OptCenter = OptCandidateCenter;
  OtherOptCenter = OtherOptCandidateCenter;
};


/** Destructor for class CandidateForTesselation.
 */
CandidateForTesselation::~CandidateForTesselation()
{
}
;

/** Checks validity of a given sphere of a candidate line.
 * Sphere must touch all candidates and the baseline endpoints and there must be no other atoms inside.
 * \param RADIUS radius of sphere
 * \param *LC LinkedCell structure with other atoms
 * \return true - sphere is valid, false - sphere contains other points
 */
bool CandidateForTesselation::CheckValidity(const double RADIUS, const LinkedCell *LC) const
{
  Info FunctionInfo(__func__);

  const double radiusSquared = RADIUS * RADIUS;
  list<const Vector *> VectorList;
  VectorList.push_back(&OptCenter);
  //VectorList.push_back(&OtherOptCenter);  // don't check the other (wrong) center

  if (!pointlist.empty())
    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
  else
    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl);
  // check baseline for OptCenter and OtherOptCenter being on sphere's surface
  for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
    for (int i = 0; i < 2; i++) {
      const double distance = fabs((*VRunner)->DistanceSquared(*BaseLine->endpoints[i]->node->node) - radiusSquared);
      if (distance > HULLEPSILON) {
        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
        return false;
      }
    }
  }

  // check Candidates for OptCenter and OtherOptCenter being on sphere's surface
  for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
    const TesselPoint *Walker = *Runner;
    for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
      const double distance = fabs((*VRunner)->DistanceSquared(*Walker->node) - radiusSquared);
      if (distance > HULLEPSILON) {
        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl);
        return false;
      } else {
        DoLog(1) && (Log() << Verbose(1) << "Candidate " << *Walker << " is inside by " << distance << "." << endl);
      }
    }
  }

  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
  bool flag = true;
  for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
    // get all points inside the sphere
    TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner));

    DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << OtherOptCenter << ":" << endl);
    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(OtherOptCenter) << "." << endl);

    // remove baseline's endpoints and candidates
    for (int i = 0; i < 2; i++) {
      DoLog(1) && (Log() << Verbose(1) << "INFO: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "." << endl);
      ListofPoints->remove(BaseLine->endpoints[i]->node);
    }
    for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) {
      DoLog(1) && (Log() << Verbose(1) << "INFO: removing candidate tesselpoint " << *(*Runner) << "." << endl);
      ListofPoints->remove(*Runner);
    }
    if (!ListofPoints->empty()) {
      DoeLog(1) && (eLog() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
      flag = false;
      DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl);
      for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
        DoeLog(1) && (eLog() << Verbose(1) << "  " << *(*Runner) << endl);
    }
    delete (ListofPoints);

    // check with animate_sphere.tcl VMD script
    if (ThirdPoint != NULL) {
      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
    } else {
      DoLog(1) && (Log() << Verbose(1) << "Check by: ... missing third point ..." << endl);
      DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter[0] << " " << OldCenter[1] << " " << OldCenter[2] << " " << (*VRunner)->at(0) << " " << (*VRunner)->at(1) << " " << (*VRunner)->at(2) << endl);
    }
  }
  return flag;
}
;

/** output operator for CandidateForTesselation.
 * \param &ost output stream
 * \param &a boundary line
 */
ostream & operator <<(ostream &ost, const CandidateForTesselation &a)
{
  ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->getName() << "," << a.BaseLine->endpoints[1]->node->getName() << "] with ";
  if (a.pointlist.empty())
    ost << "no candidate.";
  else {
    ost << "candidate";
    if (a.pointlist.size() != 1)
      ost << "s ";
    else
      ost << " ";
    for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++)
      ost << *(*Runner) << " ";
    ost << " at angle " << (a.ShortestAngle) << ".";
  }

  return ost;
}
;

// =========================================================== class TESSELATION ===========================================

/** Constructor of class Tesselation.
 */
Tesselation::Tesselation() :
  PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin())
{
  Info FunctionInfo(__func__);
}
;

/** Destructor of class Tesselation.
 * We have to free all points, lines and triangles.
 */
Tesselation::~Tesselation()
{
  Info FunctionInfo(__func__);
  DoLog(0) && (Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl);
  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
    if (runner->second != NULL) {
      delete (runner->second);
      runner->second = NULL;
    } else
      DoeLog(1) && (eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl);
  }
  DoLog(0) && (Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl);
}
;

/** PointCloud implementation of GetCenter
 * Uses PointsOnBoundary and STL stuff.
 */
Vector * Tesselation::GetCenter(ofstream *out) const
{
  Info FunctionInfo(__func__);
  Vector *Center = new Vector(0., 0., 0.);
  int num = 0;
  for (GoToFirst(); (!IsEnd()); GoToNext()) {
    (*Center) += (*GetPoint()->node);
    num++;
  }
  Center->Scale(1. / num);
  return Center;
}
;

/** PointCloud implementation of GoPoint
 * Uses PointsOnBoundary and STL stuff.
 */
TesselPoint * Tesselation::GetPoint() const
{
  Info FunctionInfo(__func__);
  return (InternalPointer->second->node);
}
;

/** PointCloud implementation of GoToNext.
 * Uses PointsOnBoundary and STL stuff.
 */
void Tesselation::GoToNext() const
{
  Info FunctionInfo(__func__);
  if (InternalPointer != PointsOnBoundary.end())
    InternalPointer++;
}
;

/** PointCloud implementation of GoToFirst.
 * Uses PointsOnBoundary and STL stuff.
 */
void Tesselation::GoToFirst() const
{
  Info FunctionInfo(__func__);
  InternalPointer = PointsOnBoundary.begin();
}
;

/** PointCloud implementation of IsEmpty.
 * Uses PointsOnBoundary and STL stuff.
 */
bool Tesselation::IsEmpty() const
{
  Info FunctionInfo(__func__);
  return (PointsOnBoundary.empty());
}
;

/** PointCloud implementation of IsLast.
 * Uses PointsOnBoundary and STL stuff.
 */
bool Tesselation::IsEnd() const
{
  Info FunctionInfo(__func__);
  return (InternalPointer == PointsOnBoundary.end());
}
;

/** Gueses first starting triangle of the convex envelope.
 * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
 * \param *out output stream for debugging
 * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
 */
void Tesselation::GuessStartingTriangle()
{
  Info FunctionInfo(__func__);
  // 4b. create a starting triangle
  // 4b1. create all distances
  DistanceMultiMap DistanceMMap;
  double distance, tmp;
  Vector PlaneVector, TrialVector;
  PointMap::iterator A, B, C; // three nodes of the first triangle
  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily

  // with A chosen, take each pair B,C and sort
  if (A != PointsOnBoundary.end()) {
    B = A;
    B++;
    for (; B != PointsOnBoundary.end(); B++) {
      C = B;
      C++;
      for (; C != PointsOnBoundary.end(); C++) {
        tmp = A->second->node->node->DistanceSquared(*B->second->node->node);
        distance = tmp * tmp;
        tmp = A->second->node->node->DistanceSquared(*C->second->node->node);
        distance += tmp * tmp;
        tmp = B->second->node->node->DistanceSquared(*C->second->node->node);
        distance += tmp * tmp;
        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
      }
    }
  }
  //    // listing distances
  //    Log() << Verbose(1) << "Listing DistanceMMap:";
  //    for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
  //      Log() << Verbose(0) << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
  //    }
  //    Log() << Verbose(0) << endl;
  // 4b2. pick three baselines forming a triangle
  // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
  DistanceMultiMap::iterator baseline = DistanceMMap.begin();
  for (; baseline != DistanceMMap.end(); baseline++) {
    // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
    // 2. next, we have to check whether all points reside on only one side of the triangle
    // 3. construct plane vector
    PlaneVector = Plane(*A->second->node->node,
                        *baseline->second.first->second->node->node,
                        *baseline->second.second->second->node->node).getNormal();
    DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl);
    // 4. loop over all points
    double sign = 0.;
    PointMap::iterator checker = PointsOnBoundary.begin();
    for (; checker != PointsOnBoundary.end(); checker++) {
      // (neglecting A,B,C)
      if ((checker == A) || (checker == baseline->second.first) || (checker == baseline->second.second))
        continue;
      // 4a. project onto plane vector
      TrialVector = (*checker->second->node->node);
      TrialVector.SubtractVector(*A->second->node->node);
      distance = TrialVector.ScalarProduct(PlaneVector);
      if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
        continue;
      DoLog(2) && (Log() << Verbose(2) << "Projection of " << checker->second->node->getName() << " yields distance of " << distance << "." << endl);
      tmp = distance / fabs(distance);
      // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
      if ((sign != 0) && (tmp != sign)) {
        // 4c. If so, break 4. loop and continue with next candidate in 1. loop
        DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->getName() << "," << baseline->second.first->second->node->getName() << "," << baseline->second.second->second->node->getName() << " leaves " << checker->second->node->getName() << " outside the convex hull." << endl);
        break;
      } else { // note the sign for later
        DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->getName() << "," << baseline->second.first->second->node->getName() << "," << baseline->second.second->second->node->getName() << " leave " << checker->second->node->getName() << " inside the convex hull." << endl);
        sign = tmp;
      }
      // 4d. Check whether the point is inside the triangle (check distance to each node
      tmp = checker->second->node->node->DistanceSquared(*A->second->node->node);
      int innerpoint = 0;
      if ((tmp < A->second->node->node->DistanceSquared(*baseline->second.first->second->node->node)) && (tmp < A->second->node->node->DistanceSquared(*baseline->second.second->second->node->node)))
        innerpoint++;
      tmp = checker->second->node->node->DistanceSquared(*baseline->second.first->second->node->node);
      if ((tmp < baseline->second.first->second->node->node->DistanceSquared(*A->second->node->node)) && (tmp < baseline->second.first->second->node->node->DistanceSquared(*baseline->second.second->second->node->node)))
        innerpoint++;
      tmp = checker->second->node->node->DistanceSquared(*baseline->second.second->second->node->node);
      if ((tmp < baseline->second.second->second->node->node->DistanceSquared(*baseline->second.first->second->node->node)) && (tmp < baseline->second.second->second->node->node->DistanceSquared(*A->second->node->node)))
        innerpoint++;
      // 4e. If so, break 4. loop and continue with next candidate in 1. loop
      if (innerpoint == 3)
        break;
    }
    // 5. come this far, all on same side? Then break 1. loop and construct triangle
    if (checker == PointsOnBoundary.end()) {
      DoLog(2) && (Log() << Verbose(2) << "Looks like we have a candidate!" << endl);
      break;
    }
  }
  if (baseline != DistanceMMap.end()) {
    BPS[0] = baseline->second.first->second;
    BPS[1] = baseline->second.second->second;
    BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    BPS[0] = A->second;
    BPS[1] = baseline->second.second->second;
    BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
    BPS[0] = baseline->second.first->second;
    BPS[1] = A->second;
    BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);

    // 4b3. insert created triangle
    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
    TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
    TrianglesOnBoundaryCount++;
    for (int i = 0; i < NDIM; i++) {
      LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
      LinesOnBoundaryCount++;
    }

    DoLog(1) && (Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl);
  } else {
    DoeLog(0) && (eLog() << Verbose(0) << "No starting triangle found." << endl);
  }
}
;

/** Tesselates the convex envelope of a cluster from a single starting triangle.
 * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
 * 2 triangles. Hence, we go through all current lines:
 * -# if the lines contains to only one triangle
 * -# We search all points in the boundary
 *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
 *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
 *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
 *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
 * \param *out output stream for debugging
 * \param *configuration for IsAngstroem
 * \param *cloud cluster of points
 */
void Tesselation::TesselateOnBoundary(const PointCloud * const cloud)
{
  Info FunctionInfo(__func__);
  bool flag;
  PointMap::iterator winner;
  class BoundaryPointSet *peak = NULL;
  double SmallestAngle, TempAngle;
  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *Center = NULL;
  LineMap::iterator LineChecker[2];

  Center = cloud->GetCenter();
  // create a first tesselation with the given BoundaryPoints
  do {
    flag = false;
    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
      if (baseline->second->triangles.size() == 1) {
        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
        SmallestAngle = M_PI;

        // get peak point with respect to this base line's only triangle
        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
        DoLog(0) && (Log() << Verbose(0) << "Current baseline is between " << *(baseline->second) << "." << endl);
        for (int i = 0; i < 3; i++)
          if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
            peak = BTS->endpoints[i];
        DoLog(1) && (Log() << Verbose(1) << " and has peak " << *peak << "." << endl);

        // prepare some auxiliary vectors
        Vector BaseLineCenter, BaseLine;
        BaseLineCenter = 0.5 * ((*baseline->second->endpoints[0]->node->node) +
                                (*baseline->second->endpoints[1]->node->node));
        BaseLine = (*baseline->second->endpoints[0]->node->node) - (*baseline->second->endpoints[1]->node->node);

        // offset to center of triangle
        CenterVector.Zero();
        for (int i = 0; i < 3; i++)
          CenterVector += BTS->getEndpoint(i);
        CenterVector.Scale(1. / 3.);
        DoLog(2) && (Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl);

        // normal vector of triangle
        NormalVector = (*Center) - CenterVector;
        BTS->GetNormalVector(NormalVector);
        NormalVector = BTS->NormalVector;
        DoLog(2) && (Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl);

        // vector in propagation direction (out of triangle)
        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
        PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
        TempVector = CenterVector - (*baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
        //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
        if (PropagationVector.ScalarProduct(TempVector) > 0) // make sure normal propagation vector points outward from baseline
          PropagationVector.Scale(-1.);
        DoLog(2) && (Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl);
        winner = PointsOnBoundary.end();

        // loop over all points and calculate angle between normal vector of new and present triangle
        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
          if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
            DoLog(1) && (Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl);

            // first check direction, so that triangles don't intersect
            VirtualNormalVector = (*target->second->node->node) - BaseLineCenter;
            VirtualNormalVector.ProjectOntoPlane(NormalVector);
            TempAngle = VirtualNormalVector.Angle(PropagationVector);
            DoLog(2) && (Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl);
            if (TempAngle > (M_PI / 2.)) { // no bends bigger than Pi/2 (90 degrees)
              DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl);
              continue;
            } else
              DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl);

            // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
            LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
            LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
            if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) {
              DoLog(2) && (Log() << Verbose(2) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl);
              continue;
            }
            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) {
              DoLog(2) && (Log() << Verbose(2) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl);
              continue;
            }

            // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
            if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
              DoLog(4) && (Log() << Verbose(4) << "Current target is peak!" << endl);
              continue;
            }

            // check for linear dependence
            TempVector = (*baseline->second->endpoints[0]->node->node) - (*target->second->node->node);
            helper = (*baseline->second->endpoints[1]->node->node) - (*target->second->node->node);
            helper.ProjectOntoPlane(TempVector);
            if (fabs(helper.NormSquared()) < MYEPSILON) {
              DoLog(2) && (Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl);
              continue;
            }

            // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
            flag = true;
            VirtualNormalVector = Plane(*(baseline->second->endpoints[0]->node->node),
                                        *(baseline->second->endpoints[1]->node->node),
                                        *(target->second->node->node)).getNormal();
            TempVector = (1./3.) * ((*baseline->second->endpoints[0]->node->node) +
                                    (*baseline->second->endpoints[1]->node->node) +
                                    (*target->second->node->node));
            TempVector -= (*Center);
            // make it always point outward
            if (VirtualNormalVector.ScalarProduct(TempVector) < 0)
              VirtualNormalVector.Scale(-1.);
            // calculate angle
            TempAngle = NormalVector.Angle(VirtualNormalVector);
            DoLog(2) && (Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl);
            if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
              SmallestAngle = TempAngle;
              winner = target;
              DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);
            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
              helper = (*target->second->node->node) - BaseLineCenter;
              helper.ProjectOntoPlane(BaseLine);
              // ...the one with the smaller angle is the better candidate
              TempVector = (*target->second->node->node) - BaseLineCenter;
              TempVector.ProjectOntoPlane(VirtualNormalVector);
              TempAngle = TempVector.Angle(helper);
              TempVector = (*winner->second->node->node) - BaseLineCenter;
              TempVector.ProjectOntoPlane(VirtualNormalVector);
              if (TempAngle < TempVector.Angle(helper)) {
                TempAngle = NormalVector.Angle(VirtualNormalVector);
                SmallestAngle = TempAngle;
                winner = target;
                DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl);
              } else
                DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl);
            } else
              DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);
          }
        } // end of loop over all boundary points

        // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
        if (winner != PointsOnBoundary.end()) {
          DoLog(0) && (Log() << Verbose(0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl);
          // create the lins of not yet present
          BLS[0] = baseline->second;
          // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
          LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first);
          LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first);
          if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create
            BPS[0] = baseline->second->endpoints[0];
            BPS[1] = winner->second;
            BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[1]));
            LinesOnBoundaryCount++;
          } else
            BLS[1] = LineChecker[0]->second;
          if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create
            BPS[0] = baseline->second->endpoints[1];
            BPS[1] = winner->second;
            BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[2]));
            LinesOnBoundaryCount++;
          } else
            BLS[2] = LineChecker[1]->second;
          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
          BTS->GetCenter(&helper);
          helper -= (*Center);
          helper *= -1;
          BTS->GetNormalVector(helper);
          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
          TrianglesOnBoundaryCount++;
        } else {
          DoeLog(2) && (eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl);
        }

        // 5d. If the set of lines is not yet empty, go to 5. and continue
      } else
        DoLog(0) && (Log() << Verbose(0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl);
  } while (flag);

  // exit
  delete (Center);
}
;

/** Inserts all points outside of the tesselated surface into it by adding new triangles.
 * \param *out output stream for debugging
 * \param *cloud cluster of points
 * \param *LC LinkedCell structure to find nearest point quickly
 * \return true - all straddling points insert, false - something went wrong
 */
bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC)
{
  Info FunctionInfo(__func__);
  Vector Intersection, Normal;
  TesselPoint *Walker = NULL;
  Vector *Center = cloud->GetCenter();
  TriangleList *triangles = NULL;
  bool AddFlag = false;
  LinkedCell *BoundaryPoints = NULL;

  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;
    }
    Walker = cloud->GetPoint();
    DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);
    // get the next triangle
    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
    BTS = triangles->front();
    if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
      DoLog(0) && (Log() << Verbose(0) << "No triangles found, probably a tesselation point itself." << endl);
      cloud->GoToNext();
      continue;
    } else {
    }
    DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);
    // get the intersection point
    if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) {
      DoLog(0) && (Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl);
      // we have the intersection, check whether in- or outside of boundary
      if ((Center->DistanceSquared(*Walker->node) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
        // inside, next!
        DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
      } else {
        // outside!
        DoLog(0) && (Log() << Verbose(0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl);
        class BoundaryLineSet *OldLines[3], *NewLines[3];
        class BoundaryPointSet *OldPoints[3], *NewPoint;
        // store the three old lines and old points
        for (int i = 0; i < 3; i++) {
          OldLines[i] = BTS->lines[i];
          OldPoints[i] = BTS->endpoints[i];
        }
        Normal = BTS->NormalVector;
        // add Walker to boundary points
        DoLog(0) && (Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl);
        AddFlag = true;
        if (AddBoundaryPoint(Walker, 0))
          NewPoint = BPS[0];
        else
          continue;
        // remove triangle
        DoLog(0) && (Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl);
        TrianglesOnBoundary.erase(BTS->Nr);
        delete (BTS);
        // create three new boundary lines
        for (int i = 0; i < 3; i++) {
          BPS[0] = NewPoint;
          BPS[1] = OldPoints[i];
          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
          DoLog(1) && (Log() << Verbose(1) << "Creating new line " << *NewLines[i] << "." << endl);
          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
          LinesOnBoundaryCount++;
        }
        // create three new triangle with new point
        for (int i = 0; i < 3; i++) { // find all baselines
          BLS[0] = OldLines[i];
          int n = 1;
          for (int j = 0; j < 3; j++) {
            if (NewLines[j]->IsConnectedTo(BLS[0])) {
              if (n > 2) {
                DoeLog(2) && (eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl);
                return false;
              } else
                BLS[n++] = NewLines[j];
            }
          }
          // create the triangle
          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
          Normal.Scale(-1.);
          BTS->GetNormalVector(Normal);
          Normal.Scale(-1.);
          DoLog(0) && (Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl);
          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
          TrianglesOnBoundaryCount++;
        }
      }
    } else { // something is wrong with FindClosestTriangleToPoint!
      DoeLog(1) && (eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl);
      return false;
    }
    cloud->GoToNext();
  }

  // exit
  delete (Center);
  return true;
}
;

/** Adds a point to the tesselation::PointsOnBoundary list.
 * \param *Walker point to add
 * \param n TesselStruct::BPS index to put pointer into
 * \return true - new point was added, false - point already present
 */
bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n)
{
  Info FunctionInfo(__func__);
  PointTestPair InsertUnique;
  BPS[n] = new class BoundaryPointSet(Walker);
  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[n]));
  if (InsertUnique.second) { // if new point was not present before, increase counter
    PointsOnBoundaryCount++;
    return true;
  } else {
    delete (BPS[n]);
    BPS[n] = InsertUnique.first->second;
    return false;
  }
}
;

/** Adds point to Tesselation::PointsOnBoundary if not yet present.
 * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
 * @param Candidate point to add
 * @param n index for this point in Tesselation::TPS array
 */
void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n)
{
  Info FunctionInfo(__func__);
  PointTestPair InsertUnique;
  TPS[n] = new class BoundaryPointSet(Candidate);
  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
  if (InsertUnique.second) { // if new point was not present before, increase counter
    PointsOnBoundaryCount++;
  } else {
    delete TPS[n];
    DoLog(0) && (Log() << Verbose(0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl);
    TPS[n] = (InsertUnique.first)->second;
  }
}
;

/** Sets point to a present Tesselation::PointsOnBoundary.
 * Tesselation::TPS is set to the existing one or NULL if not found.
 * @param Candidate point to set to
 * @param n index for this point in Tesselation::TPS array
 */
void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const
{
  Info FunctionInfo(__func__);
  PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr);
  if (FindPoint != PointsOnBoundary.end())
    TPS[n] = FindPoint->second;
  else
    TPS[n] = NULL;
}
;

/** Function tries to add line from current Points in BPS to BoundaryLineSet.
 * If successful it raises the line count and inserts the new line into the BLS,
 * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
 * @param *OptCenter desired OptCenter if there are more than one candidate line
 * @param *candidate third point of the triangle to be, for checking between multiple open line candidates
 * @param *a first endpoint
 * @param *b second endpoint
 * @param n index of Tesselation::BLS giving the line with both endpoints
 */
void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
{
  bool insertNewLine = true;
  LineMap::iterator FindLine = a->lines.find(b->node->nr);
  BoundaryLineSet *WinningLine = NULL;
  if (FindLine != a->lines.end()) {
    DoLog(1) && (Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl);

    pair<LineMap::iterator, LineMap::iterator> FindPair;
    FindPair = a->lines.equal_range(b->node->nr);

    for (FindLine = FindPair.first; (FindLine != FindPair.second) && (insertNewLine); FindLine++) {
      DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl);
      // If there is a line with less than two attached triangles, we don't need a new line.
      if (FindLine->second->triangles.size() == 1) {
        CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
        if (!Finder->second->pointlist.empty())
          DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl);
        else
          DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate." << endl);
        // get open line
        for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
          if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches
            DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << "." << endl);
            insertNewLine = false;
            WinningLine = FindLine->second;
            break;
          } else {
            DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << "." << endl);
          }
        }
      }
    }
  }

  if (insertNewLine) {
    AddNewTesselationTriangleLine(a, b, n);
  } else {
    AddExistingTesselationTriangleLine(WinningLine, n);
  }
}
;

/**
 * Adds lines from each of the current points in the BPS to BoundaryLineSet.
 * Raises the line count and inserts the new line into the BLS.
 *
 * @param *a first endpoint
 * @param *b second endpoint
 * @param n index of Tesselation::BLS giving the line with both endpoints
 */
void Tesselation::AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
{
  Info FunctionInfo(__func__);
  DoLog(0) && (Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl);
  BPS[0] = a;
  BPS[1] = b;
  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); // this also adds the line to the local maps
  // add line to global map
  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
  // increase counter
  LinesOnBoundaryCount++;
  // also add to open lines
  CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
  OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
}
;

/** Uses an existing line for a new triangle.
 * Sets Tesselation::BLS[\a n] and removes the lines from Tesselation::OpenLines.
 * \param *FindLine the line to add
 * \param n index of the line to set in Tesselation::BLS
 */
void Tesselation::AddExistingTesselationTriangleLine(class BoundaryLineSet *Line, int n)
{
  Info FunctionInfo(__func__);
  DoLog(0) && (Log() << Verbose(0) << "Using existing line " << *Line << endl);

  // set endpoints and line
  BPS[0] = Line->endpoints[0];
  BPS[1] = Line->endpoints[1];
  BLS[n] = Line;
  // remove existing line from OpenLines
  CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
  if (CandidateLine != OpenLines.end()) {
    DoLog(1) && (Log() << Verbose(1) << " Removing line from OpenLines." << endl);
    delete (CandidateLine->second);
    OpenLines.erase(CandidateLine);
  } else {
    DoeLog(1) && (eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
  }
}
;

/** Function adds triangle to global list.
 * Furthermore, the triangle receives the next free id and id counter \a TrianglesOnBoundaryCount is increased.
 */
void Tesselation::AddTesselationTriangle()
{
  Info FunctionInfo(__func__);
  DoLog(1) && (Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl);

  // add triangle to global map
  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(const int nr)
{
  Info FunctionInfo(__func__);
  DoLog(0) && (Log() << Verbose(0) << "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
}
;

/** Removes a triangle from the tesselation.
 * Removes itself from the TriangleMap's of its lines, calls for them RemoveTriangleLine() if they are no more connected.
 * Removes itself from memory.
 * \param *triangle to remove
 */
void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle)
{
  Info FunctionInfo(__func__);
  if (triangle == NULL)
    return;
  for (int i = 0; i < 3; i++) {
    if (triangle->lines[i] != NULL) {
      DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl);
      triangle->lines[i]->triangles.erase(triangle->Nr);
      if (triangle->lines[i]->triangles.empty()) {
        DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl);
        RemoveTesselationLine(triangle->lines[i]);
      } else {
        DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: " << endl);
        OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL));
        for (TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
          DoLog(0) && (Log() << Verbose(0) << "\t[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t");
        DoLog(0) && (Log() << Verbose(0) << endl);
        //        for (int j=0;j<2;j++) {
        //          Log() << Verbose(0) << "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++)
        //            Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
        //          Log() << Verbose(0) << endl;
        //        }
      }
      triangle->lines[i] = NULL; // free'd or not: disconnect
    } else
      DoeLog(1) && (eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl);
  }

  if (TrianglesOnBoundary.erase(triangle->Nr))
    DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl);
  delete (triangle);
}
;

/** Removes a line from the tesselation.
 * Removes itself from each endpoints' LineMap, then removes itself from global LinesOnBoundary list and free's the line.
 * \param *line line to remove
 */
void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line)
{
  Info FunctionInfo(__func__);
  int Numbers[2];

  if (line == NULL)
    return;
  // get other endpoint number for finding copies of same line
  if (line->endpoints[1] != NULL)
    Numbers[0] = line->endpoints[1]->Nr;
  else
    Numbers[0] = -1;
  if (line->endpoints[0] != NULL)
    Numbers[1] = line->endpoints[0]->Nr;
  else
    Numbers[1] = -1;

  for (int i = 0; i < 2; i++) {
    if (line->endpoints[i] != NULL) {
      if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
        pair<LineMap::iterator, LineMap::iterator> erasor = line->endpoints[i]->lines.equal_range(Numbers[i]);
        for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
          if ((*Runner).second == line) {
            DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);
            line->endpoints[i]->lines.erase(Runner);
            break;
          }
      } else { // there's just a single line left
        if (line->endpoints[i]->lines.erase(line->Nr))
          DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);
      }
      if (line->endpoints[i]->lines.empty()) {
        DoLog(0) && (Log() << Verbose(0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl);
        RemoveTesselationPoint(line->endpoints[i]);
      } else {
        DoLog(0) && (Log() << Verbose(0) << *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++)
          DoLog(0) && (Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t");
        DoLog(0) && (Log() << Verbose(0) << endl);
      }
      line->endpoints[i] = NULL; // free'd or not: disconnect
    } else
      DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl);
  }
  if (!line->triangles.empty())
    DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl);

  if (LinesOnBoundary.erase(line->Nr))
    DoLog(0) && (Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl);
  delete (line);
}
;

/** Removes a point from the tesselation.
 * Checks whether there are still lines connected, removes from global PointsOnBoundary list, then free's the point.
 * \note If a point should be removed, while keep the tesselated surface intact (i.e. closed), use RemovePointFromTesselatedSurface()
 * \param *point point to remove
 */
void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point)
{
  Info FunctionInfo(__func__);
  if (point == NULL)
    return;
  if (PointsOnBoundary.erase(point->Nr))
    DoLog(0) && (Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl);
  delete (point);
}
;

/** Checks validity of a given sphere of a candidate line.
 * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
 * We check CandidateForTesselation::OtherOptCenter
 * \param &CandidateLine contains other degenerated candidates which we have to subtract as well
 * \param RADIUS radius of sphere
 * \param *LC LinkedCell structure with other atoms
 * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
 */
bool Tesselation::CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) const
{
  Info FunctionInfo(__func__);

  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
  bool flag = true;

  DoLog(1) && (Log() << Verbose(1) << "Check by: draw sphere {" << CandidateLine.OtherOptCenter[0] << " " << CandidateLine.OtherOptCenter[1] << " " << CandidateLine.OtherOptCenter[2] << "} radius " << RADIUS << " resolution 30" << endl);
  // get all points inside the sphere
  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);

  DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
  for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);

  // remove triangles's endpoints
  for (int i = 0; i < 2; i++)
    ListofPoints->remove(CandidateLine.BaseLine->endpoints[i]->node);

  // remove other candidates
  for (TesselPointList::const_iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); ++Runner)
    ListofPoints->remove(*Runner);

  // check for other points
  if (!ListofPoints->empty()) {
    DoLog(1) && (Log() << Verbose(1) << "CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
    flag = false;
    DoLog(1) && (Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);
  }
  delete (ListofPoints);

  return flag;
}
;

/** 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 integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
 *                 triangles exist which is the maximum for three points
 */
int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const
{
  Info FunctionInfo(__func__);
  int adjacentTriangleCount = 0;
  class BoundaryPointSet *Points[3];

  // builds a triangle point set (Points) of the end points
  for (int i = 0; i < 3; i++) {
    PointMap::const_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::const_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;
            DoLog(1) && (Log() << Verbose(1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl);
            for (TriangleMap::const_iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
              if (FindTriangle->second->IsPresentTupel(Points)) {
                adjacentTriangleCount++;
              }
            }
            DoLog(1) && (Log() << Verbose(1) << "end." << endl);
          }
          // Only one of the triangle lines must be considered for the triangle count.
          //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
          //return adjacentTriangleCount;
        }
      }
    }
  }

  DoLog(0) && (Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << 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(TesselPoint *Candidates[3])
{
  Info FunctionInfo(__func__);
  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.
          //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
          //return adjacentTriangleCount;
        }
      }
    }
  }

  return triangle;
}
;

/** Finds the starting triangle for FindNonConvexBorder().
 * Looks at the outermost point per axis, then FindSecondPointForTesselation()
 * for the second and FindNextSuitablePointViaAngleOfSphere() for the third
 * point are called.
 * \param *out output stream for debugging
 * \param RADIUS radius of virtual rolling sphere
 * \param *LC LinkedCell structure with neighbouring TesselPoint's
 * \return true - a starting triangle has been created, false - no valid triple of points found
 */
bool Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC)
{
  Info FunctionInfo(__func__);
  int i = 0;
  TesselPoint* MaxPoint[NDIM];
  TesselPoint* Temporary;
  double maxCoordinate[NDIM];
  BoundaryLineSet *BaseLine = NULL;
  Vector helper;
  Vector Chord;
  Vector SearchDirection;
  Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
  Vector SphereCenter;
  Vector NormalVector;

  NormalVector.Zero();

  for (i = 0; i < 3; i++) {
    MaxPoint[i] = NULL;
    maxCoordinate[i] = -1;
  }

  // 1. searching topmost point with respect to each axis
  for (int i = 0; i < NDIM; i++) { // each axis
    LC->n[i] = LC->N[i] - 1; // current axis is topmost cell
    for (LC->n[(i + 1) % NDIM] = 0; LC->n[(i + 1) % NDIM] < LC->N[(i + 1) % NDIM]; LC->n[(i + 1) % NDIM]++)
      for (LC->n[(i + 2) % NDIM] = 0; LC->n[(i + 2) % NDIM] < LC->N[(i + 2) % NDIM]; LC->n[(i + 2) % NDIM]++) {
        const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
        //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
        if (List != NULL) {
          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            if ((*Runner)->node->at(i) > maxCoordinate[i]) {
              DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl);
              maxCoordinate[i] = (*Runner)->node->at(i);
              MaxPoint[i] = (*Runner);
            }
          }
        } else {
          DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
        }
      }
  }

  DoLog(1) && (Log() << Verbose(1) << "Found maximum coordinates: ");
  for (int i = 0; i < NDIM; i++)
    DoLog(0) && (Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t");
  DoLog(0) && (Log() << Verbose(0) << endl);

  BTS = NULL;
  for (int k = 0; k < NDIM; k++) {
    NormalVector.Zero();
    NormalVector[k] = 1.;
    BaseLine = new BoundaryLineSet();
    BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
    DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl);

    double ShortestAngle;
    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.

    Temporary = NULL;
    FindSecondPointForTesselation(BaseLine->endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
    if (Temporary == NULL) {
      // have we found a second point?
      delete BaseLine;
      continue;
    }
    BaseLine->endpoints[1] = new BoundaryPointSet(Temporary);

    // construct center of circle
    CircleCenter = 0.5 * ((*BaseLine->endpoints[0]->node->node) + (*BaseLine->endpoints[1]->node->node));

    // construct normal vector of circle
    CirclePlaneNormal = (*BaseLine->endpoints[0]->node->node) - (*BaseLine->endpoints[1]->node->node);

    double radius = CirclePlaneNormal.NormSquared();
    double CircleRadius = sqrt(RADIUS * RADIUS - radius / 4.);

    NormalVector.ProjectOntoPlane(CirclePlaneNormal);
    NormalVector.Normalize();
    ShortestAngle = 2. * M_PI; // This will indicate the quadrant.

    SphereCenter = (CircleRadius * NormalVector) +  CircleCenter;
    // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)

    // look in one direction of baseline for initial candidate
    SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal();  // whether we look "left" first or "right" first is not important ...

    // adding point 1 and point 2 and add the line between them
    DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl);
    DoLog(0) && (Log() << Verbose(0) << "Found second point is at " << *BaseLine->endpoints[1]->node << ".\n");

    //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
    CandidateForTesselation OptCandidates(BaseLine);
    FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
    DoLog(0) && (Log() << Verbose(0) << "List of third Points is:" << endl);
    for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) {
      DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl);
    }
    if (!OptCandidates.pointlist.empty()) {
      BTS = NULL;
      AddCandidatePolygon(OptCandidates, RADIUS, LC);
    } else {
      delete BaseLine;
      continue;
    }

    if (BTS != NULL) { // we have created one starting triangle
      delete BaseLine;
      break;
    } else {
      // remove all candidates from the list and then the list itself
      OptCandidates.pointlist.clear();
    }
    delete BaseLine;
  }

  return (BTS != NULL);
}
;

/** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates.
 * This is supposed to prevent early closing of the tesselation.
 * \param CandidateLine CandidateForTesselation with baseline and shortestangle , i.e. not \a *OptCandidate
 * \param *ThirdNode third point in triangle, not in BoundaryLineSet::endpoints
 * \param RADIUS radius of sphere
 * \param *LC LinkedCell structure
 * \return true - there is a better candidate (smaller angle than \a ShortestAngle), false - no better TesselPoint candidate found
 */
//bool Tesselation::HasOtherBaselineBetterCandidate(CandidateForTesselation &CandidateLine, const TesselPoint * const ThirdNode, double RADIUS, const LinkedCell * const LC) const
//{
//      Info FunctionInfo(__func__);
//  bool result = false;
//  Vector CircleCenter;
//  Vector CirclePlaneNormal;
//  Vector OldSphereCenter;
//  Vector SearchDirection;
//  Vector helper;
//  TesselPoint *OtherOptCandidate = NULL;
//  double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant.
//  double radius, CircleRadius;
//  BoundaryLineSet *Line = NULL;
//  BoundaryTriangleSet *T = NULL;
//
//  // check both other lines
//  PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr);
//  if (FindPoint != PointsOnBoundary.end()) {
//    for (int i=0;i<2;i++) {
//      LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr);
//      if (FindLine != (FindPoint->second)->lines.end()) {
//        Line = FindLine->second;
//        Log() << Verbose(0) << "Found line " << *Line << "." << endl;
//        if (Line->triangles.size() == 1) {
//          T = Line->triangles.begin()->second;
//          // construct center of circle
//          CircleCenter.CopyVector(Line->endpoints[0]->node->node);
//          CircleCenter.AddVector(Line->endpoints[1]->node->node);
//          CircleCenter.Scale(0.5);
//
//          // construct normal vector of circle
//          CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node);
//          CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node);
//
//          // calculate squared radius of circle
//          radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
//          if (radius/4. < RADIUS*RADIUS) {
//            CircleRadius = RADIUS*RADIUS - radius/4.;
//            CirclePlaneNormal.Normalize();
//            //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
//
//            // construct old center
//            GetCenterofCircumcircle(&OldSphereCenter, *T->endpoints[0]->node->node, *T->endpoints[1]->node->node, *T->endpoints[2]->node->node);
//            helper.CopyVector(&T->NormalVector);  // normal vector ensures that this is correct center of the two possible ones
//            radius = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
//            helper.Scale(sqrt(RADIUS*RADIUS - radius));
//            OldSphereCenter.AddVector(&helper);
//            OldSphereCenter.SubtractVector(&CircleCenter);
//            //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
//
//            // construct SearchDirection
//            SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal);
//            helper.CopyVector(Line->endpoints[0]->node->node);
//            helper.SubtractVector(ThirdNode->node);
//            if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
//              SearchDirection.Scale(-1.);
//            SearchDirection.ProjectOntoPlane(&OldSphereCenter);
//            SearchDirection.Normalize();
//            Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
//            if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
//              // rotated the wrong way!
//              DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
//            }
//
//            // add third point
//            FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, OptCandidates, ThirdNode, RADIUS, LC);
//            for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); ++it) {
//              if (((*it) == BaseRay->endpoints[0]->node) || ((*it) == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested
//                continue;
//              Log() << Verbose(0) << " Third point candidate is " << (*it)
//              << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
//              Log() << Verbose(0) << " Baseline is " << *BaseRay << endl;
//
//              // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
//              TesselPoint *PointCandidates[3];
//              PointCandidates[0] = (*it);
//              PointCandidates[1] = BaseRay->endpoints[0]->node;
//              PointCandidates[2] = BaseRay->endpoints[1]->node;
//              bool check=false;
//              int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates);
//              // If there is no triangle, add it regularly.
//              if (existentTrianglesCount == 0) {
//                SetTesselationPoint((*it), 0);
//                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
//                SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
//
//                if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
//                  OtherOptCandidate = (*it);
//                  check = true;
//                }
//              } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
//                SetTesselationPoint((*it), 0);
//                SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
//                SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
//
//                // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
//                // i.e. at least one of the three lines must be present with TriangleCount <= 1
//                if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) {
//                  OtherOptCandidate = (*it);
//                  check = true;
//                }
//              }
//
//              if (check) {
//                if (ShortestAngle > OtherShortestAngle) {
//                  Log() << Verbose(0) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl;
//                  result = true;
//                  break;
//                }
//              }
//            }
//            delete(OptCandidates);
//            if (result)
//              break;
//          } else {
//            Log() << Verbose(0) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl;
//          }
//        } else {
//          DoeLog(2) && (eLog()<< Verbose(2) << "Baseline is connected to two triangles already?" << endl);
//        }
//      } else {
//        Log() << Verbose(1) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl;
//      }
//    }
//  } else {
//    DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl);
//  }
//
//  return result;
//};

/** This function finds a triangle to a line, adjacent to an existing one.
 * @param out output stream for debugging
 * @param CandidateLine current cadndiate baseline to search from
 * @param T current triangle which \a Line is edge of
 * @param RADIUS radius of the rolling ball
 * @param N number of found triangles
 * @param *LC LinkedCell structure with neighbouring points
 */
bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
{
  Info FunctionInfo(__func__);
  Vector CircleCenter;
  Vector CirclePlaneNormal;
  Vector RelativeSphereCenter;
  Vector SearchDirection;
  Vector helper;
  BoundaryPointSet *ThirdPoint = NULL;
  LineMap::iterator testline;
  double radius, CircleRadius;

  for (int i = 0; i < 3; i++)
    if ((T.endpoints[i] != CandidateLine.BaseLine->endpoints[0]) && (T.endpoints[i] != CandidateLine.BaseLine->endpoints[1])) {
      ThirdPoint = T.endpoints[i];
      break;
    }
  DoLog(0) && (Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdPoint " << *ThirdPoint << " of triangle " << T << "." << endl);

  CandidateLine.T = &T;

  // construct center of circle
  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
                        (*CandidateLine.BaseLine->endpoints[1]->node->node));

  // construct normal vector of circle
  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
                      (*CandidateLine.BaseLine->endpoints[1]->node->node);

  // calculate squared radius of circle
  radius = CirclePlaneNormal.ScalarProduct(CirclePlaneNormal);
  if (radius / 4. < RADIUS * RADIUS) {
    // construct relative sphere center with now known CircleCenter
    RelativeSphereCenter = T.SphereCenter - CircleCenter;

    CircleRadius = RADIUS * RADIUS - radius / 4.;
    CirclePlaneNormal.Normalize();
    DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);

    DoLog(1) && (Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl);

    // construct SearchDirection and an "outward pointer"
    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
    helper = CircleCenter - (*ThirdPoint->node->node);
    if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
      SearchDirection.Scale(-1.);
    DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
    if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {
      // rotated the wrong way!
      DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
    }

    // add third point
    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdPoint, RADIUS, LC);

  } else {
    DoLog(0) && (Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl);
  }

  if (CandidateLine.pointlist.empty()) {
    DoeLog(2) && (eLog() << Verbose(2) << "Could not find a suitable candidate." << endl);
    return false;
  }
  DoLog(0) && (Log() << Verbose(0) << "Third Points are: " << endl);
  for (TesselPointList::iterator it = CandidateLine.pointlist.begin(); it != CandidateLine.pointlist.end(); ++it) {
    DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl);
  }

  return true;
}
;

/** Walks through Tesselation::OpenLines() and finds candidates for newly created ones.
 * \param *&LCList atoms in LinkedCell list
 * \param RADIUS radius of the virtual sphere
 * \return true - for all open lines without candidates so far, a candidate has been found,
 *         false - at least one open line without candidate still
 */
bool Tesselation::FindCandidatesforOpenLines(const double RADIUS, const LinkedCell *&LCList)
{
  bool TesselationFailFlag = true;
  CandidateForTesselation *baseline = NULL;
  BoundaryTriangleSet *T = NULL;

  for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) {
    baseline = Runner->second;
    if (baseline->pointlist.empty()) {
      assert((baseline->BaseLine->triangles.size() == 1) && ("Open line without exactly one attached triangle"));
      T = (((baseline->BaseLine->triangles.begin()))->second);
      DoLog(1) && (Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl);
      TesselationFailFlag = TesselationFailFlag && FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
    }
  }
  return TesselationFailFlag;
}
;

/** Adds the present line and candidate point from \a &CandidateLine to the Tesselation.
 * \param CandidateLine triangle to add
 * \param RADIUS Radius of sphere
 * \param *LC LinkedCell structure
 * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in
 * AddTesselationLine() in AddCandidateTriangle()
 */
void Tesselation::AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC)
{
  Info FunctionInfo(__func__);
  Vector Center;
  TesselPoint * const TurningPoint = CandidateLine.BaseLine->endpoints[0]->node;
  TesselPointList::iterator Runner;
  TesselPointList::iterator Sprinter;

  // fill the set of neighbours
  TesselPointSet 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);

  DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl);
  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
    DoLog(0) && (Log() << Verbose(0) << " " << **TesselRunner << endl);

  // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
  Runner = connectedClosestPoints->begin();
  Sprinter = Runner;
  Sprinter++;
  while (Sprinter != connectedClosestPoints->end()) {
    DoLog(0) && (Log() << Verbose(0) << "Current Runner is " << *(*Runner) << " and sprinter is " << *(*Sprinter) << "." << endl);

    AddTesselationPoint(TurningPoint, 0);
    AddTesselationPoint(*Runner, 1);
    AddTesselationPoint(*Sprinter, 2);

    AddCandidateTriangle(CandidateLine, Opt);

    Runner = Sprinter;
    Sprinter++;
    if (Sprinter != connectedClosestPoints->end()) {
      // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
      FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OptCenter); // Assume BTS contains last triangle
      DoLog(0) && (Log() << Verbose(0) << " There are still more triangles to add." << endl);
    }
    // pick candidates for other open lines as well
    FindCandidatesforOpenLines(RADIUS, LC);

    // check whether we add a degenerate or a normal triangle
    if (CheckDegeneracy(CandidateLine, RADIUS, LC)) {
      // add normal and degenerate triangles
      DoLog(1) && (Log() << Verbose(1) << "Triangle of endpoints " << *TPS[0] << "," << *TPS[1] << " and " << *TPS[2] << " is degenerated, adding both sides." << endl);
      AddCandidateTriangle(CandidateLine, OtherOpt);

      if (Sprinter != connectedClosestPoints->end()) {
        // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
        FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OtherOptCenter);
      }
      // pick candidates for other open lines as well
      FindCandidatesforOpenLines(RADIUS, LC);
    }
  }
  delete (connectedClosestPoints);
};

/** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate.
 * \param *Sprinter next candidate to which internal open lines are set
 * \param *OptCenter OptCenter for this candidate
 */
void Tesselation::FindDegeneratedCandidatesforOpenLines(TesselPoint * const Sprinter, const Vector * const OptCenter)
{
  Info FunctionInfo(__func__);

  pair<LineMap::iterator, LineMap::iterator> FindPair = TPS[0]->lines.equal_range(TPS[2]->node->nr);
  for (LineMap::const_iterator FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
    DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl);
    // If there is a line with less than two attached triangles, we don't need a new line.
    if (FindLine->second->triangles.size() == 1) {
      CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
      if (!Finder->second->pointlist.empty())
        DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl);
      else {
        DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter) << endl);
        Finder->second->T = BTS;  // is last triangle
        Finder->second->pointlist.push_back(Sprinter);
        Finder->second->ShortestAngle = 0.;
        Finder->second->OptCenter = *OptCenter;
      }
    }
  }
};

/** If a given \a *triangle is degenerated, this adds both sides.
 * i.e. the triangle with same BoundaryPointSet's but NormalVector in opposite direction.
 * Note that endpoints are stored in Tesselation::TPS
 * \param CandidateLine CanddiateForTesselation structure for the desired BoundaryLine
 * \param RADIUS radius of sphere
 * \param *LC pointer to LinkedCell structure
 */
void Tesselation::AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC)
{
  Info FunctionInfo(__func__);
  Vector Center;
  CandidateMap::const_iterator CandidateCheck = OpenLines.end();
  BoundaryTriangleSet *triangle = NULL;

  /// 1. Create or pick the lines for the first triangle
  DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for first triangle ..." << endl);
  for (int i = 0; i < 3; i++) {
    BLS[i] = NULL;
    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
    AddTesselationLine(&CandidateLine.OptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
  }

  /// 2. create the first triangle and NormalVector and so on
  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding first triangle with center at " << CandidateLine.OptCenter << " ..." << endl);
  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
  AddTesselationTriangle();

  // create normal vector
  BTS->GetCenter(&Center);
  Center -= CandidateLine.OptCenter;
  BTS->SphereCenter = CandidateLine.OptCenter;
  BTS->GetNormalVector(Center);
  // give some verbose output about the whole procedure
  if (CandidateLine.T != NULL)
    DoLog(0) && (Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
  else
    DoLog(0) && (Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
  triangle = BTS;

  /// 3. Gather candidates for each new line
  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding candidates to new lines ..." << endl);
  for (int i = 0; i < 3; i++) {
    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
    CandidateCheck = OpenLines.find(BLS[i]);
    if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
      if (CandidateCheck->second->T == NULL)
        CandidateCheck->second->T = triangle;
      FindNextSuitableTriangle(*(CandidateCheck->second), *CandidateCheck->second->T, RADIUS, LC);
    }
  }

  /// 4. Create or pick the lines for the second triangle
  DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for second triangle ..." << endl);
  for (int i = 0; i < 3; i++) {
    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
    AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
  }

  /// 5. create the second triangle and NormalVector and so on
  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangle with center at " << CandidateLine.OtherOptCenter << " ..." << endl);
  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
  AddTesselationTriangle();

  BTS->SphereCenter = CandidateLine.OtherOptCenter;
  // create normal vector in other direction
  BTS->GetNormalVector(triangle->NormalVector);
  BTS->NormalVector.Scale(-1.);
  // give some verbose output about the whole procedure
  if (CandidateLine.T != NULL)
    DoLog(0) && (Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
  else
    DoLog(0) && (Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);

  /// 6. Adding triangle to new lines
  DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangles to new lines ..." << endl);
  for (int i = 0; i < 3; i++) {
    DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
    CandidateCheck = OpenLines.find(BLS[i]);
    if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
      if (CandidateCheck->second->T == NULL)
        CandidateCheck->second->T = BTS;
    }
  }
}
;

/** Adds a triangle to the Tesselation structure from three given TesselPoint's.
 * Note that endpoints are in Tesselation::TPS.
 * \param CandidateLine CandidateForTesselation structure contains other information
 * \param type which opt center to add (i.e. which side) and thus which NormalVector to take
 */
void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine, enum centers type)
{
  Info FunctionInfo(__func__);
  Vector Center;
  Vector *OptCenter = (type == Opt) ? &CandidateLine.OptCenter : &CandidateLine.OtherOptCenter;

  // add the lines
  AddTesselationLine(OptCenter, TPS[2], TPS[0], TPS[1], 0);
  AddTesselationLine(OptCenter, TPS[1], TPS[0], TPS[2], 1);
  AddTesselationLine(OptCenter, TPS[0], TPS[1], TPS[2], 2);

  // add the triangles
  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
  AddTesselationTriangle();

  // create normal vector
  BTS->GetCenter(&Center);
  Center.SubtractVector(*OptCenter);
  BTS->SphereCenter = *OptCenter;
  BTS->GetNormalVector(Center);

  // give some verbose output about the whole procedure
  if (CandidateLine.T != NULL)
    DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
  else
    DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
}
;

/** Checks whether the quadragon of the two triangles connect to \a *Base is convex.
 * We look whether the closest point on \a *Base with respect to the other baseline is outside
 * of the segment formed by both endpoints (concave) or not (convex).
 * \param *out output stream for debugging
 * \param *Base line to be flipped
 * \return NULL - convex, otherwise endpoint that makes it concave
 */
class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base)
{
  Info FunctionInfo(__func__);
  class BoundaryPointSet *Spot = NULL;
  class BoundaryLineSet *OtherBase;
  Vector *ClosestPoint;

  int m = 0;
  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
      if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
        BPS[m++] = runner->second->endpoints[j];
  OtherBase = new class BoundaryLineSet(BPS, -1);

  DoLog(1) && (Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl);
  DoLog(1) && (Log() << Verbose(1) << "INFO: Other base line is " << *OtherBase << "." << endl);

  // get the closest point on each line to the other line
  ClosestPoint = GetClosestPointBetweenLine(Base, OtherBase);

  // delete the temporary other base line
  delete (OtherBase);

  // get the distance vector from Base line to OtherBase line
  Vector DistanceToIntersection[2], BaseLine;
  double distance[2];
  BaseLine = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
  for (int i = 0; i < 2; i++) {
    DistanceToIntersection[i] = (*ClosestPoint) - (*Base->endpoints[i]->node->node);
    distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
  }
  delete (ClosestPoint);
  if ((distance[0] * distance[1]) > 0) { // have same sign?
    DoLog(1) && (Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl);
    if (distance[0] < distance[1]) {
      Spot = Base->endpoints[0];
    } else {
      Spot = Base->endpoints[1];
    }
    return Spot;
  } else { // different sign, i.e. we are in between
    DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl);
    return NULL;
  }

}
;

void Tesselation::PrintAllBoundaryPoints(ofstream *out) const
{
  Info FunctionInfo(__func__);
  // print all lines
  DoLog(0) && (Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl);
  for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin(); PointRunner != PointsOnBoundary.end(); PointRunner++)
    DoLog(0) && (Log() << Verbose(0) << *(PointRunner->second) << endl);
}
;

void Tesselation::PrintAllBoundaryLines(ofstream *out) const
{
  Info FunctionInfo(__func__);
  // print all lines
  DoLog(0) && (Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl);
  for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
    DoLog(0) && (Log() << Verbose(0) << *(LineRunner->second) << endl);
}
;

void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const
{
  Info FunctionInfo(__func__);
  // print all triangles
  DoLog(0) && (Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl);
  for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
    DoLog(0) && (Log() << Verbose(0) << *(TriangleRunner->second) << endl);
}
;

/** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher".
 * \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(class BoundaryLineSet *Base)
{
  Info FunctionInfo(__func__);
  class BoundaryLineSet *OtherBase;
  Vector *ClosestPoint[2];
  double volume;

  int m = 0;
  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
      if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
        BPS[m++] = runner->second->endpoints[j];
  OtherBase = new class BoundaryLineSet(BPS, -1);

  DoLog(0) && (Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl);
  DoLog(0) && (Log() << Verbose(0) << "INFO: Other base line is " << *OtherBase << "." << endl);

  // get the closest point on each line to the other line
  ClosestPoint[0] = GetClosestPointBetweenLine(Base, OtherBase);
  ClosestPoint[1] = GetClosestPointBetweenLine(OtherBase, Base);

  // get the distance vector from Base line to OtherBase line
  Vector Distance = (*ClosestPoint[1]) - (*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
  delete (ClosestPoint[0]);
  delete (ClosestPoint[1]);
  delete (OtherBase);

  if (Distance.NormSquared() < MYEPSILON) { // check for intersection
    DoLog(0) && (Log() << Verbose(0) << "REJECT: Both lines have an intersection: Nothing to do." << endl);
    return false;
  } else { // check for sign against BaseLineNormal
    Vector BaseLineNormal;
    BaseLineNormal.Zero();
    if (Base->triangles.size() < 2) {
      DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
      return 0.;
    }
    for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
    DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
      BaseLineNormal += (runner->second->NormalVector);
    }
    BaseLineNormal.Scale(1. / 2.);

    if (Distance.ScalarProduct(BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
      DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl);
      // calculate volume summand as a general tetraeder
      return volume;
    } else { // Base higher than OtherBase -> do nothing
      DoLog(0) && (Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl);
      return 0.;
    }
  }
}
;

/** For a given baseline and its two connected triangles, flips the baseline.
 * I.e. we create the new baseline between the other two endpoints of these four
 * endpoints and reconstruct the two triangles accordingly.
 * \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(class BoundaryLineSet *Base)
{
  Info FunctionInfo(__func__);
  class BoundaryLineSet *OldLines[4], *NewLine;
  class BoundaryPointSet *OldPoints[2];
  Vector BaseLineNormal;
  int OldTriangleNrs[2], OldBaseLineNr;
  int i, m;

  // calculate NormalVector for later use
  BaseLineNormal.Zero();
  if (Base->triangles.size() < 2) {
    DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
    return NULL;
  }
  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
    DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
    BaseLineNormal += (runner->second->NormalVector);
  }
  BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()

  // get the two triangles
  // gather four endpoints and four lines
  for (int j = 0; j < 4; j++)
    OldLines[j] = NULL;
  for (int j = 0; j < 2; j++)
    OldPoints[j] = NULL;
  i = 0;
  m = 0;
  DoLog(0) && (Log() << Verbose(0) << "The four old lines are: ");
  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
      if (runner->second->lines[j] != Base) { // pick not the central baseline
        OldLines[i++] = runner->second->lines[j];
        DoLog(0) && (Log() << Verbose(0) << *runner->second->lines[j] << "\t");
      }
  DoLog(0) && (Log() << Verbose(0) << endl);
  DoLog(0) && (Log() << Verbose(0) << "The two old points are: ");
  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
    for (int j = 0; j < 3; j++) // all of their endpoints and baselines
      if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints
        OldPoints[m++] = runner->second->endpoints[j];
        DoLog(0) && (Log() << Verbose(0) << *runner->second->endpoints[j] << "\t");
      }
  DoLog(0) && (Log() << Verbose(0) << endl);

  // check whether everything is in place to create new lines and triangles
  if (i < 4) {
    DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
    return NULL;
  }
  for (int j = 0; j < 4; j++)
    if (OldLines[j] == NULL) {
      DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
      return NULL;
    }
  for (int j = 0; j < 2; j++)
    if (OldPoints[j] == NULL) {
      DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl);
      return NULL;
    }

  // remove triangles and baseline removes itself
  DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting baseline " << *Base << " from global list." << endl);
  OldBaseLineNr = Base->Nr;
  m = 0;
  // first obtain all triangle to delete ... (otherwise we pull the carpet (Base) from under the for-loop's feet)
  list <BoundaryTriangleSet *> TrianglesOfBase;
  for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); ++runner)
    TrianglesOfBase.push_back(runner->second);
  // .. then delete each triangle (which deletes the line as well)
  for (list <BoundaryTriangleSet *>::iterator runner = TrianglesOfBase.begin(); !TrianglesOfBase.empty(); runner = TrianglesOfBase.begin()) {
    DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting triangle " << *(*runner) << "." << endl);
    OldTriangleNrs[m++] = (*runner)->Nr;
    RemoveTesselationTriangle((*runner));
    TrianglesOfBase.erase(runner);
  }

  // construct new baseline (with same number as old one)
  BPS[0] = OldPoints[0];
  BPS[1] = OldPoints[1];
  NewLine = new class BoundaryLineSet(BPS, OldBaseLineNr);
  LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
  DoLog(0) && (Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl);

  // construct new triangles with flipped baseline
  i = -1;
  if (OldLines[0]->IsConnectedTo(OldLines[2]))
    i = 2;
  if (OldLines[0]->IsConnectedTo(OldLines[3]))
    i = 3;
  if (i != -1) {
    BLS[0] = OldLines[0];
    BLS[1] = OldLines[i];
    BLS[2] = NewLine;
    BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
    BTS->GetNormalVector(BaseLineNormal);
    AddTesselationTriangle(OldTriangleNrs[0]);
    DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);

    BLS[0] = (i == 2 ? OldLines[3] : OldLines[2]);
    BLS[1] = OldLines[1];
    BLS[2] = NewLine;
    BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
    BTS->GetNormalVector(BaseLineNormal);
    AddTesselationTriangle(OldTriangleNrs[1]);
    DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);
  } else {
    DoeLog(0) && (eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl);
    return NULL;
  }

  return NewLine;
}
;

/** Finds the second point of starting triangle.
 * \param *a first node
 * \param Oben vector indicating the outside
 * \param OptCandidate reference to recommended candidate on return
 * \param Storage[3] array storing angles and other candidate information
 * \param RADIUS radius of virtual sphere
 * \param *LC LinkedCell structure with neighbouring points
 */
void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC)
{
  Info FunctionInfo(__func__);
  Vector AngleCheck;
  class TesselPoint* Candidate = NULL;
  double norm = -1.;
  double angle = 0.;
  int N[NDIM];
  int Nlower[NDIM];
  int Nupper[NDIM];

  if (LC->SetIndexToNode(a)) { // get cell for the starting point
    for (int i = 0; i < NDIM; i++) // store indices of this cell
      N[i] = LC->n[i];
  } else {
    DoeLog(1) && (eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl);
    return;
  }
  // then go through the current and all neighbouring cells and check the contained points for possible candidates
  for (int i = 0; i < NDIM; i++) {
    Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
    Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
  }
  DoLog(0) && (Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << 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 LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
        //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
        if (List != NULL) {
          for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            Candidate = (*Runner);
            // check if we only have one unique point yet ...
            if (a != Candidate) {
              // Calculate center of the circle with radius RADIUS through points a and Candidate
              Vector OrthogonalizedOben, aCandidate, Center;
              double distance, scaleFactor;

              OrthogonalizedOben = Oben;
              aCandidate = (*a->node) - (*Candidate->node);
              OrthogonalizedOben.ProjectOntoPlane(aCandidate);
              OrthogonalizedOben.Normalize();
              distance = 0.5 * aCandidate.Norm();
              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
              OrthogonalizedOben.Scale(scaleFactor);

              Center = 0.5 * ((*Candidate->node) + (*a->node));
              Center += OrthogonalizedOben;

              AngleCheck = Center - (*a->node);
              norm = aCandidate.Norm();
              // second point shall have smallest angle with respect to Oben vector
              if (norm < RADIUS * 2.) {
                angle = AngleCheck.Angle(Oben);
                if (angle < Storage[0]) {
                  //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
                  DoLog(1) && (Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n");
                  OptCandidate = Candidate;
                  Storage[0] = angle;
                  //Log() << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
                } else {
                  //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl;
                }
              } else {
                //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
              }
            } else {
              //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
            }
          }
        } else {
          DoLog(0) && (Log() << Verbose(0) << "Linked cell list is empty." << endl);
        }
      }
}
;

/** This recursive function finds a third point, to form a triangle with two given ones.
 * Note that this function is for the starting triangle.
 * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
 * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
 * the center of the sphere is still fixed up to a single parameter. The band of possible values
 * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
 * us the "null" on this circle, the new center of the candidate point will be some way along this
 * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
 * by the normal vector of the base triangle that always points outwards by construction.
 * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
 * We construct the normal vector that defines the plane this circle lies in, it is just in the
 * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
 * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
 * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
 * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
 * both.
 * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
 * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
 * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
 * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
 * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
 * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
 * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa FindStartingTriangle())
 * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
 * @param CandidateLine CandidateForTesselation with the current base line and list of candidates and ShortestAngle
 * @param ThirdPoint third point to avoid in search
 * @param RADIUS radius of sphere
 * @param *LC LinkedCell structure with neighbouring points
 */
void Tesselation::FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
{
  Info FunctionInfo(__func__);
  Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
  Vector SphereCenter;
  Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
  Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
  Vector NewNormalVector; // normal vector of the Candidate's triangle
  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
  Vector RelativeOldSphereCenter;
  Vector NewPlaneCenter;
  double CircleRadius; // radius of this circle
  double radius;
  double otherradius;
  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
  TesselPoint *Candidate = NULL;

  DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl);

  // copy old center
  CandidateLine.OldCenter = OldSphereCenter;
  CandidateLine.ThirdPoint = ThirdPoint;
  CandidateLine.pointlist.clear();

  // construct center of circle
  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
                        (*CandidateLine.BaseLine->endpoints[1]->node->node));

  // construct normal vector of circle
  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
                      (*CandidateLine.BaseLine->endpoints[1]->node->node);

  RelativeOldSphereCenter = OldSphereCenter - CircleCenter;

  // calculate squared radius TesselPoint *ThirdPoint,f circle
  radius = CirclePlaneNormal.NormSquared() / 4.;
  if (radius < RADIUS * RADIUS) {
    CircleRadius = RADIUS * RADIUS - radius;
    CirclePlaneNormal.Normalize();
    DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);

    // test whether old center is on the band's plane
    if (fabs(RelativeOldSphereCenter.ScalarProduct(CirclePlaneNormal)) > HULLEPSILON) {
      DoeLog(1) && (eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(CirclePlaneNormal)) << "!" << endl);
      RelativeOldSphereCenter.ProjectOntoPlane(CirclePlaneNormal);
    }
    radius = RelativeOldSphereCenter.NormSquared();
    if (fabs(radius - CircleRadius) < HULLEPSILON) {
      DoLog(1) && (Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl);

      // check SearchDirection
      DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
      if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
        DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);
      }

      // get cell for the starting point
      if (LC->SetIndexToVector(&CircleCenter)) {
        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;
      } else {
        DoeLog(1) && (eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl);
        return;
      }
      // then go through the current and all neighbouring cells and check the contained points for possible candidates
      //Log() << Verbose(1) << "LC Intervals:";
      for (int i = 0; i < NDIM; i++) {
        Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
        Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
        //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
      }
      //Log() << Verbose(0) << 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 LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
            //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
            if (List != NULL) {
              for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
                Candidate = (*Runner);

                // check for three unique points
                DoLog(2) && (Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl);
                if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node)) {

                  // find center on the plane
                  GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
                  DoLog(1) && (Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl);

                  try {
                    NewNormalVector = Plane(*(CandidateLine.BaseLine->endpoints[0]->node->node),
                                            *(CandidateLine.BaseLine->endpoints[1]->node->node),
                                            *(Candidate->node)).getNormal();
                    DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl);
                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(NewPlaneCenter);
                    DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
                    DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
                    DoLog(1) && (Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl);
                    if (radius < RADIUS * RADIUS) {
                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(NewPlaneCenter);
                      if (fabs(radius - otherradius) < HULLEPSILON) {
                        // construct both new centers
                        NewSphereCenter = NewPlaneCenter;
                        OtherNewSphereCenter= NewPlaneCenter;
                        helper = NewNormalVector;
                        helper.Scale(sqrt(RADIUS * RADIUS - radius));
                        DoLog(2) && (Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl);
                        NewSphereCenter += helper;
                        DoLog(2) && (Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl);
                        // OtherNewSphereCenter is created by the same vector just in the other direction
                        helper.Scale(-1.);
                        OtherNewSphereCenter += helper;
                        DoLog(2) && (Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl);
                        alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                        Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                        if ((ThirdPoint != NULL) && (Candidate == ThirdPoint->node)) { // in that case only the other circlecenter is valid
                          if (OldSphereCenter.DistanceSquared(NewSphereCenter) < OldSphereCenter.DistanceSquared(OtherNewSphereCenter))
                            alpha = Otheralpha;
                        } else
                          alpha = min(alpha, Otheralpha);
                        // if there is a better candidate, drop the current list and add the new candidate
                        // otherwise ignore the new candidate and keep the list
                        if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
                          if (fabs(alpha - Otheralpha) > MYEPSILON) {
                            CandidateLine.OptCenter = NewSphereCenter;
                            CandidateLine.OtherOptCenter = OtherNewSphereCenter;
                          } else {
                            CandidateLine.OptCenter = OtherNewSphereCenter;
                            CandidateLine.OtherOptCenter = NewSphereCenter;
                          }
                          // if there is an equal candidate, add it to the list without clearing the list
                          if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
                            CandidateLine.pointlist.push_back(Candidate);
                            DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl);
                          } else {
                            // remove all candidates from the list and then the list itself
                            CandidateLine.pointlist.clear();
                            CandidateLine.pointlist.push_back(Candidate);
                            DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl);
                          }
                          CandidateLine.ShortestAngle = alpha;
                          DoLog(0) && (Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl);
                        } else {
                          if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
                            DoLog(1) && (Log() << Verbose(1) << "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl);
                          } else {
                            DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl);
                          }
                        }
                      } else {
                        DoLog(1) && (Log() << Verbose(1) << "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius) << endl);
                      }
                    } else {
                      DoLog(1) && (Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl);
                    }
                  }
                  catch (LinearDependenceException &excp){
                    Log() << Verbose(1) << excp;
                    Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
                  }
                } else {
                  if (ThirdPoint != NULL) {
                    DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << "." << endl);
                  } else {
                    DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl);
                  }
                }
              }
            }
          }
    } else {
      DoeLog(1) && (eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
    }
  } else {
    if (ThirdPoint != NULL)
      DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!" << endl);
    else
      DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl);
  }

  DoLog(1) && (Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl);
  if (CandidateLine.pointlist.size() > 1) {
    CandidateLine.pointlist.unique();
    CandidateLine.pointlist.sort(); //SortCandidates);
  }

  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
    DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl);
    performCriticalExit();
  }
}
;

/** Finds the endpoint two lines are sharing.
 * \param *line1 first line
 * \param *line2 second line
 * \return point which is shared or NULL if none
 */
class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const
{
  Info FunctionInfo(__func__);
  const BoundaryLineSet * lines[2] = { line1, line2 };
  class BoundaryPointSet *node = NULL;
  PointMap OrderMap;
  PointTestPair OrderTest;
  for (int i = 0; i < 2; i++)
    // for both lines
    for (int j = 0; j < 2; j++) { // for both endpoints
      OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
      if (!OrderTest.second) { // if insertion fails, we have common endpoint
        node = OrderTest.first->second;
        DoLog(1) && (Log() << Verbose(1) << "Common endpoint of lines " << *line1 << " and " << *line2 << " is: " << *node << "." << endl);
        j = 2;
        i = 2;
        break;
      }
    }
  return node;
}
;

/** 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()) {
    DoeLog(1) && (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];
  DoLog(1) && (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 LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
        if (List != NULL) {
          for (LinkedCell::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));
              DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
            }
          }
        } else {
          DoeLog(1) && (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()) {
    DoeLog(1) && (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) {
    DoeLog(1) && (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
  DoLog(1) && (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 = 0.5 * ((*(LineRunner->second)->endpoints[0]->node->node) +
                      (*(LineRunner->second)->endpoints[1]->node->node));
      Center = (*x) - helper;
      BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
                 (*(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 = (*x) - (*(LineRunner->second)->endpoints[0]->node->node) - Center;
        const double lengthA = helper.ScalarProduct(BaseLine);
        helper = (*x) - (*(LineRunner->second)->endpoints[1]->node->node) - Center;
        const double lengthB = helper.ScalarProduct(BaseLine);
        if (lengthB * lengthA < 0) { // if have different sign
          ClosestLine = LineRunner->second;
          MinDistance = distance;
          DoLog(1) && (Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl);
        } else {
          DoLog(1) && (Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl);
        }
      } else {
        DoLog(1) && (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) {
    DoLog(0) && (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) {
    DoeLog(1) && (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
  DoLog(1) && (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 = (*(LineRunner->second)->endpoints[0]->node->node) -
                 (*(LineRunner->second)->endpoints[1]->node->node);
      const double lengthBase = BaseLine.NormSquared();

      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[0]->node->node);
      const double lengthEndA = BaseLineIntersection.NormSquared();

      BaseLineIntersection = (*x) - (*(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;
          DoLog(1) && (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);
          DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl);
        } else { // line is worse
          DoLog(1) && (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 = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
        Center = BaseLineIntersection;
        Center.ProjectOntoPlane(BaseLine);
        BaseLineIntersection -= Center;
        const double distance = BaseLineIntersection.NormSquared();
        if (Center.NormSquared() > BaseLine.NormSquared()) {
          DoeLog(0) && (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;
          DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl);
        } else {
          DoLog(2) && (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()) {
    DoLog(0) && (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;
}
;

/** Finds closest triangle to a point.
 * This basically just takes care of the degenerate case, which is not handled in FindClosestTrianglesToPoint().
 * \param *out output stream for debugging
 * \param *x Vector to look from
 * \param &distance contains found distance on return
 * \return list of BoundaryTriangleSet of nearest triangles or NULL.
 */
class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
{
  Info FunctionInfo(__func__);
  class BoundaryTriangleSet *result = NULL;
  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
  TriangleList candidates;
  Vector Center;
  Vector helper;

  if ((triangles == NULL) || (triangles->empty()))
    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 = (*x) - Center;
    const double Alignment = helper.Angle((*Runner)->NormalVector);
    if (Alignment < MinAlignment) {
      result = *Runner;
      MinAlignment = Alignment;
      DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl);
    } else {
      DoLog(1) && (Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << 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
{
  Info FunctionInfo(__func__);
  TriangleIntersectionList Intersections(&Point, this, LC);

  return Intersections.IsInside();
}
;

/** 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!
 *
 * @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__);
  Vector Center;
  Vector helper;
  Vector DistanceToCenter;
  Vector Intersection;
  double distance = 0.;

  if (triangle == NULL) {// is boundary point or only point in point cloud?
    DoLog(1) && (Log() << Verbose(1) << "No triangle given!" << endl);
    return -1.;
  } else {
    DoLog(1) && (Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl);
  }

  triangle->GetCenter(&Center);
  DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
  DistanceToCenter = Center - Point;
  DoLog(2) && (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 = Point + triangle->NormalVector; // points outside
    Center = Point - triangle->NormalVector; // points towards MolCenter
    DoLog(1) && (Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl);
    if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
      DoLog(1) && (Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl);
      return 0.;
    } else {
      DoLog(1) && (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);
    DoLog(1) && (Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl);

    // then check direction to boundary
    if (DistanceToCenter.ScalarProduct(triangle->NormalVector) > MYEPSILON) {
      DoLog(1) && (Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl);
      return -distance;
    } else {
      DoLog(1) && (Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl);
      return +distance;
    }
  }
}
;

/** Calculates minimum distance from \a&Point 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::GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const
{
  Info FunctionInfo(__func__);
  TriangleIntersectionList Intersections(&Point, this, LC);

  return Intersections.GetSmallestDistance();
}
;

/** Calculates minimum distance from \a&Point 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
 */
BoundaryTriangleSet * Tesselation::GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const
{
  Info FunctionInfo(__func__);
  TriangleIntersectionList Intersections(&Point, this, LC);

  return Intersections.GetClosestTriangle();
}
;

/** Gets all points connected to the provided point by triangulation lines.
 *
 * @param *Point of which get all connected points
 *
 * @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;
  class BoundaryPointSet *ReferencePoint = NULL;
  TesselPoint* current;
  bool takePoint = false;
  // find the respective boundary point
  PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
  if (PointRunner != PointsOnBoundary.end()) {
    ReferencePoint = PointRunner->second;
  } else {
    DoeLog(2) && (eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
    ReferencePoint = NULL;
  }

  // little trick so that we look just through lines connect to the BoundaryPoint
  // OR fall-back to look through all lines if there is no such BoundaryPoint
  const LineMap *Lines;
  ;
  if (ReferencePoint != NULL)
    Lines = &(ReferencePoint->lines);
  else
    Lines = &LinesOnBoundary;
  LineMap::const_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) {
      DoLog(1) && (Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl);
      connectedPoints->insert(current);
    }

    findLines++;
  }

  if (connectedPoints->empty()) { // if have not found any points
    DoeLog(1) && (eLog() << Verbose(1) << "We have not found any connected points to " << *Point << "." << endl);
    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.
 * 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 *out output stream for debugging
 * @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::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
{
  Info FunctionInfo(__func__);
  map<double, TesselPoint*> anglesOfPoints;
  TesselPointList *connectedCircle = new TesselPointList;
  Vector PlaneNormal;
  Vector AngleZero;
  Vector OrthogonalVector;
  Vector helper;
  const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL };
  TriangleList *triangles = NULL;

  if (SetOfNeighbours == NULL) {
    DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);
    delete (connectedCircle);
    return NULL;
  }

  // calculate central point
  triangles = FindTriangles(TrianglePoints);
  if ((triangles != NULL) && (!triangles->empty())) {
    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
      PlaneNormal += (*Runner)->NormalVector;
  } else {
    DoeLog(0) && (eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);
    performCriticalExit();
  }
  PlaneNormal.Scale(1.0 / triangles->size());
  DoLog(1) && (Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl);
  PlaneNormal.Normalize();

  // construct one orthogonal vector
  if (Reference != NULL) {
    AngleZero = (*Reference) - (*Point->node);
    AngleZero.ProjectOntoPlane(PlaneNormal);
  }
  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON)) {
    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
    AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
    AngleZero.ProjectOntoPlane(PlaneNormal);
    if (AngleZero.NormSquared() < MYEPSILON) {
      DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
      performCriticalExit();
    }
  }
  DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);
  if (AngleZero.NormSquared() > MYEPSILON)
    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
  else
    OrthogonalVector.MakeNormalTo(PlaneNormal);
  DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);

  // go through all connected points and calculate angle
  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
    helper = (*(*listRunner)->node) - (*Point->node);
    helper.ProjectOntoPlane(PlaneNormal);
    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
    DoLog(0) && (Log() << Verbose(0) << "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);
  }

  return connectedCircle;
}

/** 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) {
    DoeLog(2) && (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;
  }

  DoLog(1) && (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 = Plane(*((*TesselA)->node),
                   *((*TesselB)->node),
                   *((*TesselC)->node)).getNormal();
    DoLog(0) && (Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl);
    counter++;
    TesselA++;
    TesselB++;
    TesselC++;
    PlaneNormal += 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();
  DoLog(1) && (Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl);

  // construct one orthogonal vector
  if (Reference != NULL) {
    AngleZero = (*Reference) - (*Point->node);
    AngleZero.ProjectOntoPlane(PlaneNormal);
  }
  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);
    AngleZero = (*(*SetOfNeighbours->begin())->node) - (*Point->node);
    AngleZero.ProjectOntoPlane(PlaneNormal);
    if (AngleZero.NormSquared() < MYEPSILON) {
      DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
      performCriticalExit();
    }
  }
  DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);
  if (AngleZero.NormSquared() > MYEPSILON)
    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
  else
    OrthogonalVector.MakeNormalTo(PlaneNormal);
  DoLog(1) && (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 = (*(*listRunner)->node) - (*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;
    DoLog(0) && (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) {
      DoeLog(0) && (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.
 *
 * @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
 */
ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
{
  Info FunctionInfo(__func__);
  map<double, TesselPoint*> anglesOfPoints;
  list<TesselPointList *> *ListOfPaths = new list<TesselPointList *> ;
  TesselPointList *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::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
  if (PointRunner != PointsOnBoundary.end()) {
    ReferencePoint = PointRunner->second;
  } else {
    DoeLog(1) && (eLog() << Verbose(1) << "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()) {
        DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl);
      } else if (!LineRunner->second) {
        LineRunner->second = true;
        connectedPath = new TesselPointList;
        triangle = NULL;
        CurrentLine = runner->second;
        StartLine = CurrentLine;
        CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
        DoLog(1) && (Log() << Verbose(1) << "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl);
        do {
          // push current one
          DoLog(1) && (Log() << Verbose(1) << "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++) {
            DoLog(1) && (Log() << Verbose(1) << "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;
                  DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl);
                  break;
                } else {
                  DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl);
                  triangle = NULL;
                }
              } else {
                DoeLog(1) && (eLog() << Verbose(1) << "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];
              DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl);
              break;
            }
          }
          LineRunner = TouchedLine.find(CurrentLine);
          if (LineRunner == TouchedLine.end())
            DoeLog(1) && (eLog() << Verbose(1) << "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
        DoLog(1) && (Log() << Verbose(1) << "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 {
        DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl);
      }
    }
  } else {
    DoeLog(1) && (eLog() << Verbose(1) << "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
 */
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;
  int count = 0;
  TesselPointList::iterator CircleRunner;
  TesselPointList::iterator CircleStart;

  for (list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
    connectedPath = *ListRunner;

    DoLog(1) && (Log() << Verbose(1) << "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
    TesselPointList::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
        DoLog(1) && (Log() << Verbose(1) << count + 1 << ". closed path consists of: ");
        newPath = new TesselPointList;
        TesselPointList::iterator CircleSprinter = Marker;
        for (; CircleSprinter != CircleRunner; CircleSprinter++) {
          newPath->push_back(*CircleSprinter);
          DoLog(0) && (Log() << Verbose(0) << (**CircleSprinter) << " <-> ");
        }
        DoLog(0) && (Log() << Verbose(0) << ".." << endl);
        count++;
        Marker = CircleRunner;

        // add to list
        ListofClosedPaths->push_back(newPath);
      }
    }
  }
  DoLog(1) && (Log() << Verbose(1) << "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
 */
TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
{
  Info FunctionInfo(__func__);
  TriangleSet *connectedTriangles = new TriangleSet;

  if (Point == NULL) {
    DoeLog(1) && (eLog() << Verbose(1) << "Point given is NULL." << endl);
  } else {
    // go through its lines and insert all triangles
    for (LineMap::const_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).
 * \param *out output stream for debugging
 * \param *point point to be removed
 * \return volume added to the volume inside the tesselated surface by the removal
 */
double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point)
{
  class BoundaryLineSet *line = NULL;
  class BoundaryTriangleSet *triangle = NULL;
  Vector OldPoint, NormalVector;
  double volume = 0;
  int count = 0;

  if (point == NULL) {
    DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl);
    return 0.;
  } else
    DoLog(0) && (Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl);

  // copy old location for the volume
  OldPoint = (*point->node->node);

  // get list of connected points
  if (point->lines.empty()) {
    DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl);
    return 0.;
  }

  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
  TesselPointList *connectedPath = NULL;

  // gather all triangles
  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
    count += LineRunner->second->triangles.size();
  TriangleMap Candidates;
  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    line = LineRunner->second;
    for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
      triangle = TriangleRunner->second;
      Candidates.insert(TrianglePair(triangle->Nr, triangle));
    }
  }

  // remove all triangles
  count = 0;
  NormalVector.Zero();
  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
    DoLog(1) && (Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl);
    NormalVector -= Runner->second->NormalVector; // has to point inward
    RemoveTesselationTriangle(Runner->second);
    count++;
  }
  DoLog(1) && (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;
  double angle;
  double smallestangle;
  Vector Point, Reference, OrthogonalVector;
  if (count > 2) { // less than three triangles, then nothing will be created
    class TesselPoint *TriangleCandidates[3];
    count = 0;
    for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
      if (ListAdvance != ListOfClosedPaths->end())
        ListAdvance++;

      connectedPath = *ListRunner;
      // re-create all triangles by going through connected points list
      LineList NewLines;
      for (; !connectedPath->empty();) {
        // search middle node with widest angle to next neighbours
        EndNode = connectedPath->end();
        smallestangle = 0.;
        for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
          DoLog(1) && (Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);
          // construct vectors to next and previous neighbour
          StartNode = MiddleNode;
          if (StartNode == connectedPath->begin())
            StartNode = connectedPath->end();
          StartNode--;
          //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
          Point = (*(*StartNode)->node) - (*(*MiddleNode)->node);
          StartNode = MiddleNode;
          StartNode++;
          if (StartNode == connectedPath->end())
            StartNode = connectedPath->begin();
          //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
          Reference = (*(*StartNode)->node) - (*(*MiddleNode)->node);
          OrthogonalVector = (*(*MiddleNode)->node) - OldPoint;
          OrthogonalVector.MakeNormalTo(Reference);
          angle = GetAngle(Point, Reference, OrthogonalVector);
          //if (angle < M_PI)  // no wrong-sided triangles, please?
          if (fabs(angle - M_PI) < fabs(smallestangle - M_PI)) { // get straightest angle (i.e. construct those triangles with smallest area first)
            smallestangle = angle;
            EndNode = MiddleNode;
          }
        }
        MiddleNode = EndNode;
        if (MiddleNode == connectedPath->end()) {
          DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl);
          performCriticalExit();
        }
        StartNode = MiddleNode;
        if (StartNode == connectedPath->begin())
          StartNode = connectedPath->end();
        StartNode--;
        EndNode++;
        if (EndNode == connectedPath->end())
          EndNode = connectedPath->begin();
        DoLog(2) && (Log() << Verbose(2) << "INFO: StartNode is " << **StartNode << "." << endl);
        DoLog(2) && (Log() << Verbose(2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);
        DoLog(2) && (Log() << Verbose(2) << "INFO: EndNode is " << **EndNode << "." << endl);
        DoLog(1) && (Log() << Verbose(1) << "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << "." << endl);
        TriangleCandidates[0] = *StartNode;
        TriangleCandidates[1] = *MiddleNode;
        TriangleCandidates[2] = *EndNode;
        triangle = GetPresentTriangle(TriangleCandidates);
        if (triangle != NULL) {
          DoeLog(0) && (eLog() << Verbose(0) << "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;
        }
        DoLog(3) && (Log() << Verbose(3) << "Adding new triangle points." << endl);
        AddTesselationPoint(*StartNode, 0);
        AddTesselationPoint(*MiddleNode, 1);
        AddTesselationPoint(*EndNode, 2);
        DoLog(3) && (Log() << Verbose(3) << "Adding new triangle lines." << endl);
        AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
        AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
        NewLines.push_back(BLS[1]);
        AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
        BTS->GetNormalVector(NormalVector);
        AddTesselationTriangle();
        // calculate volume summand as a general tetraeder
        volume += CalculateVolumeofGeneralTetraeder(*TPS[0]->node->node, *TPS[1]->node->node, *TPS[2]->node->node, OldPoint);
        // advance number
        count++;

        // prepare nodes for next triangle
        StartNode = EndNode;
        DoLog(2) && (Log() << Verbose(2) << "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!
          DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl);
          performCriticalExit();
        } else {
          MiddleNode = StartNode;
          MiddleNode++;
          if (MiddleNode == connectedPath->end())
            MiddleNode = connectedPath->begin();
          EndNode = MiddleNode;
          EndNode++;
          if (EndNode == connectedPath->end())
            EndNode = connectedPath->begin();
        }
      }
      // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
      if (NewLines.size() > 1) {
        LineList::iterator Candidate;
        class BoundaryLineSet *OtherBase = NULL;
        double tmp, maxgain;
        do {
          maxgain = 0;
          for (LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
            tmp = PickFarthestofTwoBaselines(*Runner);
            if (maxgain < tmp) {
              maxgain = tmp;
              Candidate = Runner;
            }
          }
          if (maxgain != 0) {
            volume += maxgain;
            DoLog(1) && (Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl);
            OtherBase = FlipBaseline(*Candidate);
            NewLines.erase(Candidate);
            NewLines.push_back(OtherBase);
          }
        } while (maxgain != 0.);
      }

      ListOfClosedPaths->remove(connectedPath);
      delete (connectedPath);
    }
    DoLog(0) && (Log() << Verbose(0) << count << " triangles were created." << endl);
  } else {
    while (!ListOfClosedPaths->empty()) {
      ListRunner = ListOfClosedPaths->begin();
      connectedPath = *ListRunner;
      ListOfClosedPaths->remove(connectedPath);
      delete (connectedPath);
    }
    DoLog(0) && (Log() << Verbose(0) << "No need to create any triangles." << endl);
  }
  delete (ListOfClosedPaths);

  DoLog(0) && (Log() << Verbose(0) << "Removed volume is " << volume << "." << endl);

  return volume;
}
;

/**
 * Finds triangles belonging to the three provided points.
 *
 * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
 *
 * @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;
  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++;
      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);
                  }
                }
              }
              // 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:
      DoeLog(0) && (eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl);
      performCriticalExit();
      break;
  }

  return result;
}

struct BoundaryLineSetCompare
{
  bool operator()(const BoundaryLineSet * const a, const BoundaryLineSet * const b)
  {
    int lowerNra = -1;
    int lowerNrb = -1;

    if (a->endpoints[0] < a->endpoints[1])
      lowerNra = 0;
    else
      lowerNra = 1;

    if (b->endpoints[0] < b->endpoints[1])
      lowerNrb = 0;
    else
      lowerNrb = 1;

    if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])
      return true;
    else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
      return false;
    else { // both lower-numbered endpoints are the same ...
      if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2])
        return true;
      else if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2])
        return false;
    }
    return false;
  }
  ;
};

#define UniqueLines set < class BoundaryLineSet *, BoundaryLineSetCompare>

/**
 * 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
 */
IndexToIndex * Tesselation::FindAllDegeneratedLines()
{
  Info FunctionInfo(__func__);
  UniqueLines AllLines;
  IndexToIndex * DegeneratedLines = new IndexToIndex;

  // sanity check
  if (LinesOnBoundary.empty()) {
    DoeLog(2) && (eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.");
    return DegeneratedLines;
  }
  LineMap::iterator LineRunner1;
  pair<UniqueLines::iterator, bool> tester;
  for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
    tester = AllLines.insert(LineRunner1->second);
    if (!tester.second) { // found degenerated line
      DegeneratedLines->insert(pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr));
      DegeneratedLines->insert(pair<int, int> ((*tester.first)->Nr, LineRunner1->second->Nr));
    }
  }

  AllLines.clear();

  DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl);
  IndexToIndex::iterator it;
  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
    const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
    const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second);
    if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end())
      DoLog(0) && (Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl);
    else
      DoeLog(1) && (eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl);
  }

  return DegeneratedLines;
}

/**
 * Finds all degenerated triangles within the tesselation structure.
 *
 * @return map of keys of degenerated triangle pairs, each triangle occurs twice
 *         in the list, once as key and once as value
 */
IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
{
  Info FunctionInfo(__func__);
  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
  TriangleMap::iterator TriangleRunner1, TriangleRunner2;
  LineMap::iterator Liner;
  class BoundaryLineSet *line1 = NULL, *line2 = NULL;

  for (IndexToIndex::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));
        }
      }
    }
  }
  delete (DegeneratedLines);

  DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl);
  IndexToIndex::iterator it;
  for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
    DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);

  return DegeneratedTriangles;
}

/**
 * Purges degenerated triangles from the tesselation structure if they are not
 * necessary to keep a single point within the structure.
 */
void Tesselation::RemoveDegeneratedTriangles()
{
  Info FunctionInfo(__func__);
  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
  TriangleMap::iterator finder;
  BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
  int count = 0;

  for (IndexToIndex::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;

    bool trianglesShareLine = false;
    for (int i = 0; i < 3; ++i)
      for (int j = 0; j < 3; ++j)
        trianglesShareLine = trianglesShareLine || triangle->lines[i] == partnerTriangle->lines[j];

    if (trianglesShareLine && (triangle->endpoints[1]->LinesCount > 2) && (triangle->endpoints[2]->LinesCount > 2) && (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);
      DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl);
      RemoveTesselationTriangle(triangle);
      count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
      DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl);
      RemoveTesselationTriangle(partnerTriangle);
    } else {
      DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle << " and its partner " << *partnerTriangle << " because it is essential for at" << " least one of the endpoints to be kept in the tesselation structure." << endl);
    }
  }
  delete (DegeneratedTriangles);
  if (count > 0)
    LastTriangle = NULL;

  DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << 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(class TesselPoint *point, LinkedCell *LC)
{
  Info FunctionInfo(__func__);
  // find nearest boundary point
  class TesselPoint *BackupPoint = NULL;
  class TesselPoint *NearestPoint = FindClosestTesselPoint(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 {
    DoeLog(1) && (eLog() << Verbose(1) << "I cannot find the boundary point." << endl);
    return;
  }
  DoLog(0) && (Log() << Verbose(0) << "Nearest point on boundary is " << NearestPoint->getName() << "." << 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 = (*Runner->second->endpoints[0]->node->node) -
               (*Runner->second->endpoints[1]->node->node);
    CenterToPoint = 0.5 * ((*Runner->second->endpoints[0]->node->node) +
                           (*Runner->second->endpoints[1]->node->node));
    CenterToPoint -= (*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)
  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);
  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
  AddTesselationPoint(point, 2);
  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);
  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
  BTS->GetNormalVector(TempTriangle->NormalVector);
  BTS->NormalVector.Scale(-1.);
  DoLog(1) && (Log() << Verbose(1) << "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
  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);
  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
  AddTesselationPoint(point, 2);
  DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);
  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
  BTS->GetNormalVector(TempTriangle->NormalVector);
  DoLog(1) && (Log() << Verbose(1) << "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]) {
        DoeLog(0) && (eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl);
        performCriticalExit();
      }
      BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle));
      TempTriangle->lines[nr] = BTS->lines[i];
      break;
    }
  }
}
;

/** 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(const char *filename, const PointCloud * const cloud)
{
  Info FunctionInfo(__func__);
  ofstream *tempstream = NULL;
  string NameofTempFile;
  string NumberName;

  if (LastTriangle != NULL) {
    stringstream sstr;
    sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
    NumberName = sstr.str();
    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);
      DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");
      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
      WriteTecplotFile(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);
      DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");
      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
      WriteRaster3dFile(tempstream, this, cloud);
      IncludeSphereinRaster3D(tempstream, this, cloud);
      tempstream->close();
      tempstream->flush();
      delete (tempstream);
    }
  }
  if (DoTecplotOutput || DoRaster3DOutput)
    TriangleFilesWritten++;
}
;

struct BoundaryPolygonSetCompare
{
  bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const
  {
    if (s1->endpoints.size() < s2->endpoints.size())
      return true;
    else if (s1->endpoints.size() > s2->endpoints.size())
      return false;
    else { // equality of number of endpoints
      PointSet::const_iterator Walker1 = s1->endpoints.begin();
      PointSet::const_iterator Walker2 = s2->endpoints.begin();
      while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
        if ((*Walker1)->Nr < (*Walker2)->Nr)
          return true;
        else if ((*Walker1)->Nr > (*Walker2)->Nr)
          return false;
        Walker1++;
        Walker2++;
      }
      return false;
    }
  }
};

#define UniquePolygonSet set < BoundaryPolygonSet *, BoundaryPolygonSetCompare>

/** Finds all degenerated polygons and calls ReTesselateDegeneratedPolygon()/
 * \return number of polygons found
 */
int Tesselation::CorrectAllDegeneratedPolygons()
{
  Info FunctionInfo(__func__);
  /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
  set<BoundaryPointSet *> EndpointCandidateList;
  pair<set<BoundaryPointSet *>::iterator, bool> InsertionTester;
  pair<map<int, Vector *>::iterator, bool> TriangleInsertionTester;
  for (PointMap::const_iterator Runner = PointsOnBoundary.begin(); Runner != PointsOnBoundary.end(); Runner++) {
    DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Runner->second << "." << endl);
    map<int, Vector *> TriangleVectors;
    // gather all NormalVectors
    DoLog(1) && (Log() << Verbose(1) << "Gathering triangles ..." << endl);
    for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++)
      for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
        if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {
          TriangleInsertionTester = TriangleVectors.insert(pair<int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)));
          if (TriangleInsertionTester.second)
            DoLog(1) && (Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl);
        } else {
          DoLog(1) && (Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl);
        }
      }
    // check whether there are two that are parallel
    DoLog(1) && (Log() << Verbose(1) << "Finding two parallel triangles ..." << endl);
    for (map<int, Vector *>::iterator VectorWalker = TriangleVectors.begin(); VectorWalker != TriangleVectors.end(); VectorWalker++)
      for (map<int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++)
        if (VectorWalker != VectorRunner) { // skip equals
          const double SCP = VectorWalker->second->ScalarProduct(*VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles
          DoLog(1) && (Log() << Verbose(1) << "Checking " << *VectorWalker->second << " against " << *VectorRunner->second << ": " << SCP << endl);
          if (fabs(SCP + 1.) < ParallelEpsilon) {
            InsertionTester = EndpointCandidateList.insert((Runner->second));
            if (InsertionTester.second)
              DoLog(0) && (Log() << Verbose(0) << " Adding " << *Runner->second << " to endpoint candidate list." << endl);
            // and break out of both loops
            VectorWalker = TriangleVectors.end();
            VectorRunner = TriangleVectors.end();
            break;
          }
        }
  }
  delete DegeneratedTriangles;

  /// 3. Find connected endpoint candidates and put them into a polygon
  UniquePolygonSet ListofDegeneratedPolygons;
  BoundaryPointSet *Walker = NULL;
  BoundaryPointSet *OtherWalker = NULL;
  BoundaryPolygonSet *Current = NULL;
  stack<BoundaryPointSet*> ToCheckConnecteds;
  while (!EndpointCandidateList.empty()) {
    Walker = *(EndpointCandidateList.begin());
    if (Current == NULL) { // create a new polygon with current candidate
      DoLog(0) && (Log() << Verbose(0) << "Starting new polygon set at point " << *Walker << endl);
      Current = new BoundaryPolygonSet;
      Current->endpoints.insert(Walker);
      EndpointCandidateList.erase(Walker);
      ToCheckConnecteds.push(Walker);
    }

    // go through to-check stack
    while (!ToCheckConnecteds.empty()) {
      Walker = ToCheckConnecteds.top(); // fetch ...
      ToCheckConnecteds.pop(); // ... and remove
      for (LineMap::const_iterator LineWalker = Walker->lines.begin(); LineWalker != Walker->lines.end(); LineWalker++) {
        OtherWalker = (LineWalker->second)->GetOtherEndpoint(Walker);
        DoLog(1) && (Log() << Verbose(1) << "Checking " << *OtherWalker << endl);
        set<BoundaryPointSet *>::iterator Finder = EndpointCandidateList.find(OtherWalker);
        if (Finder != EndpointCandidateList.end()) { // found a connected partner
          DoLog(1) && (Log() << Verbose(1) << " Adding to polygon." << endl);
          Current->endpoints.insert(OtherWalker);
          EndpointCandidateList.erase(Finder); // remove from candidates
          ToCheckConnecteds.push(OtherWalker); // but check its partners too
        } else {
          DoLog(1) && (Log() << Verbose(1) << " is not connected to " << *Walker << endl);
        }
      }
    }

    DoLog(0) && (Log() << Verbose(0) << "Final polygon is " << *Current << endl);
    ListofDegeneratedPolygons.insert(Current);
    Current = NULL;
  }

  const int counter = ListofDegeneratedPolygons.size();

  DoLog(0) && (Log() << Verbose(0) << "The following " << counter << " degenerated polygons have been found: " << endl);
  for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++)
    DoLog(0) && (Log() << Verbose(0) << " " << **PolygonRunner << endl);

  /// 4. Go through all these degenerated polygons
  for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++) {
    stack<int> TriangleNrs;
    Vector NormalVector;
    /// 4a. Gather all triangles of this polygon
    TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();

    // check whether number is bigger than 2, otherwise it's just a simply degenerated one and nothing to do.
    if (T->size() == 2) {
      DoLog(1) && (Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl);
      delete (T);
      continue;
    }

    // check whether number is even
    // If this case occurs, we have to think about it!
    // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
    // connections to either polygon ...
    if (T->size() % 2 != 0) {
      DoeLog(0) && (eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl);
      performCriticalExit();
    }
    TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator
    /// 4a. Get NormalVector for one side (this is "front")
    NormalVector = (*TriangleWalker)->NormalVector;
    DoLog(1) && (Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl);
    TriangleWalker++;
    TriangleSet::iterator TriangleSprinter = TriangleWalker; // is the inner advanced iterator
    /// 4b. Remove all triangles whose NormalVector is in opposite direction (i.e. "back")
    BoundaryTriangleSet *triangle = NULL;
    while (TriangleSprinter != T->end()) {
      TriangleWalker = TriangleSprinter;
      triangle = *TriangleWalker;
      TriangleSprinter++;
      DoLog(1) && (Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl);
      if (triangle->NormalVector.ScalarProduct(NormalVector) < 0) { // if from other side, then delete and remove from list
        DoLog(1) && (Log() << Verbose(1) << " Removing ... " << endl);
        TriangleNrs.push(triangle->Nr);
        T->erase(TriangleWalker);
        RemoveTesselationTriangle(triangle);
      } else
        DoLog(1) && (Log() << Verbose(1) << " Keeping ... " << endl);
    }
    /// 4c. Copy all "front" triangles but with inverse NormalVector
    TriangleWalker = T->begin();
    while (TriangleWalker != T->end()) { // go through all front triangles
      DoLog(1) && (Log() << Verbose(1) << " Re-creating triangle " << **TriangleWalker << " with NormalVector " << (*TriangleWalker)->NormalVector << endl);
      for (int i = 0; i < 3; i++)
        AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i);
      AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
      AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
      AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
      if (TriangleNrs.empty())
        DoeLog(0) && (eLog() << Verbose(0) << "No more free triangle numbers!" << endl);
      BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ...
      AddTesselationTriangle(); // ... and add
      TriangleNrs.pop();
      BTS->NormalVector = -1 * (*TriangleWalker)->NormalVector;
      TriangleWalker++;
    }
    if (!TriangleNrs.empty()) {
      DoeLog(0) && (eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl);
    }
    delete (T); // remove the triangleset
  }
  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
  DoLog(0) && (Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl);
  IndexToIndex::iterator it;
  for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
    DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);
  delete (SimplyDegeneratedTriangles);
  /// 5. exit
  UniquePolygonSet::iterator PolygonRunner;
  while (!ListofDegeneratedPolygons.empty()) {
    PolygonRunner = ListofDegeneratedPolygons.begin();
    delete (*PolygonRunner);
    ListofDegeneratedPolygons.erase(PolygonRunner);
  }

  return counter;
}
;