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

File:

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

src/tesselation.cpp

@@ -11 +11 @@
 #include <iomanip>
 
-#include "BoundaryPointSet.hpp"
-#include "BoundaryLineSet.hpp"
-#include "BoundaryTriangleSet.hpp"
-#include "BoundaryPolygonSet.hpp"
-#include "TesselPoint.hpp"
-#include "CandidateForTesselation.hpp"
-#include "PointCloud.hpp"
-#include "helpers.hpp"
-#include "info.hpp"
+#include "Helpers/helpers.hpp"
+#include "Helpers/Info.hpp"
 #include "linkedcell.hpp"
-#include "log.hpp"
+#include "Helpers/Log.hpp"
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
 #include "triangleintersectionlist.hpp"
-#include "vector.hpp"
-#include "Line.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/Line.hpp"
 #include "vector_ops.hpp"
-#include "verbose.hpp"
-#include "Plane.hpp"
+#include "Helpers/Verbose.hpp"
+#include "LinearAlgebra/Plane.hpp"
 #include "Exceptions/LinearDependenceException.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__);
+  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 ===========================================
 
 /** Constructor of class Tesselation.
@@ -71 +1254 @@
   int num = 0;
   for (GoToFirst(); (!IsEnd()); GoToNext()) {
-    (*Center) += (GetPoint()->getPosition());
+    (*Center) += (*GetPoint()->node);
     num++;
   }
@@ -154 +1337 @@
       C++;
       for (; C != PointsOnBoundary.end(); C++) {
-        tmp = A->second->node->DistanceSquared(B->second->node->getPosition());
+        tmp = A->second->node->node->DistanceSquared(*B->second->node->node);
         distance = tmp * tmp;
-        tmp = A->second->node->DistanceSquared(C->second->node->getPosition());
+        tmp = A->second->node->node->DistanceSquared(*C->second->node->node);
         distance += tmp * tmp;
-        tmp = B->second->node->DistanceSquared(C->second->node->getPosition());
+        tmp = B->second->node->node->DistanceSquared(*C->second->node->node);
         distance += tmp * tmp;
         DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
@@ -177 +1360 @@
     // 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->getPosition(),
-                        baseline->second.first->second->node->getPosition(),
-                        baseline->second.second->second->node->getPosition()).getNormal();
+    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
@@ -189 +1372 @@
         continue;
       // 4a. project onto plane vector
-      TrialVector = (checker->second->node->getPosition() - A->second->node->getPosition());
+      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
@@ -205 +1389 @@
       }
       // 4d. Check whether the point is inside the triangle (check distance to each node
-      tmp = checker->second->node->DistanceSquared(A->second->node->getPosition());
+      tmp = checker->second->node->node->DistanceSquared(*A->second->node->node);
       int innerpoint = 0;
-      if ((tmp < A->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < A->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
+      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->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())))
+      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->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())))
+      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
@@ -294 +1478 @@
         // prepare some auxiliary vectors
         Vector BaseLineCenter, BaseLine;
-        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());
+        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
@@ -314 +1498 @@
         // 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->getPosition()); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+        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
@@ -327 +1511 @@
 
             // first check direction, so that triangles don't intersect
-            VirtualNormalVector = (target->second->node->getPosition()) - BaseLineCenter;
+            VirtualNormalVector = (*target->second->node->node) - BaseLineCenter;
             VirtualNormalVector.ProjectOntoPlane(NormalVector);
             TempAngle = VirtualNormalVector.Angle(PropagationVector);
@@ -356 +1540 @@
 
             // check for linear dependence
-            TempVector = (baseline->second->endpoints[0]->node->getPosition()) - (target->second->node->getPosition());
-            helper = (baseline->second->endpoints[1]->node->getPosition()) - (target->second->node->getPosition());
+            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) {
@@ -366 +1550 @@
             // 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->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()));
+            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
@@ -385 +1569 @@
             } 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->getPosition()) - BaseLineCenter;
+              helper = (*target->second->node->node) - BaseLineCenter;
               helper.ProjectOntoPlane(BaseLine);
               // ...the one with the smaller angle is the better candidate
-              TempVector = (target->second->node->getPosition()) - BaseLineCenter;
+              TempVector = (*target->second->node->node) - BaseLineCenter;
               TempVector.ProjectOntoPlane(VirtualNormalVector);
               TempAngle = TempVector.Angle(helper);
-              TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
+              TempVector = (*winner->second->node->node) - BaseLineCenter;
               TempVector.ProjectOntoPlane(VirtualNormalVector);
               if (TempAngle < TempVector.Angle(helper)) {
@@ -430 +1614 @@
             BLS[2] = LineChecker[1]->second;
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          BTS->GetCenter(helper);
+          BTS->GetCenter(&helper);
           helper -= (*Center);
           helper *= -1;
@@ -478 +1662 @@
     DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);
     // get the next triangle
-    triangles = FindClosestTrianglesToVector(Walker->getPosition(), BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     if (triangles != NULL)
       BTS = triangles->front();
@@ -492 +1676 @@
     DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);
     // get the intersection point
-    if (BTS->GetIntersectionInsideTriangle(*Center, Walker->getPosition(), Intersection)) {
+    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->getPosition()) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
+      if ((Center->DistanceSquared(*Walker->node) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
         // inside, next!
         DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
@@ -902 +2086 @@
   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)->distance(CandidateLine.OtherOptCenter) << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);
 
   // remove triangles's endpoints
@@ -918 +2102 @@
     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)->distance(CandidateLine.OtherOptCenter) << "." << endl);
+      DoLog(1) && (Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->distance(CandidateLine.OtherOptCenter) << "." << endl);
   }
   delete (ListofPoints);
@@ -1073 +2257 @@
         if (List != NULL) {
           for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            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]);
+            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]);
               MaxPoint[map[0]] = (*Runner);
             }
@@ -1111 +2295 @@
 
     // construct center of circle
-    CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
+    CircleCenter = 0.5 * ((*BaseLine->endpoints[0]->node->node) + (*BaseLine->endpoints[1]->node->node));
 
     // construct normal vector of circle
-    CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
+    CirclePlaneNormal = (*BaseLine->endpoints[0]->node->node) - (*BaseLine->endpoints[1]->node->node);
 
     double radius = CirclePlaneNormal.NormSquared();
@@ -1331 +2515 @@
 
   // construct center of circle
-  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
-                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
+  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->getPosition()) -
-                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
+  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
+                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // calculate squared radius of circle
@@ -1352 +2536 @@
     // construct SearchDirection and an "outward pointer"
     SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
-    helper = CircleCenter - (ThirdPoint->node->getPosition());
+    helper = CircleCenter - (*ThirdPoint->node->node);
     if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
@@ -1427 +2611 @@
   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->getPosition());
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
 
   DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl);
@@ -1528 +2712 @@
 
   // create normal vector
-  BTS->GetCenter(Center);
+  BTS->GetCenter(&Center);
   Center -= CandidateLine.OptCenter;
   BTS->SphereCenter = CandidateLine.OptCenter;
@@ -1607 +2791 @@
 
   // create normal vector
-  BTS->GetCenter(Center);
+  BTS->GetCenter(&Center);
   Center.SubtractVector(*OptCenter);
   BTS->SphereCenter = *OptCenter;
@@ -1653 +2837 @@
   Vector DistanceToIntersection[2], BaseLine;
   double distance[2];
-  BaseLine = (Base->endpoints[1]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
+  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->getPosition());
+    DistanceToIntersection[i] = (*ClosestPoint) - (*Base->endpoints[i]->node->node);
     distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
   }
@@ -1735 +2919 @@
 
   // calculate volume
-  volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->getPosition(), OtherBase->endpoints[0]->node->getPosition(), OtherBase->endpoints[1]->node->getPosition(), Base->endpoints[0]->node->getPosition());
+  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
@@ -1941 +3125 @@
 
               OrthogonalizedOben = Oben;
-              aCandidate = (a->getPosition()) - (Candidate->getPosition());
+              aCandidate = (*a->node) - (*Candidate->node);
               OrthogonalizedOben.ProjectOntoPlane(aCandidate);
               OrthogonalizedOben.Normalize();
@@ -1948 +3132 @@
               OrthogonalizedOben.Scale(scaleFactor);
 
-              Center = 0.5 * ((Candidate->getPosition()) + (a->getPosition()));
+              Center = 0.5 * ((*Candidate->node) + (*a->node));
               Center += OrthogonalizedOben;
 
-              AngleCheck = Center - (a->getPosition());
+              AngleCheck = Center - (*a->node);
               norm = aCandidate.Norm();
               // second point shall have smallest angle with respect to Oben vector
@@ -2036 +3220 @@
 
   // construct center of circle
-  CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
-                        (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
+  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->getPosition()) -
-                      (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
+  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
+                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
 
   RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
@@ -2068 +3252 @@
 
       // 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];
@@ -2098 +3282 @@
 
                   // find center on the plane
-                  GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
+                  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->getPosition()),
-                                            (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
-                                            (Candidate->getPosition())).getNormal();
+                    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->DistanceSquared(NewPlaneCenter);
+                    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->DistanceSquared(NewPlaneCenter);
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(NewPlaneCenter);
                       if (fabs(radius - otherradius) < HULLEPSILON) {
                         // construct both new centers
@@ -2238 +3422 @@
  * \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__);
@@ -2266 +3450 @@
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end()) {
-              points->insert(DistanceToPointPair(FindPoint->second->node->DistanceSquared(x), FindPoint->second));
+              points->insert(DistanceToPointPair(FindPoint->second->node->node->DistanceSquared(*x), FindPoint->second));
               DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
             }
@@ -2290 +3474 @@
  * \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__);
@@ -2301 +3485 @@
 
   // 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.;
@@ -2310 +3494 @@
     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->getPosition()) +
-                      ((LineRunner->second)->endpoints[1]->node->getPosition()));
-      Center = (x) - helper;
-      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
-                 ((LineRunner->second)->endpoints[1]->node->getPosition());
+      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->getPosition()) - Center;
+        helper = (*x) - (*(LineRunner->second)->endpoints[0]->node->node) - Center;
         const double lengthA = helper.ScalarProduct(BaseLine);
-        helper = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition()) - Center;
+        helper = (*x) - (*(LineRunner->second)->endpoints[1]->node->node) - Center;
         const double lengthB = helper.ScalarProduct(BaseLine);
         if (lengthB * lengthA < 0) { // if have different sign
@@ -2350 +3534 @@
  * \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__);
@@ -2361 +3545 @@
 
   // 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;
@@ -2371 +3555 @@
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
 
-      BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
-                 ((LineRunner->second)->endpoints[1]->node->getPosition());
+      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->getPosition());
+      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[0]->node->node);
       const double lengthEndA = BaseLineIntersection.NormSquared();
 
-      BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
+      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
       const double lengthEndB = BaseLineIntersection.NormSquared();
 
@@ -2396 +3580 @@
       } else { // intersection is closer, calculate
         // calculate closest point on line to desired point
-        BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
+        BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
         Center = BaseLineIntersection;
         Center.ProjectOntoPlane(BaseLine);
@@ -2437 +3621 @@
  * \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__);
@@ -2452 +3636 @@
   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) {
@@ -2479 +3663 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(Point, this, LC);
+  TriangleIntersectionList Intersections(&Point, this, LC);
 
   return Intersections.IsInside();
@@ -2519 +3703 @@
   }
 
-  triangle->GetCenter(Center);
+  triangle->GetCenter(&Center);
   DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
   DistanceToCenter = Center - Point;
@@ -2530 +3714 @@
     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.;
@@ -2539 +3723 @@
   } 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);
 
@@ -2563 +3747 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(Point, this, LC);
+  TriangleIntersectionList Intersections(&Point, this, LC);
 
   return Intersections.GetSmallestDistance();
@@ -2578 +3762 @@
 {
   Info FunctionInfo(__func__);
-  TriangleIntersectionList Intersections(Point, this, LC);
+  TriangleIntersectionList Intersections(&Point, this, LC);
 
   return Intersections.GetClosestTriangle();
@@ -2655 +3839 @@
  * @return list of the all points linked to the provided one
  */
-TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
   Info FunctionInfo(__func__);
@@ -2687 +3871 @@
 
   // construct one orthogonal vector
-  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());
+  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) {
@@ -2707 +3893 @@
   // go through all connected points and calculate angle
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
+    helper = (*(*listRunner)->node) - (*Point->node);
     helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
@@ -2729 +3915 @@
  * @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 (0,0,0) for no preference
+ * @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 &Reference) const
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
   Info FunctionInfo(__func__);
@@ -2756 +3942 @@
   }
 
-  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();
@@ -2766 +3952 @@
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
-    helper = Plane(((*TesselA)->getPosition()),
-                   ((*TesselB)->getPosition()),
-                   ((*TesselC)->getPosition())).getNormal();
+    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++;
@@ -2788 +3974 @@
 
   // construct one orthogonal vector
-  if (!Reference.IsZero()) {
-    AngleZero = (Reference) - (Point->getPosition());
+  if (Reference != NULL) {
+    AngleZero = (*Reference) - (*Point->node);
     AngleZero.ProjectOntoPlane(PlaneNormal);
   }
-  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());
+  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) {
@@ -2811 +3997 @@
   pair<map<double, TesselPoint*>::iterator, bool> InserterTest;
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper = ((*listRunner)->getPosition()) - (Point->getPosition());
+    helper = (*(*listRunner)->node) - (*Point->node);
     helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
@@ -3056 +4242 @@
 
   // copy old location for the volume
-  OldPoint = (point->node->getPosition());
+  OldPoint = (*point->node->node);
 
   // get list of connected points
@@ -3119 +4305 @@
           StartNode--;
           //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
-          Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+          Point = (*(*StartNode)->node) - (*(*MiddleNode)->node);
           StartNode = MiddleNode;
           StartNode++;
@@ -3125 +4311 @@
             StartNode = connectedPath->begin();
           //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
-          Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
-          OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
+          Reference = (*(*StartNode)->node) - (*(*MiddleNode)->node);
+          OrthogonalVector = (*(*MiddleNode)->node) - OldPoint;
           OrthogonalVector.MakeNormalTo(Reference);
           angle = GetAngle(Point, Reference, OrthogonalVector);
@@ -3181 +4367 @@
         AddTesselationTriangle();
         // calculate volume summand as a general tetraeder
-        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
+        volume += CalculateVolumeofGeneralTetraeder(*TPS[0]->node->node, *TPS[1]->node->node, *TPS[2]->node->node, OldPoint);
         // advance number
         count++;
@@ -3566 +4752 @@
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->getPosition(), BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
@@ -3587 +4773 @@
   class BoundaryLineSet *BestLine = NULL;
   for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
-    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());
+    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.)) {