Changes in src/tesselation.cpp [952f38:97b825]

File:

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

src/tesselation.cpp

@@ -13 +13 @@
 #include "Helpers/helpers.hpp"
 #include "Helpers/Info.hpp"
+#include "BoundaryPointSet.hpp"
+#include "BoundaryLineSet.hpp"
+#include "BoundaryTriangleSet.hpp"
+#include "BoundaryPolygonSet.hpp"
+#include "TesselPoint.hpp"
+#include "CandidateForTesselation.hpp"
+#include "PointCloud.hpp"
 #include "linkedcell.hpp"
 #include "Helpers/Log.hpp"
@@ -20 +27 @@
 #include "LinearAlgebra/Vector.hpp"
 #include "LinearAlgebra/Line.hpp"
-#include "vector_ops.hpp"
+#include "LinearAlgebra/vector_ops.hpp"
 #include "Helpers/Verbose.hpp"
 #include "LinearAlgebra/Plane.hpp"
@@ -28 +35 @@
 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__);
-  double angle = CalculateConvexity();
-  if (angle > -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;
-  }
-}
-
-
-/** Calculates the angle between two triangles with respect to their normal vector.
- * We sum the two angles of each height vector with respect to the center of the baseline.
- * \return angle > 0 then convex, if < 0 then concave
- */
-double BoundaryLineSet::CalculateConvexity() 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);
-  return (angle - M_PI);
-}
-
-/** 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::lines, 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;
-}
-;
-
-/** Returns other triangle of the line.
- * \param *point other endpoint
- * \return NULL - if triangle not contained in BoundaryLineSet::triangles, or pointer to BoundaryTriangleSet otherwise
- */
-class BoundaryTriangleSet *BoundaryLineSet::GetOtherTriangle(const BoundaryTriangleSet * const triangle) const
-{
-  Info FunctionInfo(__func__);
-  if (triangles.size() == 2) {
-    for (TriangleMap::const_iterator TriangleRunner = triangles.begin(); TriangleRunner != triangles.end(); ++TriangleRunner)
-      if (TriangleRunner->second != triangle)
-        return TriangleRunner->second;
-  }
-  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;
-}
-;
-
-/** Returns the baseline which does not contain the given boundary point \a *point.
- * \param *point endpoint which is neither endpoint of the desired line
- * \return pointer to desired third baseline
- */
-class BoundaryLineSet *BoundaryTriangleSet::GetThirdLine(const BoundaryPointSet * const point) const
-{
-  Info FunctionInfo(__func__);
-  // sanity check
-  if (!ContainsBoundaryPoint(point))
-    return NULL;
-  for (int i = 0; i < 3; i++)
-    if (!lines[i]->ContainsBoundaryPoint(point))
-      return lines[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 " << (*VRunner) << ":" << endl);
-    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(*(*VRunner)) << "." << 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) << " at distance " << setprecision(13) << (*Runner)->node->distance(*(*VRunner)) << setprecision(6) << "." << endl);
-
-      // check with animate_sphere.tcl VMD script
-      if (ThirdPoint != NULL) {
-        DoeLog(1) && (eLog() << 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 {
-        DoeLog(1) && (eLog() << Verbose(1) << "Check by: ... missing third point ..." << endl);
-        DoeLog(1) && (eLog() << 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);
-      }
-    }
-    delete (ListofPoints);
-
-  }
-  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 ===========================================
+const char *TecplotSuffix=".dat";
+const char *Raster3DSuffix=".r3d";
+const char *VRMLSUffix=".wrl";
+
+const double ParallelEpsilon=1e-3;
+const double Tesselation::HULLEPSILON = 1e-9;
 
 /** Constructor of class Tesselation.
  */
 Tesselation::Tesselation() :
-  PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin())
+  PointsOnBoundaryCount(0),
+  LinesOnBoundaryCount(0),
+  TrianglesOnBoundaryCount(0),
+  LastTriangle(NULL),
+  TriangleFilesWritten(0),
+  InternalPointer(PointsOnBoundary.begin())
 {
   Info FunctionInfo(__func__);
@@ -1254 +83 @@
   int num = 0;
   for (GoToFirst(); (!IsEnd()); GoToNext()) {
-    (*Center) += (*GetPoint()->node);
+    (*Center) += (GetPoint()->getPosition());
     num++;
   }
@@ -1337 +166 @@
       C++;
       for (; C != PointsOnBoundary.end(); C++) {
-        tmp = A->second->node->node->DistanceSquared(*B->second->node->node);
+        tmp = A->second->node->DistanceSquared(B->second->node->getPosition());
         distance = tmp * tmp;
-        tmp = A->second->node->node->DistanceSquared(*C->second->node->node);
+        tmp = A->second->node->DistanceSquared(C->second->node->getPosition());
         distance += tmp * tmp;
-        tmp = B->second->node->node->DistanceSquared(*C->second->node->node);
+        tmp = B->second->node->DistanceSquared(C->second->node->getPosition());
         distance += tmp * tmp;
         DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
@@ -1360 +189 @@
     // 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();
+    PlaneVector = Plane(A->second->node->getPosition(),
+                        baseline->second.first->second->node->getPosition(),
+                        baseline->second.second->second->node->getPosition()).getNormal();
     DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl);
     // 4. loop over all points
@@ -1372 +201 @@
         continue;
       // 4a. project onto plane vector
-      TrialVector = (*checker->second->node->node);
-      TrialVector.SubtractVector(*A->second->node->node);
+      TrialVector = (checker->second->node->getPosition() - A->second->node->getPosition());
       distance = TrialVector.ScalarProduct(PlaneVector);
       if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
@@ -1389 +217 @@
       }
       // 4d. Check whether the point is inside the triangle (check distance to each node
-      tmp = checker->second->node->node->DistanceSquared(*A->second->node->node);
+      tmp = checker->second->node->DistanceSquared(A->second->node->getPosition());
       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)))
+      if ((tmp < A->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < A->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
         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)))
+      tmp = checker->second->node->DistanceSquared(baseline->second.first->second->node->getPosition());
+      if ((tmp < baseline->second.first->second->node->DistanceSquared(A->second->node->getPosition())) && (tmp < baseline->second.first->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
         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)))
+      tmp = checker->second->node->DistanceSquared(baseline->second.second->second->node->getPosition());
+      if ((tmp < baseline->second.second->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < baseline->second.second->second->node->DistanceSquared(A->second->node->getPosition())))
         innerpoint++;
       // 4e. If so, break 4. loop and continue with next candidate in 1. loop
@@ -1478 +306 @@
         // 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);
+        BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                (baseline->second->endpoints[1]->node->getPosition()));
+        BaseLine = (baseline->second->endpoints[0]->node->getPosition()) - (baseline->second->endpoints[1]->node->getPosition());
 
         // offset to center of triangle
@@ -1498 +326 @@
         // 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!
+        TempVector = CenterVector - (baseline->second->endpoints[0]->node->getPosition()); // 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
@@ -1511 +339 @@
 
             // first check direction, so that triangles don't intersect
-            VirtualNormalVector = (*target->second->node->node) - BaseLineCenter;
+            VirtualNormalVector = (target->second->node->getPosition()) - BaseLineCenter;
             VirtualNormalVector.ProjectOntoPlane(NormalVector);
             TempAngle = VirtualNormalVector.Angle(PropagationVector);
@@ -1540 +368 @@
 
             // 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);
+            TempVector = (baseline->second->endpoints[0]->node->getPosition()) - (target->second->node->getPosition());
+            helper = (baseline->second->endpoints[1]->node->getPosition()) - (target->second->node->getPosition());
             helper.ProjectOntoPlane(TempVector);
             if (fabs(helper.NormSquared()) < MYEPSILON) {
@@ -1550 +378 @@
             // 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));
+            VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()),
+                                        (baseline->second->endpoints[1]->node->getPosition()),
+                                        (target->second->node->getPosition())).getNormal();
+            TempVector = (1./3.) * ((baseline->second->endpoints[0]->node->getPosition()) +
+                                    (baseline->second->endpoints[1]->node->getPosition()) +
+                                    (target->second->node->getPosition()));
             TempVector -= (*Center);
             // make it always point outward
@@ -1569 +397 @@
             } 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 = (target->second->node->getPosition()) - BaseLineCenter;
               helper.ProjectOntoPlane(BaseLine);
               // ...the one with the smaller angle is the better candidate
-              TempVector = (*target->second->node->node) - BaseLineCenter;
+              TempVector = (target->second->node->getPosition()) - BaseLineCenter;
               TempVector.ProjectOntoPlane(VirtualNormalVector);
               TempAngle = TempVector.Angle(helper);
-              TempVector = (*winner->second->node->node) - BaseLineCenter;
+              TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
               TempVector.ProjectOntoPlane(VirtualNormalVector);
               if (TempAngle < TempVector.Angle(helper)) {
@@ -1614 +442 @@
             BLS[2] = LineChecker[1]->second;
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          BTS->GetCenter(&helper);
+          BTS->GetCenter(helper);
           helper -= (*Center);
           helper *= -1;
@@ -1662 +490 @@
     DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);
     // get the next triangle
-    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->getPosition(), BoundaryPoints);
     if (triangles != NULL)
       BTS = triangles->front();
@@ -1676 +504 @@
     DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);
     // get the intersection point
-    if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) {
+    if (BTS->GetIntersectionInsideTriangle(*Center, Walker->getPosition(), 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) {
+      if ((Center->DistanceSquared(Walker->getPosition()) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
         // inside, next!
         DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
@@ -2086 +914 @@
   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);
+    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
 
   // remove triangles's endpoints
@@ -2102 +930 @@
     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);
+      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
   }
   delete (ListofPoints);
@@ -2257 +1085 @@
         if (List != NULL) {
           for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            if ((*Runner)->node->at(map[0]) > maxCoordinate[map[0]]) {
-              DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << map[0] << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl);
-              maxCoordinate[map[0]] = (*Runner)->node->at(map[0]);
+            if ((*Runner)->at(map[0]) > maxCoordinate[map[0]]) {
+              DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << map[0] << " node is " << *(*Runner) << " at " << (*Runner)->getPosition() << "." << endl);
+              maxCoordinate[map[0]] = (*Runner)->at(map[0]);
               MaxPoint[map[0]] = (*Runner);
             }
@@ -2295 +1123 @@
 
     // construct center of circle
-    CircleCenter = 0.5 * ((*BaseLine->endpoints[0]->node->node) + (*BaseLine->endpoints[1]->node->node));
+    CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
 
     // construct normal vector of circle
-    CirclePlaneNormal = (*BaseLine->endpoints[0]->node->node) - (*BaseLine->endpoints[1]->node->node);
+    CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
 
     double radius = CirclePlaneNormal.NormSquared();
@@ -2515 +1343 @@
 
   // construct center of circle
-  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
-                        (*CandidateLine.BaseLine->endpoints[1]->node->node));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
 
   // construct normal vector of circle
-  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
-                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
 
   // calculate squared radius of circle
@@ -2536 +1364 @@
     // construct SearchDirection and an "outward pointer"
     SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
-    helper = CircleCenter - (*ThirdPoint->node->node);
+    helper = CircleCenter - (ThirdPoint->node->getPosition());
     if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
@@ -2611 +1439 @@
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
-  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->getPosition());
 
   DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl);
@@ -2712 +1540 @@
 
   // create normal vector
-  BTS->GetCenter(&Center);
+  BTS->GetCenter(Center);
   Center -= CandidateLine.OptCenter;
   BTS->SphereCenter = CandidateLine.OptCenter;
@@ -2791 +1619 @@
 
   // create normal vector
-  BTS->GetCenter(&Center);
+  BTS->GetCenter(Center);
   Center.SubtractVector(*OptCenter);
   BTS->SphereCenter = *OptCenter;
@@ -2837 +1665 @@
   Vector DistanceToIntersection[2], BaseLine;
   double distance[2];
-  BaseLine = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
+  BaseLine = (Base->endpoints[1]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
   for (int i = 0; i < 2; i++) {
-    DistanceToIntersection[i] = (*ClosestPoint) - (*Base->endpoints[i]->node->node);
+    DistanceToIntersection[i] = (*ClosestPoint) - (Base->endpoints[i]->node->getPosition());
     distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
   }
@@ -2919 +1747 @@
 
   // calculate volume
-  volume = CalculateVolumeofGeneralTetraeder(*Base->endpoints[1]->node->node, *OtherBase->endpoints[0]->node->node, *OtherBase->endpoints[1]->node->node, *Base->endpoints[0]->node->node);
+  volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->getPosition(), OtherBase->endpoints[0]->node->getPosition(), OtherBase->endpoints[1]->node->getPosition(), Base->endpoints[0]->node->getPosition());
 
   // delete the temporary other base line and the closest points
@@ -3125 +1953 @@
 
               OrthogonalizedOben = Oben;
-              aCandidate = (*a->node) - (*Candidate->node);
+              aCandidate = (a->getPosition()) - (Candidate->getPosition());
               OrthogonalizedOben.ProjectOntoPlane(aCandidate);
               OrthogonalizedOben.Normalize();
@@ -3132 +1960 @@
               OrthogonalizedOben.Scale(scaleFactor);
 
-              Center = 0.5 * ((*Candidate->node) + (*a->node));
+              Center = 0.5 * ((Candidate->getPosition()) + (a->getPosition()));
               Center += OrthogonalizedOben;
 
-              AngleCheck = Center - (*a->node);
+              AngleCheck = Center - (a->getPosition());
               norm = aCandidate.Norm();
               // second point shall have smallest angle with respect to Oben vector
@@ -3220 +2048 @@
 
   // construct center of circle
-  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
-                        (*CandidateLine.BaseLine->endpoints[1]->node->node));
+  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
+                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
 
   // construct normal vector of circle
-  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
-                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
+  CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
+                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
 
   RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
@@ -3252 +2080 @@
 
       // get cell for the starting point
-      if (LC->SetIndexToVector(&CircleCenter)) {
+      if (LC->SetIndexToVector(CircleCenter)) {
         for (int i = 0; i < NDIM; i++) // store indices of this cell
           N[i] = LC->n[i];
@@ -3282 +2110 @@
 
                   // find center on the plane
-                  GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
+                  GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
                   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();
+                    NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()),
+                                            (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
+                                            (Candidate->getPosition())).getNormal();
                     DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl);
-                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(NewPlaneCenter);
+                    radius = CandidateLine.BaseLine->endpoints[0]->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);
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter);
                       if (fabs(radius - otherradius) < HULLEPSILON) {
                         // construct both new centers
@@ -3309 +2137 @@
                         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);
+                        alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
+                        Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
                         if ((ThirdPoint != NULL) && (Candidate == ThirdPoint->node)) { // in that case only the other circlecenter is valid
                           if (OldSphereCenter.DistanceSquared(NewSphereCenter) < OldSphereCenter.DistanceSquared(OtherNewSphereCenter))
@@ -3422 +2250 @@
  * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
  */
-DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector &x, const LinkedCell* LC) const
 {
   Info FunctionInfo(__func__);
@@ -3450 +2278 @@
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end()) {
-              points->insert(DistanceToPointPair(FindPoint->second->node->node->DistanceSquared(*x), FindPoint->second));
+              points->insert(DistanceToPointPair(FindPoint->second->node->DistanceSquared(x), FindPoint->second));
               DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
             }
@@ -3474 +2302 @@
  * \return closest BoundaryLineSet or NULL in degenerate case.
  */
-BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector &x, const LinkedCell* LC) const
 {
   Info FunctionInfo(__func__);
@@ -3485 +2313 @@
 
   // for each point, check its lines, remember closest
-  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl);
+  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << x << " ... " << endl);
   BoundaryLineSet *ClosestLine = NULL;
   double MinDistance = -1.;
@@ -3494 +2322 @@
     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);
+      helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) +
+                      ((LineRunner->second)->endpoints[1]->node->getPosition()));
+      Center = (x) - helper;
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
       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;
+        helper = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition()) - Center;
         const double lengthA = helper.ScalarProduct(BaseLine);
-        helper = (*x) - (*(LineRunner->second)->endpoints[1]->node->node) - Center;
+        helper = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition()) - Center;
         const double lengthB = helper.ScalarProduct(BaseLine);
         if (lengthB * lengthA < 0) { // if have different sign
@@ -3534 +2362 @@
  * \return BoundaryTriangleSet of nearest triangle or NULL.
  */
-TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector &x, const LinkedCell* LC) const
 {
   Info FunctionInfo(__func__);
@@ -3545 +2373 @@
 
   // for each point, check its lines, remember closest
-  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl);
+  DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << x << " ... " << endl);
   LineSet ClosestLines;
   double MinDistance = 1e+16;
@@ -3555 +2383 @@
     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);
+      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
+                 ((LineRunner->second)->endpoints[1]->node->getPosition());
       const double lengthBase = BaseLine.NormSquared();
 
-      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[0]->node->node);
+      BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition());
       const double lengthEndA = BaseLineIntersection.NormSquared();
 
-      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
+      BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
       const double lengthEndB = BaseLineIntersection.NormSquared();
 
@@ -3580 +2408 @@
       } else { // intersection is closer, calculate
         // calculate closest point on line to desired point
-        BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
+        BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
         Center = BaseLineIntersection;
         Center.ProjectOntoPlane(BaseLine);
@@ -3621 +2449 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector &x, const LinkedCell* LC) const
 {
   Info FunctionInfo(__func__);
@@ -3636 +2464 @@
   double MinAlignment = 2. * M_PI;
   for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
-    (*Runner)->GetCenter(&Center);
-    helper = (*x) - Center;
+    (*Runner)->GetCenter(Center);
+    helper = (x) - Center;
     const double Alignment = helper.Angle((*Runner)->NormalVector);
     if (Alignment < MinAlignment) {
@@ -3663 +2491 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point, this, LC);
+  TriangleIntersectionList Intersections(Point, this, LC);
 
   return Intersections.IsInside();
@@ -3703 +2531 @@
   }
 
-  triangle->GetCenter(&Center);
+  triangle->GetCenter(Center);
   DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
   DistanceToCenter = Center - Point;
@@ -3714 +2542 @@
     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)) {
+    if (triangle->GetIntersectionInsideTriangle(Center, DistanceToCenter, Intersection)) {
       DoLog(1) && (Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl);
       return 0.;
@@ -3723 +2551 @@
   } else {
     // calculate smallest distance
-    distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection);
+    distance = triangle->GetClosestPointInsideTriangle(Point, Intersection);
     DoLog(1) && (Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl);
 
@@ -3747 +2575 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point, this, LC);
+  TriangleIntersectionList Intersections(Point, this, LC);
 
   return Intersections.GetSmallestDistance();
@@ -3762 +2590 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(&Point, this, LC);
+  TriangleIntersectionList Intersections(Point, this, LC);
 
   return Intersections.GetClosestTriangle();
@@ -3839 +2667 @@
  * @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
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
 {
   Info FunctionInfo(__func__);
@@ -3871 +2699 @@
 
   // 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 = (Reference) - (Point->getPosition());
+  AngleZero.ProjectOntoPlane(PlaneNormal);
+  if ((AngleZero.NormSquared() < MYEPSILON)) {
+    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
+    AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
@@ -3893 +2719 @@
   // go through all connected points and calculate angle
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper = (*(*listRunner)->node) - (*Point->node);
+    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
     helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
@@ -3915 +2741 @@
  * @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
+ * @param *Reference Reference vector for zero angle or (0,0,0) 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
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
 {
   Info FunctionInfo(__func__);
@@ -3942 +2768 @@
   }
 
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl);
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << Reference << "." << endl);
   // calculate central point
   TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
@@ -3952 +2778 @@
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
-    helper = Plane(*((*TesselA)->node),
-                   *((*TesselB)->node),
-                   *((*TesselC)->node)).getNormal();
+    helper = Plane(((*TesselA)->getPosition()),
+                   ((*TesselB)->getPosition()),
+                   ((*TesselC)->getPosition())).getNormal();
     DoLog(0) && (Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl);
     counter++;
@@ -3974 +2800 @@
 
   // construct one orthogonal vector
-  if (Reference != NULL) {
-    AngleZero = (*Reference) - (*Point->node);
+  if (!Reference.IsZero()) {
+    AngleZero = (Reference) - (Point->getPosition());
     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);
+  if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON )) {
+    DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
+    AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
     AngleZero.ProjectOntoPlane(PlaneNormal);
     if (AngleZero.NormSquared() < MYEPSILON) {
@@ -3997 +2823 @@
   pair<map<double, TesselPoint*>::iterator, bool> InserterTest;
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper = (*(*listRunner)->node) - (*Point->node);
+    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
     helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
@@ -4242 +3068 @@
 
   // copy old location for the volume
-  OldPoint = (*point->node->node);
+  OldPoint = (point->node->getPosition());
 
   // get list of connected points
@@ -4305 +3131 @@
           StartNode--;
           //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
-          Point = (*(*StartNode)->node) - (*(*MiddleNode)->node);
+          Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
           StartNode = MiddleNode;
           StartNode++;
@@ -4311 +3137 @@
             StartNode = connectedPath->begin();
           //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
-          Reference = (*(*StartNode)->node) - (*(*MiddleNode)->node);
-          OrthogonalVector = (*(*MiddleNode)->node) - OldPoint;
+          Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+          OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
           OrthogonalVector.MakeNormalTo(Reference);
           angle = GetAngle(Point, Reference, OrthogonalVector);
@@ -4367 +3193 @@
         AddTesselationTriangle();
         // calculate volume summand as a general tetraeder
-        volume += CalculateVolumeofGeneralTetraeder(*TPS[0]->node->node, *TPS[1]->node->node, *TPS[2]->node->node, OldPoint);
+        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
         // advance number
         count++;
@@ -4752 +3578 @@
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->getPosition(), BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
@@ -4773 +3599 @@
   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);
+    BaseLine = (Runner->second->endpoints[0]->node->getPosition()) -
+               (Runner->second->endpoints[1]->node->getPosition());
+    CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) +
+                           (Runner->second->endpoints[1]->node->getPosition()));
+    CenterToPoint -= (point->getPosition());
     angle = CenterToPoint.Angle(BaseLine);
     if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {