Index: molecuilder/src/tesselation.cpp
===================================================================
--- molecuilder/src/tesselation.cpp	(revision 342fbc9863df30e9e725bc9bc001015a1f173ebd)
+++ molecuilder/src/tesselation.cpp	(revision 086db0fff858be2b9fe495455008d5a091ef3595)
@@ -654,40 +654,4 @@
 };
 
-/** Checks validity of a given sphere of a candidate line.
- * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
- * \param *OtherOptCenter center of the other triangle
- * \param RADIUS radius of sphere
- * \param *LC LinkedCell structure with other atoms
- * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
- */
-bool BoundaryTriangleSet::CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const
-{
-  Info FunctionInfo(__func__);
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
-  bool flag = true;
-
-  // get all points inside the sphere
-  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, OtherOptCenter);
-
-  // remove triangles's endpoints
-  for (int i=0;i<3;i++)
-    ListofPoints->remove(endpoints[i]->node);
-
-  // check for other points
-  if (!ListofPoints->empty()) {
-    Log() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl;
-    flag = false;
-    Log() << Verbose(1) << "External atoms inside of sphere at " << *OtherOptCenter << ":" << endl;
-    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
-      Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
-  }
-  delete(ListofPoints);
-
-  return flag;
-};
-
-
-
 /** Returns the endpoint which is not contained in the given \a *line.
  * \param *line baseline defining two endpoints
@@ -1078,5 +1042,5 @@
   Info FunctionInfo(__func__);
 
-  const double radiusSquared = RADIUS;
+  const double radiusSquared = RADIUS*RADIUS;
   list<const Vector *> VectorList;
   VectorList.push_back(&OptCenter);
@@ -1089,9 +1053,11 @@
   // 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++)
-      if (fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->node) - radiusSquared) < MYEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << "." << endl);
+    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 " << setprecision(13) << distance << "." << endl);
         return false;
       }
+    }
   }
 
@@ -1100,6 +1066,7 @@
     const TesselPoint *Walker = *Runner;
     for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) {
-      if (fabs((*VRunner)->DistanceSquared(Walker->node) - radiusSquared) < MYEPSILON) {
-        DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << Walker << " is out of sphere at " << *(*VRunner) << "." << endl);
+      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 " << setprecision(13) << distance << "." << endl);
         return false;
       }
@@ -1827,4 +1794,5 @@
  * If successful it raises the line count and inserts the new line into the BLS,
  * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
+ * @param *OptCenter desired OptCenter if there are more than one candidate line
  * @param *candidate third point of the triangle to be, for checking between multiple open line candidates
  * @param *a first endpoint
@@ -1832,5 +1800,5 @@
  * @param n index of Tesselation::BLS giving the line with both endpoints
  */
-void Tesselation::AddTesselationLine(const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) {
+void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) {
   bool insertNewLine = true;
 
@@ -1847,9 +1815,14 @@
       // If there is a line with less than two attached triangles, we don't need a new line.
       if (FindLine->second->triangles.size() == 1) {
+        CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
+        if (!Finder->second->pointlist.empty())
+          Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl;
+        else
+          Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate." << endl;
         // get open line
-        CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
-        insertNewLine = false;
-        WinningLine = FindLine->second;
-        if ((!Finder->second->pointlist.empty()) && (*(Finder->second->pointlist.begin()) == candidate->node)) { // stop searching if candidate matches
+        if ((!Finder->second->pointlist.empty()) && (*(Finder->second->pointlist.begin()) == candidate->node)
+            && (OptCenter == NULL || *OptCenter == Finder->second->OptCenter)) { // stop searching if candidate matches
+          insertNewLine = false;
+          WinningLine = FindLine->second;
           break;
         }
@@ -2063,4 +2036,47 @@
 };
 
+
+/** Checks validity of a given sphere of a candidate line.
+ * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
+ * Note that endpoints are stored in Tesselation::TPS.
+ * \param *OtherOptCenter center of the other triangle
+ * \param RADIUS radius of sphere
+ * \param *LC LinkedCell structure with other atoms
+ * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
+ */
+bool Tesselation::CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const
+{
+  Info FunctionInfo(__func__);
+
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
+  bool flag = true;
+
+  cout << Verbose(1) << "Check by: draw sphere {" << OtherOptCenter->x[0] << " " << OtherOptCenter->x[1] << " " << OtherOptCenter->x[2] << "} radius " << RADIUS << " resolution 30" << endl;
+
+  // get all points inside the sphere
+  TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, OtherOptCenter);
+
+  Log() << Verbose(1) << "The following atoms are inside sphere at " << *OtherOptCenter << ":" << endl;
+  for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
+    Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
+
+  // remove triangles's endpoints
+  for (int i=0;i<3;i++)
+    ListofPoints->remove(TPS[i]->node);
+
+  // check for other points
+  if (!ListofPoints->empty()) {
+    Log() << Verbose(1) << "CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere." << endl;
+    flag = false;
+    Log() << Verbose(1) << "External atoms inside of sphere at " << *OtherOptCenter << ":" << endl;
+    for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
+      Log() << Verbose(1) << "  " << *(*Runner) << " with distance " << (*Runner)->node->Distance(OtherOptCenter) << "." << endl;
+  }
+  delete(ListofPoints);
+
+  return flag;
+};
+
+
 /** Checks whether the triangle consisting of the three points is already present.
  * Searches for the points in Tesselation::PointsOnBoundary and checks their
@@ -2182,5 +2198,5 @@
   TesselPoint* Temporary;
   double maxCoordinate[NDIM];
-  BoundaryLineSet BaseLine;
+  BoundaryLineSet *BaseLine = NULL;
   Vector helper;
   Vector Chord;
@@ -2228,23 +2244,24 @@
     NormalVector.Zero();
     NormalVector.x[k] = 1.;
-    BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
-    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
+    BaseLine = new BoundaryLineSet();
+    BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
+    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl;
 
     double ShortestAngle;
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
 
-    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(BaseLine->endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     if (Temporary == NULL)  // have we found a second point?
       continue;
-    BaseLine.endpoints[1] = new BoundaryPointSet(Temporary);
+    BaseLine->endpoints[1] = new BoundaryPointSet(Temporary);
 
     // construct center of circle
-    CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node);
-    CircleCenter.AddVector(BaseLine.endpoints[1]->node->node);
+    CircleCenter.CopyVector(BaseLine->endpoints[0]->node->node);
+    CircleCenter.AddVector(BaseLine->endpoints[1]->node->node);
     CircleCenter.Scale(0.5);
 
     // construct normal vector of circle
-    CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node);
-    CirclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node);
+    CirclePlaneNormal.CopyVector(BaseLine->endpoints[0]->node->node);
+    CirclePlaneNormal.SubtractVector(BaseLine->endpoints[1]->node->node);
 
     double radius = CirclePlaneNormal.NormSquared();
@@ -2264,9 +2281,9 @@
 
     // adding point 1 and point 2 and add the line between them
-    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
-    Log() << Verbose(0) << "Found second point is at " << *BaseLine.endpoints[1]->node << ".\n";
+    Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl;
+    Log() << Verbose(0) << "Found second point is at " << *BaseLine->endpoints[1]->node << ".\n";
 
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    CandidateForTesselation OptCandidates(&BaseLine);
+    CandidateForTesselation OptCandidates(BaseLine);
     FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
     Log() << Verbose(0) << "List of third Points is:" << endl;
@@ -2274,9 +2291,11 @@
         Log() << Verbose(0) << " " << *(*it) << endl;
     }
-
-    BTS = NULL;
-    AddCandidatePolygon(OptCandidates, RADIUS, LC);
-//    delete(BaseLine.endpoints[0]);
-//    delete(BaseLine.endpoints[1]);
+    if (!OptCandidates.pointlist.empty()) {
+      BTS = NULL;
+      AddCandidatePolygon(OptCandidates, RADIUS, LC);
+    } else {
+      delete BaseLine;
+      continue;
+    }
 
     if (BTS != NULL) // we have created one starting triangle
@@ -2286,4 +2305,5 @@
       OptCandidates.pointlist.clear();
     }
+    delete BaseLine;
   }
 };
@@ -2436,5 +2456,5 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
 {
 	Info FunctionInfo(__func__);
@@ -2627,49 +2647,16 @@
   Sprinter++;
   while(Sprinter != connectedClosestPoints->end()) {
-    // add this triangle
-    AddCandidateTriangle(TurningPoint, *Runner, *Sprinter);
-
-    // store values from tesselation procedure therein
-    BTS->SetTopNode(CandidateLine.T);
-    if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
-      BTS->AngleFromTop = CandidateLine.ShortestAngle;
-    else
-      BTS->AngleFromTop = -1.;
-
-    // create normal vector
-    BTS->GetCenter(&Center);
-    Center.SubtractVector(&CandidateLine.OptCenter);
-    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
-    BTS->GetNormalVector(Center);
-
-    // give some verbose output about the whole procedure
-    if (CandidateLine.T != NULL)
-      Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
-    else
-      Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
-
-    // check whether its degenerate and if so, add the other side triangle as well
-    if (BTS->CheckDegeneracy(&CandidateLine.OtherOptCenter, RADIUS, LC)) {
-      Log() << Verbose(1) << "Triangle " << *BTS << " is degenerated, adding other side." << endl;
-      BoundaryTriangleSet * const triangle = BTS;
-      AddDegeneratedTriangle(BTS, &CandidateLine.OtherOptCenter, RADIUS, LC);
-      // check whether other side triangle has been created
-      if (BTS != NULL) {
-        BTS->SetTopNode(CandidateLine.T);
-        if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
-          BTS->AngleFromTop = CandidateLine.ShortestAngle;
-        else
-          BTS->AngleFromTop = -1.;
-        BTS->SphereCenter.CopyVector(&CandidateLine.OtherOptCenter);
-        // normal vector points other way round
-        BTS->NormalVector.CopyVector(&triangle->NormalVector);
-        BTS->NormalVector.Scale(-1.);
-
-        // give some verbose output about the whole procedure
-        if (CandidateLine.T != NULL)
-          Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
-        else
-          Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
-      }
+    AddTesselationPoint(TurningPoint, 0);
+    AddTesselationPoint(*Runner, 1);
+    AddTesselationPoint(*Sprinter, 2);
+
+    // check whether we add a degenerate or a normal triangle
+    if (CheckDegeneracy(&CandidateLine.OtherOptCenter, RADIUS, LC)) {
+      // add normal and degenerate triangles
+      Log() << Verbose(1) << "Triangle of endpoints " << *TPS[0] << "," << *TPS[1] << " and " << *TPS[2] << " is degenerated, adding both sides." << endl;
+      BoundaryTriangleSet * const triangle = AddDegeneratedTriangle(CandidateLine, RADIUS, LC);
+    } else {
+      // add this triangle
+      AddCandidateTriangle(CandidateLine);
     }
 
@@ -2683,95 +2670,218 @@
 };
 
-/** If a given \a *triangle is degenerated, this adds the "other side".
+/** If a given \a *triangle is degenerated, this adds both sides.
  * i.e. the triangle with same BoundaryPointSet's but NormalVector in opposite direction.
- * \param *triangle degenerated triangle
- * \param *OtherOptCenter center of the sphere for the other triangle
+ * Note that endpoints are stored in Tesselation::TPS
+ * \param CandidateLine CanddiateForTesselation structure for the desired BoundaryLine
  * \param RADIUS radius of sphere
  * \param *LC pointer to LinkedCell structure
  * \return pointer to created other triangle
  */
-void Tesselation::AddDegeneratedTriangle(BoundaryTriangleSet *&triangle, Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC)
+BoundaryTriangleSet * const  Tesselation::AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC)
 {
   Info FunctionInfo(__func__);
-  TesselPoint *ThirdNode = NULL;
-  BoundaryTriangleSet *otherTriangle = NULL;
-  /// 1. For each of the open sides of \a *triangle find best candidate
+  BoundaryTriangleSet *triangle = NULL;
+  pair<LineMap::iterator,LineMap::iterator> FindPair[3];
+  int LineCase[3];
+  int count[3];
+  Vector Center;
+  CandidateMap::const_iterator CandidateCheck = OpenLines.end();
+
+  /// 1. For each new pair of endpoints, find and count the number of present lines
+  Log() << Verbose(0) << "INFO: Counting present open lines ..." << endl;
   for (int i=0;i<3;i++) {
     BLS[i] = NULL;
-    // go through all possible lines (also degenerate ones) and find ones to be used
-    pair<LineMap::iterator,LineMap::iterator> FindPair = triangle->lines[i]->endpoints[0]->lines.equal_range(triangle->lines[i]->endpoints[1]->Nr);
-    for (LineMap::iterator FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
-      if ((*FindLine).second->triangles.size() == 1) {
-        Log() << Verbose(1) << "Found line " << *(*FindLine).second << " from triangle " << *((*FindLine).second->triangles.begin()->second) << "." << endl;
-        // check whether best Candidate for this line from the degenerate triangle is the third endpoint
-        CandidateForTesselation CandidateLine((*FindLine).second);
-        otherTriangle = ((*FindLine).second->triangles.begin())->second;  // this is unique and present, above we check that there is exactly one
-        ThirdNode = otherTriangle->GetThirdEndpoint((*FindLine).second)->node;
-        FindNextSuitableTriangle(CandidateLine, *triangle, RADIUS, LC);
-
-        /// 2. If best candidate is third endpoint,
-        for (TesselPointList::const_iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); ++Runner) {
-          Log() << Verbose(1) << "Found candidate " << *(*Runner) << endl;
-          if ((*Runner) == ThirdNode) { // seems like we found a line to connect to
-            /// sets Tesselation::BLS[i] accordingly and removes the line from Tesselation::OpenLines
-            AddExistingTesselationTriangleLine((*FindLine).second,i);
-            break;
-          }
-        }
+    LineCase[i] = -1;
+    // find already present lines
+    FindPair[i] = TPS[(i+0)%3]->lines.equal_range(TPS[(i+1)%3]->Nr);
+    // count their number
+    count[i] = 0;
+    Log() << Verbose(1) << "For line between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    for (LineMap::iterator FindLine = FindPair[i].first; FindLine != FindPair[i].second; FindLine++)
+      if (FindLine->second->triangles.size() == 1) {
+        Log() << Verbose(1) << "  open line " << *FindLine->second << "." << endl;
+        count[i]++;
       }
-    }
-    /// otherwise we create a new BoundaryLineSet for this Tesselation::BLS[i] and add it to OpenLines
-  }
-  if ((BLS[0] == NULL) && (BLS[1] == NULL) && (BLS[2] == NULL)) {
-    Log() << Verbose(0) << "This degenerate triangle would not share any line with its other side, rejecting." << endl;
-    BTS = NULL;
-  } else {
-    // created new lines for the other egdes of the other side triangle
-    for (int i=0;i<3;i++) {
-      if (BLS[i] == NULL) {
-        AddNewTesselationTriangleLine(triangle->lines[i]->endpoints[0],triangle->lines[i]->endpoints[1],i);
-      }
-    }
-    /// 3. Create the triangle
-    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-    AddTesselationTriangle();
-
-    // add triangle to the CandidateForTesselation of its open lines
-    CandidateMap::iterator Finder = OpenLines.end();
-    for (int i=0;i<3;i++) {
-      Finder = OpenLines.find(BTS->lines[i]);
-      if (Finder != OpenLines.end()) {
-        Log() << Verbose(1) << "INFO: Found open line of this triangle " << *(Finder->first) << endl;
-        if (Finder->second->T == NULL) {
-          Finder->second->T = BTS;
-          Log() << Verbose(1) << "ACCEPT: Setting source triangle of open line " << *(Finder->first) << " to the newly created one." << endl;
-        }
-      } else
-        Log() << Verbose(1) << "REJECT: line " << *BTS->lines[i] << " of this triangle is not found to be open." << endl;
-    }
-  }
+    Log() << Verbose(1) << "  SUM: there are " << count[i] << " line(s) present."  << endl;
+  }
+
+  /// 2. Create or pick the lines for the first triangle
+  // for each count we have a different case:
+  // case 0: no triangles at this line and not closed
+  // case 1: no triangles at new line is closed
+  // case 2: one line with one triangle attached, one will used this line
+  // case 3: one line with two triangles attached, new line will be used by both degenerate triangles
+  // case 4: two lines with one triangle each, each new triangle will pick one
+  Log() << Verbose(0) << "INFO: Creating/Picking lines for first triangle ..." << endl;
+  for (int i=0;i<3;i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    switch (count[i]) {
+      case 0:
+        AddNewTesselationTriangleLine(TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      case 1:
+//        // check whether present line has right candidate
+//        CandidateCheck = OpenLines.find(FindPair[i].first->second);
+//        if (CandidateCheck != OpenLines.end()) {
+//          if (!CandidateCheck->second->pointlist.empty()) {
+//            for (TesselPointList::const_iterator Runner = CandidateCheck->second->pointlist.begin(); Runner != CandidateCheck->second->pointlist.end(); ++Runner) {
+//              Log() << Verbose(1) << "Found candidate " << *(*Runner) << " with Candidate center is " << CandidateCheck->second->OptCenter << ", desired opt center is " << CandidateLine.OptCenter << "." << endl;
+//              if (((*Runner) == TPS[(i+2)%3]->node) && (CandidateLine.OptCenter == CandidateCheck->second->OptCenter)) { // seems like we found a line to connect to
+//                AddExistingTesselationTriangleLine(CandidateCheck->first, i);
+//                break;  // have to break as a AddExistingTesselationTriangleLine changes OpenLines
+//              }
+//            }
+//          } else {
+//            Log() << Verbose(1) << "Open line between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << " has no candidate." << endl;
+//          }
+//        } else {
+//          Log() << Verbose(1) << "No open lines present for line between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << "." << endl;
+//        }
+//        if (BLS[i] == NULL) {
+//          AddNewTesselationTriangleLine(TPS[(i+0)%3], TPS[(i+1)%3], i);
+//        }
+        AddTesselationLine(&CandidateLine.OptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      default:
+      case 2:
+        AddTesselationLine(&CandidateLine.OptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+//        DoeLog(0) && (eLog() << Verbose(0) << "There are more than two BoundaryPointSets between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << "." << endl);
+//        performCriticalExit();
+//        break;
+    }
+  }
+
+  /// 3. create the first triangle and NormalVector and so on
+  Log() << Verbose(0) << "INFO: Adding first triangle with center at " << CandidateLine.OptCenter << " ..." << endl;
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTesselationTriangle();
+  // check whether triangle has been created
+  BTS->SetTopNode(CandidateLine.T);
+  if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
+    BTS->AngleFromTop = CandidateLine.ShortestAngle;
+  else
+    BTS->AngleFromTop = -1.;
+  // create normal vector
+  BTS->GetCenter(&Center);
+  Center.SubtractVector(&CandidateLine.OptCenter);
+  BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+  BTS->GetNormalVector(Center);
+  // give some verbose output about the whole procedure
+  if (CandidateLine.T != NULL)
+    Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
+  else
+    Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
+  triangle = BTS;
+
+  /// 4. Gather candidates for each new line
+  Log() << Verbose(0) << "INFO: Adding candidates to new lines ..." << endl;
+  for (int i=0;i<3;i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    CandidateCheck = OpenLines.find(BLS[i]);
+    if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
+      if (CandidateCheck->second->T == NULL)
+        CandidateCheck->second->T = triangle;
+      FindNextSuitableTriangle(*(CandidateCheck->second), *CandidateCheck->second->T, RADIUS, LC);
+    }
+  }
+
+  /// 5. Create or pick the lines for the second triangle
+  Log() << Verbose(0) << "INFO: Creating/Picking lines for second triangle ..." << endl;
+  for (int i=0;i<3;i++) {
+    Log() << Verbose(0) << "Current line is between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << ":" << endl;
+    switch (count[i]) {
+      case 0:
+        AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      case 1:
+//        // check whether present line has right candidate
+//        CandidateCheck = OpenLines.find(BLS[i]);
+//        BLS[i] = NULL;
+//        if (CandidateCheck == OpenLines.end()) {  // first triangle has closed present line
+//          AddNewTesselationTriangleLine(TPS[(i+0)%3], TPS[(i+1)%3], i);
+//        } else {
+//          for (TesselPointList::const_iterator Runner = CandidateCheck->second->pointlist.begin(); Runner != CandidateCheck->second->pointlist.end(); ++Runner) {
+//            Log() << Verbose(1) << "Found candidate " << *(*Runner) << " with Candidate center is " << CandidateCheck->second->OptCenter << ", desired other opt center is " << CandidateLine.OtherOptCenter << "." << endl;
+//            if (((*Runner) == TPS[(i+2)%3]->node) && (CandidateLine.OtherOptCenter == CandidateCheck->second->OptCenter)) { // seems like we found a line to connect to
+//              AddExistingTesselationTriangleLine(CandidateCheck->first, i);
+//              break;  // have to break as a AddExistingTesselationTriangleLine changes OpenLines
+//            }
+//          }
+//        }
+//        if (BLS[i] == NULL) {
+//          DoeLog(0) && (eLog() << Verbose(0) << "BoundaryLine between " << TPS[(i+0)%3] << " and " << TPS[(i+1)%3] << " has not the correct candidate, i.e. the third endpoint: " << *TPS[(i+2)%3] << endl);
+//          performCriticalExit();
+//        }
+        AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+      case 2:
+      default:
+        AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i+2)%3], TPS[(i+0)%3], TPS[(i+1)%3], i);
+        break;
+//        DoeLog(0) && (eLog() << Verbose(0) << "There are more than two BoundaryPointSets between " << *TPS[(i+0)%3] << " and " << *TPS[(i+1)%3] << "." << endl);
+//        performCriticalExit();
+//        break;
+    }
+  }
+
+  /// 6. create the second triangle and NormalVector and so on
+  Log() << Verbose(0) << "INFO: Adding second triangle with center at " << CandidateLine.OtherOptCenter << " ..." << endl;
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTesselationTriangle();
+  // check whether triangle has been created
+  BTS->SetTopNode(CandidateLine.T);
+  if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
+    BTS->AngleFromTop = CandidateLine.ShortestAngle;
+  else
+    BTS->AngleFromTop = -1.;
+  BTS->SphereCenter.CopyVector(&CandidateLine.OtherOptCenter);
+  // create normal vector in other direction
+  BTS->GetNormalVector(&triangle->NormalVector);
+  BTS->NormalVector.Scale(-1.);
+  // give some verbose output about the whole procedure
+  if (CandidateLine.T != NULL)
+    Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
+  else
+    Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
+
+  return triangle;
 };
 
 /** Adds a triangle to the Tesselation structure from three given TesselPoint's.
- * \param *first first TesselPoint
- * \param *second second TesselPoint
- * \param *third third TesselPoint
- */
-void Tesselation::AddCandidateTriangle(TesselPoint *first, TesselPoint *second, TesselPoint *third)
+ * Note that endpoints are in Tesselation::TPS.
+ * \param CandidateLine CandidateForTesselation structure contains other information
+ */
+void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine)
 {
   Info FunctionInfo(__func__);
-  // add the points
-  AddTesselationPoint(first, 0);
-  AddTesselationPoint(second, 1);
-  AddTesselationPoint(third, 2);
+  Vector Center;
 
   // add the lines
-  AddTesselationLine(TPS[2], TPS[0], TPS[1], 0);
-  AddTesselationLine(TPS[1], TPS[0], TPS[2], 1);
-  AddTesselationLine(TPS[0], TPS[1], TPS[2], 2);
+  AddTesselationLine(&CandidateLine.OptCenter, TPS[2], TPS[0], TPS[1], 0);
+  AddTesselationLine(&CandidateLine.OptCenter, TPS[1], TPS[0], TPS[2], 1);
+  AddTesselationLine(&CandidateLine.OptCenter, TPS[0], TPS[1], TPS[2], 2);
 
   // add the triangles
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   AddTesselationTriangle();
+  // store values from tesselation procedure therein
+  BTS->SetTopNode(CandidateLine.T);
+  if (CandidateLine.T != NULL)  // start triangle has angle from top of -1
+    BTS->AngleFromTop = CandidateLine.ShortestAngle;
+  else
+    BTS->AngleFromTop = -1.;
+
+  // create normal vector
+  BTS->GetCenter(&Center);
+  Center.SubtractVector(&CandidateLine.OptCenter);
+  BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+  BTS->GetNormalVector(Center);
+
+  // give some verbose output about the whole procedure
+  if (CandidateLine.T != NULL)
+    Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl;
+  else
+    Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl;
 };
 
@@ -3163,5 +3273,5 @@
  * @param *LC LinkedCell structure with neighbouring points
  */
-void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
+void Tesselation::FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
 {
 	Info FunctionInfo(__func__);
@@ -3344,5 +3454,5 @@
   }
 
-  if (!CandidateLine.CheckValidity(RADIUS, LC)) {
+  if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
     DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl);
     performCriticalExit();
@@ -4349,8 +4459,8 @@
         AddTesselationPoint(*EndNode, 2);
         Log() << Verbose(3) << "Adding new triangle lines."<< endl;
-        AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-        AddTesselationLine(NULL, TPS[0], TPS[2], 1);
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+        AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
         NewLines.push_back(BLS[1]);
-        AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+        AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
         BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
         BTS->GetNormalVector(NormalVector);
@@ -4801,7 +4911,7 @@
   AddTesselationPoint(point, 2);
   Log() << Verbose(2) << "Adding new triangle lines."<< endl;
-  AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-  AddTesselationLine(NULL, TPS[0], TPS[2], 1);
-  AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   BTS->GetNormalVector(TempTriangle->NormalVector);
@@ -4816,7 +4926,7 @@
   AddTesselationPoint(point, 2);
   Log() << Verbose(2) << "Adding new triangle lines."<< endl;
-  AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-  AddTesselationLine(NULL, TPS[0], TPS[2], 1);
-  AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+  AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+  AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   BTS->GetNormalVector(TempTriangle->NormalVector);
@@ -5051,7 +5161,7 @@
       for (int i = 0; i < 3; i++)
         AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i);
-      AddTesselationLine(NULL, TPS[0], TPS[1], 0);
-      AddTesselationLine(NULL, TPS[0], TPS[2], 1);
-      AddTesselationLine(NULL, TPS[1], TPS[2], 2);
+      AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+      AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
+      AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
       if (TriangleNrs.empty())
         DoeLog(0) && (eLog()<< Verbose(0) << "No more free triangle numbers!" << endl);
Index: molecuilder/src/tesselation.hpp
===================================================================
--- molecuilder/src/tesselation.hpp	(revision 342fbc9863df30e9e725bc9bc001015a1f173ebd)
+++ molecuilder/src/tesselation.hpp	(revision 086db0fff858be2b9fe495455008d5a091ef3595)
@@ -163,5 +163,4 @@
     bool IsPresentTupel(const BoundaryPointSet * const Points[3]) const;
     bool IsPresentTupel(const BoundaryTriangleSet * const T) const;
-    bool CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const;
     void SetTopNode(const BoundaryTriangleSet * const topnode);
 
@@ -283,15 +282,16 @@
     void AddTesselationPoint(TesselPoint* Candidate, const int n);
     void SetTesselationPoint(TesselPoint* Candidate, const int n) const;
-    void AddTesselationLine(const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
+    void AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
     void AddExistingTesselationTriangleLine(class BoundaryLineSet *FindLine, int n);
     void AddTesselationTriangle();
     void AddTesselationTriangle(const int nr);
-    void AddCandidateTriangle(TesselPoint *first, TesselPoint *second, TesselPoint *third);
-    void AddDegeneratedTriangle(BoundaryTriangleSet *&triangle, Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC);
+    void AddCandidateTriangle(CandidateForTesselation &CandidateLine);
+    BoundaryTriangleSet * const AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC);
     void AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC);
     void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle);
     void RemoveTesselationLine(class BoundaryLineSet *line);
     void RemoveTesselationPoint(class BoundaryPointSet *point);
+    bool CheckDegeneracy(Vector *OtherOptCenter, const double RADIUS, const LinkedCell *LC) const;
 
 
@@ -299,6 +299,6 @@
     void FindStartingTriangle(const double RADIUS, const LinkedCell *LC);
     void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC);
-    void FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const;
-    bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
+    void FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const;
+    bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
     int CheckPresenceOfTriangle(class TesselPoint *Candidates[3]) const;
     class BoundaryTriangleSet * GetPresentTriangle(TesselPoint *Candidates[3]);
