Context Navigation

← Previous Change
Next Change →

Changeset ab1932 for src/boundary.cpp

Timestamp:

Aug 3, 2009, 8:21:05 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

03e57a

Parents:

0dbddc (diff), edb93c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'TesselationRefactoring' into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/molecules.hpp

All of Saskia Metzler's new function were transfered from boundary.cpp to tesselation.cpp and the changes due to TesselPoint, LinkedCell and so on incorporated.

File:

: 1 edited

src/boundary.cpp (modified) (15 diffs)

Legend:

: Unmodified
: Added
: Removed

src/boundary.cpp

-              r0dbddc
+              rab1932
+/** \file boundary.hpp
+ *
+ * Implementations and super-function for envelopes
+ */
 #include "boundary.hpp"
+#include "linkedcell.hpp"
+#include "molecules.hpp"
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_permutation.h>
+#define DEBUG 1
+#define DoSingleStepOutput 0
+#define DoTecplotOutput 1
+#define DoRaster3DOutput 1
+#define DoVRMLOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
+#define HULLEPSILON 1e-7
+// ======================================== Points on Boundary =================================
+BoundaryPointSet::BoundaryPointSet()
+{
+  LinesCount = 0;
+  Nr = -1;
+}
+;
+BoundaryPointSet::BoundaryPointSet(atom *Walker)
+{
+  node = Walker;
+  LinesCount = 0;
+  Nr = Walker->nr;
+}
+;
+BoundaryPointSet::~BoundaryPointSet()
+{
+  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  if (!lines.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
+  node = NULL;
+}
+;
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+{
+  cout << Verbose(6) << "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++;
+}
+;
+ostream &
+operator <<(ostream &ost, BoundaryPointSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  return ost;
+}
+;
+// ======================================== Lines on Boundary =================================
+BoundaryLineSet::BoundaryLineSet()
+{
+  for (int i = 0; i < 2; i++)
+    endpoints[i] = NULL;
+  TrianglesCount = 0;
+  Nr = -1;
+}
+;
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
+{
+  // 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); //
+  // clear triangles list
+  TrianglesCount = 0;
+  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+}
+;
+BoundaryLineSet::~BoundaryLineSet()
+{
+  int Numbers[2];
+  Numbers[0] = endpoints[1]->Nr;
+  Numbers[1] = endpoints[0]->Nr;
+  for (int i = 0; i < 2; i++) {
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+    // 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) {
+        endpoints[i]->lines.erase(Runner);
+        break;
+      }
+    if (endpoints[i]->lines.empty()) {
+      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+      if (endpoints[i] != NULL) {
+        delete(endpoints[i]);
+        endpoints[i] = NULL;
+      } else
+        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
+  if (!triangles.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+}
+;
+void
+BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+{
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
+      << endl;
+  triangles.insert(TrianglePair(triangle->Nr, triangle));
+  TrianglesCount++;
+}
+;
+/** 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(class BoundaryLineSet *line)
+{
+  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 normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
+ * 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(ofstream *out)
+{
+  Vector BaseLineNormal;
+  double angle = 0;
+  // get the two triangles
+  if (TrianglesCount != 2) {
+    *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
+    return false;
+  }
+  // have a normal vector on the base line pointing outwards
+  BaseLineNormal.Zero();
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    BaseLineNormal.AddVector(&runner->second->NormalVector);
+  BaseLineNormal.Normalize();
+  // and calculate the sum of the angles with this normal vector and each of the triangle ones'
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    angle += BaseLineNormal.Angle(&runner->second->NormalVector);
+  if ((angle - M_PI) > -MYEPSILON)
+    return true;
+  else
+    return false;
+}
+/** Checks whether point is any of the two endpoints this line contains.
+ * \param *point point to test
+ * \return true - point is of the line, false - is not
+ */
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+  for(int i=0;i<2;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
+ostream &
+operator <<(ostream &ost, BoundaryLineSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "]";
+  return ost;
+}
+;
+// ======================================== Triangles on Boundary =================================
+BoundaryTriangleSet::BoundaryTriangleSet()
+{
+  for (int i = 0; i < 3; i++)
+    {
+      endpoints[i] = NULL;
+      lines[i] = NULL;
+    }
+  Nr = -1;
+}
+;
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
+{
+  // set number
+  Nr = number;
+  // set lines
+  cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
+  for (int i = 0; i < 3; i++)
+    {
+      lines[i] = line[i];
+      lines[i]->AddTriangle(this);
+    }
+  // get ascending order of endpoints
+  map<int, class BoundaryPointSet *> 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;
+  cout << Verbose(6) << " with end points ";
+  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
+      != OrderMap.end(); runner++)
+    {
+      endpoints[Counter] = runner->second;
+      cout << " " << *endpoints[Counter];
+      Counter++;
+    }
+  if (Counter < 3)
+    {
+      cerr << "ERROR! We have a triangle with only two distinct endpoints!"
+          << endl;
+      //exit(1);
+    }
+  cout << "." << endl;
+}
+;
+BoundaryTriangleSet::~BoundaryTriangleSet()
+{
+  for (int i = 0; i < 3; i++) {
+    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+    lines[i]->triangles.erase(Nr);
+    if (lines[i]->triangles.empty()) {
+      if (lines[i] != NULL) {
+        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+        delete (lines[i]);
+        lines[i] = NULL;
+      } else
+        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << 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(Vector &OtherVector)
+{
+  // get normal vector
+  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
+      &endpoints[2]->node->x);
+  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
+  if (NormalVector.Projection(&OtherVector) > 0)
+    NormalVector.Scale(-1.);
+};
+/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+ * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+ * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
+ * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
+ * smaller than the first line.
+ * \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(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
+{
+  Vector CrossPoint;
+  Vector helper;
+  int i=0;
+  if (Intersection->GetIntersectionWithPlane(out, &NormalVector, &endpoints[0]->node->x, MolCenter, x)) {
+    *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
+    return false;
+  }
+  // Calculate cross point between one baseline and the line from the third endpoint to intersection
+  do {
+    CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, &endpoints[i%3]->node->x, &endpoints[(i+1)%3]->node->x, &endpoints[(i+2)%3]->node->x, Intersection);
+    helper.CopyVector(&endpoints[(i+1)%3]->node->x);
+    helper.SubtractVector(&endpoints[i%3]->node->x);
+    i++;
+    if (i>3)
+      break;
+  } while (CrossPoint.NormSquared() < MYEPSILON);
+  if (i>3) {
+    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
+    exit(255);
+  }
+  CrossPoint.SubtractVector(&endpoints[i%3]->node->x);
+  // check whether intersection is inside or not by comparing length of intersection and length of cross point
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
+    return true;
+  } else { // outside!
+    Intersection->Zero();
+    return false;
+  }
+};
+/** 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(class BoundaryLineSet *line)
+{
+  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(class BoundaryPointSet *point)
+{
+  for(int i=0;i<3;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
+ostream &
+operator <<(ostream &ost, BoundaryTriangleSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  return ost;
+}
+;
+// ============================ CandidateForTesselation =============================
+CandidateForTesselation::CandidateForTesselation(
+  atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
+) {
+  point = candidate;
+  BaseLine = line;
+  OptCenter.CopyVector(&OptCandidateCenter);
+  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
+}
+CandidateForTesselation::~CandidateForTesselation() {
+  point = NULL;
+  BaseLine = NULL;
+}
+#include<gsl/gsl_poly.h>
 // ========================================== F U N C T I O N S =================================
-/** Finds the endpoint two lines are sharing.
- * \param *line1 first line
- * \param *line2 second line
- * \return point which is shared or NULL if none
- */
-class BoundaryPointSet *
-GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
+{
-  class BoundaryLineSet * lines[2] =
-    { line1, line2 };
-  class BoundaryPointSet *node = NULL;
-  map<int, class BoundaryPointSet *> OrderMap;
-  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
-  for (int i = 0; i < 2; i++)
-    // for both lines
-    for (int j = 0; j < 2; j++)
-      { // for both endpoints
-        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
-            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
-        if (!OrderTest.second)
-          { // if insertion fails, we have common endpoint
-            node = OrderTest.first->second;
-            cout << Verbose(5) << "Common endpoint of lines " << *line1
-                << " and " << *line2 << " is: " << *node << "." << endl;
-            j = 2;
-            i = 2;
-            break;
+          }
+      }
-  return node;
+}
+;
-/** Determines the boundary points of a cluster.
- * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
- * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
- * center and first and last point in the triple, it is thrown out.
- * \param *out output stream for debugging
- * \param *mol molecule structure representing the cluster
- */
-Boundaries *
-GetBoundaryPoints(ofstream *out, molecule *mol)
+{
-  atom *Walker = NULL;
-  PointMap PointsOnBoundary;
-  LineMap LinesOnBoundary;
-  TriangleMap TrianglesOnBoundary;
-  Vector *MolCenter = mol->DetermineCenterOfAll(out);
-  Vector helper;
-  *out << Verbose(1) << "Finding all boundary points." << endl;
-  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
-  BoundariesTestPair BoundaryTestPair;
-  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
-  double radius, angle;
-  // 3a. Go through every axis
-  for (int axis = 0; axis < NDIM; axis++) {
-    AxisVector.Zero();
-    AngleReferenceVector.Zero();
-    AngleReferenceNormalVector.Zero();
-    AxisVector.x[axis] = 1.;
-    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
-    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
-    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
-    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
-    Walker = mol->start;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      Vector ProjectedVector;
-      ProjectedVector.CopyVector(&Walker->x);
-      ProjectedVector.SubtractVector(MolCenter);
-      ProjectedVector.ProjectOntoPlane(&AxisVector);
-      // correct for negative side
-      radius = ProjectedVector.NormSquared();
-      if (fabs(radius) > MYEPSILON)
-        angle = ProjectedVector.Angle(&AngleReferenceVector);
-      else
-        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
-      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
-      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
-        angle = 2. * M_PI - angle;
+      }
-      *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
-      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
-      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
-        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
-        *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
-        *out << Verbose(2) << "New vector: " << *Walker << endl;
-        double tmp = ProjectedVector.NormSquared();
-        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
-          BoundaryTestPair.first->second.first = tmp;
-          BoundaryTestPair.first->second.second = Walker;
-          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
-        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
-          helper.CopyVector(&Walker->x);
-          helper.SubtractVector(MolCenter);
-          tmp = helper.NormSquared();
-          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
-          helper.SubtractVector(MolCenter);
-          if (helper.NormSquared() < tmp) {
-            BoundaryTestPair.first->second.second = Walker;
-            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
-          } else {
-            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
+          }
-        } else {
-          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
+        }
+      }
+    }
-    // printing all inserted for debugging
-    //    {
-    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
-    //      int i=0;
-    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-    //        if (runner != BoundaryPoints[axis].begin())
-    //          *out << ", " << i << ": " << *runner->second.second;
-    //        else
-    //          *out << i << ": " << *runner->second.second;
-    //        i++;
-    //      }
-    //      *out << endl;
-    //    }
-    // 3c. throw out points whose distance is less than the mean of left and right neighbours
-    bool flag = false;
-    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
-    do { // do as long as we still throw one out per round
-      flag = false;
-      Boundaries::iterator left = BoundaryPoints[axis].end();
-      Boundaries::iterator right = BoundaryPoints[axis].end();
-      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-        // set neighbours correctly
-        if (runner == BoundaryPoints[axis].begin()) {
-          left = BoundaryPoints[axis].end();
-        } else {
-          left = runner;
+        }
-        left--;
-        right = runner;
-        right++;
-        if (right == BoundaryPoints[axis].end()) {
-          right = BoundaryPoints[axis].begin();
+        }
-        // check distance
-        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+        {
-          Vector SideA, SideB, SideC, SideH;
-          SideA.CopyVector(&left->second.second->x);
-          SideA.SubtractVector(MolCenter);
-          SideA.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideA: ";
-          //          SideA.Output(out);
-          //          *out << endl;
-          SideB.CopyVector(&right->second.second->x);
-          SideB.SubtractVector(MolCenter);
-          SideB.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideB: ";
-          //          SideB.Output(out);
-          //          *out << endl;
-          SideC.CopyVector(&left->second.second->x);
-          SideC.SubtractVector(&right->second.second->x);
-          SideC.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideC: ";
-          //          SideC.Output(out);
-          //          *out << endl;
-          SideH.CopyVector(&runner->second.second->x);
-          SideH.SubtractVector(MolCenter);
-          SideH.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideH: ";
-          //          SideH.Output(out);
-          //          *out << endl;
-          // calculate each length
-          double a = SideA.Norm();
-          //double b = SideB.Norm();
-          //double c = SideC.Norm();
-          double h = SideH.Norm();
-          // calculate the angles
-          double alpha = SideA.Angle(&SideH);
-          double beta = SideA.Angle(&SideC);
-          double gamma = SideB.Angle(&SideH);
-          double delta = SideC.Angle(&SideH);
-          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
-          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
-          *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
-          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
-            // throw out point
-            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
-            BoundaryPoints[axis].erase(runner);
-            flag = true;
+          }
+        }
+      }
-    } while (flag);
+  }
-  delete(MolCenter);
-  return BoundaryPoints;
+}
+;
 /** Determines greatest diameters of a cluster defined by its convex envelope.
 …
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = BoundaryPtr;
+  if (BoundaryPoints == NULL)
+    {
+      BoundaryFreeFlag = true;
+      BoundaryPoints = GetBoundaryPoints(out, mol);
+    }
+  else
+    {
+      *out << Verbose(1) << "Using given boundary points set." << endl;
+    }
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
 …
       for (i=0;i<3;i++) { // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
           *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
         *vrmlfile << "\t";
+      }
 …
       for (i=0;i<3;i++) {  // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
           *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
         *rasterfile << "\t";
+      }
 …
     *out << Verbose(2) << "The following triangles were created:";
     int Counter = 1;
     atom *Walker = NULL;
+    TesselPoint *Walker = NULL;
     for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
       Walker = target->second->node;
       LookupList[Walker->nr] = Counter++;
       *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
+      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
+    }
     *tecplot << endl;
 …
+  }
+}
+/** Determines the boundary points of a cluster.
+ * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
+ * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
+ * center and first and last point in the triple, it is thrown out.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure representing the cluster
+ */
+Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol)
+{
+  atom *Walker = NULL;
+  PointMap PointsOnBoundary;
+  LineMap LinesOnBoundary;
+  TriangleMap TrianglesOnBoundary;
+  Vector *MolCenter = mol->DetermineCenterOfAll(out);
+  Vector helper;
+  *out << Verbose(1) << "Finding all boundary points." << endl;
+  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
+  BoundariesTestPair BoundaryTestPair;
+  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
+  double radius, angle;
+  // 3a. Go through every axis
+  for (int axis = 0; axis < NDIM; axis++) {
+    AxisVector.Zero();
+    AngleReferenceVector.Zero();
+    AngleReferenceNormalVector.Zero();
+    AxisVector.x[axis] = 1.;
+    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
+    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
+    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+    Walker = mol->start;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      Vector ProjectedVector;
+      ProjectedVector.CopyVector(&Walker->x);
+      ProjectedVector.SubtractVector(MolCenter);
+      ProjectedVector.ProjectOntoPlane(&AxisVector);
+      // correct for negative side
+      radius = ProjectedVector.NormSquared();
+      if (fabs(radius) > MYEPSILON)
+        angle = ProjectedVector.Angle(&AngleReferenceVector);
+      else
+        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
+        angle = 2. * M_PI - angle;
+      }
+      *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
+      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
+      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
+        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
+        *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
+        *out << Verbose(2) << "New vector: " << *Walker << endl;
+        double tmp = ProjectedVector.NormSquared();
+        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
+          BoundaryTestPair.first->second.first = tmp;
+          BoundaryTestPair.first->second.second = Walker;
+          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
+        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
+          helper.CopyVector(&Walker->x);
+          helper.SubtractVector(MolCenter);
+          tmp = helper.NormSquared();
+          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
+          helper.SubtractVector(MolCenter);
+          if (helper.NormSquared() < tmp) {
+            BoundaryTestPair.first->second.second = Walker;
+            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
+          } else {
+            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
+          }
+        } else {
+          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
+        }
+      }
+    }
+    // printing all inserted for debugging
+    //    {
+    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+    //      int i=0;
+    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+    //        if (runner != BoundaryPoints[axis].begin())
+    //          *out << ", " << i << ": " << *runner->second.second;
+    //        else
+    //          *out << i << ": " << *runner->second.second;
+    //        i++;
+    //      }
+    //      *out << endl;
+    //    }
+    // 3c. throw out points whose distance is less than the mean of left and right neighbours
+    bool flag = false;
+    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
+    do { // do as long as we still throw one out per round
+      flag = false;
+      Boundaries::iterator left = BoundaryPoints[axis].end();
+      Boundaries::iterator right = BoundaryPoints[axis].end();
+      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        // set neighbours correctly
+        if (runner == BoundaryPoints[axis].begin()) {
+          left = BoundaryPoints[axis].end();
+        } else {
+          left = runner;
+        }
+        left--;
+        right = runner;
+        right++;
+        if (right == BoundaryPoints[axis].end()) {
+          right = BoundaryPoints[axis].begin();
+        }
+        // check distance
+        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+        {
+          Vector SideA, SideB, SideC, SideH;
+          SideA.CopyVector(&left->second.second->x);
+          SideA.SubtractVector(MolCenter);
+          SideA.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideA: ";
+          //          SideA.Output(out);
+          //          *out << endl;
+          SideB.CopyVector(&right->second.second->x);
+          SideB.SubtractVector(MolCenter);
+          SideB.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideB: ";
+          //          SideB.Output(out);
+          //          *out << endl;
+          SideC.CopyVector(&left->second.second->x);
+          SideC.SubtractVector(&right->second.second->x);
+          SideC.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideC: ";
+          //          SideC.Output(out);
+          //          *out << endl;
+          SideH.CopyVector(&runner->second.second->x);
+          SideH.SubtractVector(MolCenter);
+          SideH.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideH: ";
+          //          SideH.Output(out);
+          //          *out << endl;
+          // calculate each length
+          double a = SideA.Norm();
+          //double b = SideB.Norm();
+          //double c = SideC.Norm();
+          double h = SideH.Norm();
+          // calculate the angles
+          double alpha = SideA.Angle(&SideH);
+          double beta = SideA.Angle(&SideC);
+          double gamma = SideB.Angle(&SideH);
+          double delta = SideC.Angle(&SideH);
+          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
+          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+          *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
+            // throw out point
+            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+            BoundaryPoints[axis].erase(runner);
+            flag = true;
+          }
+        }
+      }
+    } while (flag);
+  }
+  delete(MolCenter);
+  return BoundaryPoints;
+};
 /** Tesselates the convex boundary by finding all boundary points.
 …
       << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
       << endl;
+  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
+      != TesselStruct->TrianglesOnBoundary.end(); runner++)
+  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
     { // go through every triangle, calculate volume of its pyramid with CoG as peak
+      x.CopyVector(&runner->second->endpoints[0]->node->x);
+      x.SubtractVector(&runner->second->endpoints[1]->node->x);
+      y.CopyVector(&runner->second->endpoints[0]->node->x);
+      y.SubtractVector(&runner->second->endpoints[2]->node->x);
+      a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
+          &runner->second->endpoints[1]->node->x));
+      b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
+          &runner->second->endpoints[2]->node->x));
+      c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
+          &runner->second->endpoints[1]->node->x));
+      x.CopyVector(runner->second->endpoints[0]->node->node);
+      x.SubtractVector(runner->second->endpoints[1]->node->node);
+      y.CopyVector(runner->second->endpoints[0]->node->node);
+      y.SubtractVector(runner->second->endpoints[2]->node->node);
+      a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
+      b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
+      c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
       G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
+      x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
+          &runner->second->endpoints[1]->node->x,
+          &runner->second->endpoints[2]->node->x);
+      x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
+      x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
+      x.Scale(runner->second->endpoints[1]->node->node->Projection(&x));
       h = x.Norm(); // distance of CoG to triangle
       PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
 …
         FillIt = true;
         for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
           //FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
+          FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
+        }
 …
-// =========================================================== class TESSELATION ===========================================
-/** Constructor of class Tesselation.
- */
-Tesselation::Tesselation()
+{
-  PointsOnBoundaryCount = 0;
-  LinesOnBoundaryCount = 0;
-  TrianglesOnBoundaryCount = 0;
-  TriangleFilesWritten = 0;
+}
+;
-/** Constructor of class Tesselation.
- * We have to free all points, lines and triangles.
- */
-Tesselation::~Tesselation()
+{
-  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
-  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
-    if (runner->second != NULL) {
-      delete (runner->second);
-      runner->second = NULL;
-    } else
-      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+  }
+}
+;
-/** Gueses first starting triangle of the convex envelope.
- * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
- * \param *out output stream for debugging
- * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
- */
-void
-Tesselation::GuessStartingTriangle(ofstream *out)
+{
-  // 4b. create a starting triangle
-  // 4b1. create all distances
-  DistanceMultiMap DistanceMMap;
-  double distance, tmp;
-  Vector PlaneVector, TrialVector;
-  PointMap::iterator A, B, C; // three nodes of the first triangle
-  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
-  // with A chosen, take each pair B,C and sort
-  if (A != PointsOnBoundary.end())
+    {
-      B = A;
-      B++;
-      for (; B != PointsOnBoundary.end(); B++)
+        {
-          C = B;
-          C++;
-          for (; C != PointsOnBoundary.end(); C++)
+            {
-              tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
-              distance = tmp * tmp;
-              tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
-              distance += tmp * tmp;
-              tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
-              distance += tmp * tmp;
-              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
-                  PointMap::iterator, PointMap::iterator> (B, C)));
+            }
+        }
+    }
-  //    // listing distances
-  //    *out << Verbose(1) << "Listing DistanceMMap:";
-  //    for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
-  //      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
-  //    }
-  //    *out << endl;
-  // 4b2. pick three baselines forming a triangle
-  // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
-  DistanceMultiMap::iterator baseline = DistanceMMap.begin();
-  for (; baseline != DistanceMMap.end(); baseline++)
+    {
-      // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
-      // 2. next, we have to check whether all points reside on only one side of the triangle
-      // 3. construct plane vector
-      PlaneVector.MakeNormalVector(&A->second->node->x,
-          &baseline->second.first->second->node->x,
-          &baseline->second.second->second->node->x);
-      *out << Verbose(2) << "Plane vector of candidate triangle is ";
-      PlaneVector.Output(out);
-      *out << endl;
-      // 4. loop over all points
-      double sign = 0.;
-      PointMap::iterator checker = PointsOnBoundary.begin();
-      for (; checker != PointsOnBoundary.end(); checker++)
+        {
-          // (neglecting A,B,C)
-          if ((checker == A) || (checker == baseline->second.first) || (checker
-              == baseline->second.second))
-            continue;
-          // 4a. project onto plane vector
-          TrialVector.CopyVector(&checker->second->node->x);
-          TrialVector.SubtractVector(&A->second->node->x);
-          distance = TrialVector.Projection(&PlaneVector);
-          if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
-            continue;
-          *out << Verbose(3) << "Projection of " << checker->second->node->Name
-              << " yields distance of " << distance << "." << endl;
-          tmp = distance / fabs(distance);
-          // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
-          if ((sign != 0) && (tmp != sign))
+            {
-              // 4c. If so, break 4. loop and continue with next candidate in 1. loop
-              *out << Verbose(2) << "Current candidates: "
-                  << A->second->node->Name << ","
-                  << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leaves "
-                  << checker->second->node->Name << " outside the convex hull."
-                  << endl;
-              break;
+            }
-          else
-            { // note the sign for later
-              *out << Verbose(2) << "Current candidates: "
-                  << A->second->node->Name << ","
-                  << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leave "
-                  << checker->second->node->Name << " inside the convex hull."
-                  << endl;
-              sign = tmp;
+            }
-          // 4d. Check whether the point is inside the triangle (check distance to each node
-          tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
-          int innerpoint = 0;
-          if ((tmp < A->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x)) && (tmp
-              < A->second->node->x.DistanceSquared(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x);
-          if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
-              &A->second->node->x)) && (tmp
-              < baseline->second.first->second->node->x.DistanceSquared(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.DistanceSquared(
-              &baseline->second.second->second->node->x);
-          if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x)) && (tmp
-              < baseline->second.second->second->node->x.DistanceSquared(
-                  &A->second->node->x)))
-            innerpoint++;
-          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
-          if (innerpoint == 3)
-            break;
+        }
-      // 5. come this far, all on same side? Then break 1. loop and construct triangle
-      if (checker == PointsOnBoundary.end())
+        {
-          *out << "Looks like we have a candidate!" << endl;
-          break;
+        }
+    }
-  if (baseline != DistanceMMap.end())
+    {
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = A->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = A->second;
-      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      // 4b3. insert created triangle
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-      TrianglesOnBoundaryCount++;
-      for (int i = 0; i < NDIM; i++)
+        {
-          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
-          LinesOnBoundaryCount++;
+        }
-      *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+    }
-  else
+    {
-      *out << Verbose(1) << "No starting triangle found." << endl;
-      exit(255);
+    }
+}
+;
-/** Tesselates the convex envelope of a cluster from a single starting triangle.
- * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
- * 2 triangles. Hence, we go through all current lines:
- * -# if the lines contains to only one triangle
- * -# We search all points in the boundary
- *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
- *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
- *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
- *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
- * \param *out output stream for debugging
- * \param *configuration for IsAngstroem
- * \param *mol the cluster as a molecule structure
- */
-void Tesselation::TesselateOnBoundary(ofstream *out, molecule *mol)
+{
-  bool flag;
-  PointMap::iterator winner;
-  class BoundaryPointSet *peak = NULL;
-  double SmallestAngle, TempAngle;
-  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *MolCenter = NULL;
-  LineMap::iterator LineChecker[2];
-  MolCenter = mol->DetermineCenterOfAll(out);
-  // create a first tesselation with the given BoundaryPoints
-  do {
-    flag = false;
-    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
-      if (baseline->second->TrianglesCount == 1) {
-        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
-        SmallestAngle = M_PI;
-        // get peak point with respect to this base line's only triangle
-        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
-        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
-        for (int i = 0; i < 3; i++)
-          if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
-            peak = BTS->endpoints[i];
-        *out << Verbose(3) << " and has peak " << *peak << "." << endl;
-        // prepare some auxiliary vectors
-        Vector BaseLineCenter, BaseLine;
-        BaseLineCenter.CopyVector(&baseline->second->endpoints[0]->node->x);
-        BaseLineCenter.AddVector(&baseline->second->endpoints[1]->node->x);
-        BaseLineCenter.Scale(1. / 2.); // points now to center of base line
-        BaseLine.CopyVector(&baseline->second->endpoints[0]->node->x);
-        BaseLine.SubtractVector(&baseline->second->endpoints[1]->node->x);
-        // offset to center of triangle
-        CenterVector.Zero();
-        for (int i = 0; i < 3; i++)
-          CenterVector.AddVector(&BTS->endpoints[i]->node->x);
-        CenterVector.Scale(1. / 3.);
-        *out << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl;
-        // normal vector of triangle
-        NormalVector.CopyVector(MolCenter);
-        NormalVector.SubtractVector(&CenterVector);
-        BTS->GetNormalVector(NormalVector);
-        NormalVector.CopyVector(&BTS->NormalVector);
-        *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
-        // vector in propagation direction (out of triangle)
-        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
-        PropagationVector.MakeNormalVector(&BaseLine, &NormalVector);
-        TempVector.CopyVector(&CenterVector);
-        TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
-        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
-        if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
-          PropagationVector.Scale(-1.);
-        *out << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl;
-        winner = PointsOnBoundary.end();
-        // loop over all points and calculate angle between normal vector of new and present triangle
-        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
-          if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
-            *out << Verbose(3) << "Target point is " << *(target->second) << ":" << endl;
-            // first check direction, so that triangles don't intersect
-            VirtualNormalVector.CopyVector(&target->second->node->x);
-            VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
-            VirtualNormalVector.ProjectOntoPlane(&NormalVector);
-            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
-            *out << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
-            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
-              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
-              continue;
-            } else
-              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
-            // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
-            LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
-            LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
-            if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->TrianglesCount == 2))) {
-              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
-              continue;
+            }
-            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) {
-              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
-              continue;
+            }
-            // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
-            if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
-              *out << Verbose(4) << "Current target is peak!" << endl;
-              continue;
+            }
-            // check for linear dependence
-            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            TempVector.SubtractVector(&target->second->node->x);
-            helper.CopyVector(&baseline->second->endpoints[1]->node->x);
-            helper.SubtractVector(&target->second->node->x);
-            helper.ProjectOntoPlane(&TempVector);
-            if (fabs(helper.NormSquared()) < MYEPSILON) {
-              *out << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
-              continue;
+            }
-            // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
-            flag = true;
-            VirtualNormalVector.MakeNormalVector(&baseline->second->endpoints[0]->node->x, &baseline->second->endpoints[1]->node->x, &target->second->node->x);
-            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            TempVector.AddVector(&baseline->second->endpoints[1]->node->x);
-            TempVector.AddVector(&target->second->node->x);
-            TempVector.Scale(1./3.);
-            TempVector.SubtractVector(MolCenter);
-            // make it always point outward
-            if (VirtualNormalVector.Projection(&TempVector) < 0)
-              VirtualNormalVector.Scale(-1.);
-            // calculate angle
-            TempAngle = NormalVector.Angle(&VirtualNormalVector);
-            *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
-            if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
-              SmallestAngle = TempAngle;
-              winner = target;
-              *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
-            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
-              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
-              helper.CopyVector(&target->second->node->x);
-              helper.SubtractVector(&BaseLineCenter);
-              helper.ProjectOntoPlane(&BaseLine);
-              // ...the one with the smaller angle is the better candidate
-              TempVector.CopyVector(&target->second->node->x);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              TempAngle = TempVector.Angle(&helper);
-              TempVector.CopyVector(&winner->second->node->x);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              if (TempAngle < TempVector.Angle(&helper)) {
-                TempAngle = NormalVector.Angle(&VirtualNormalVector);
-                SmallestAngle = TempAngle;
-                winner = target;
-                *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
-              } else
-                *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
-            } else
-              *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
+          }
-        } // end of loop over all boundary points
-        // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
-        if (winner != PointsOnBoundary.end()) {
-          *out << Verbose(2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;
-          // create the lins of not yet present
-          BLS[0] = baseline->second;
-          // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
-          LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first);
-          LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first);
-          if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create
-            BPS[0] = baseline->second->endpoints[0];
-            BPS[1] = winner->second;
-            BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[1]));
-            LinesOnBoundaryCount++;
-          } else
-            BLS[1] = LineChecker[0]->second;
-          if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create
-            BPS[0] = baseline->second->endpoints[1];
-            BPS[1] = winner->second;
-            BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[2]));
-            LinesOnBoundaryCount++;
-          } else
-            BLS[2] = LineChecker[1]->second;
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-          TrianglesOnBoundaryCount++;
-        } else {
-          *out << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
+        }
-        // 5d. If the set of lines is not yet empty, go to 5. and continue
-      } else
-        *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl;
-  } while (flag);
-  // exit
-  delete(MolCenter);
-};
-/** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
- * \param *out output stream for debugging
- * \param *mol molecule structure with atoms
- * \return true - all straddling points insert, false - something went wrong
- */
-bool Tesselation::InsertStraddlingPoints(ofstream *out, molecule *mol)
+{
-  Vector Intersection;
-  atom *Walker = mol->start;
-  Vector *MolCenter = mol->DetermineCenterOfAll(out);
-  while (Walker->next != mol->end) {  // we only have to go once through all points, as boundary can become only bigger
-    // get the next triangle
-    BTS = FindClosestTriangleToPoint(out, &Walker->x);
-    if (BTS == NULL) {
-      *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
-      return false;
+    }
-    // get the intersection point
-    if (BTS->GetIntersectionInsideTriangle(out, MolCenter, &Walker->x, &Intersection)) {
-      // we have the intersection, check whether in- or outside of boundary
-      if ((MolCenter->DistanceSquared(&Walker->x) - MolCenter->DistanceSquared(&Intersection)) < -MYEPSILON) {
-        // inside, next!
-        *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
-      } else {
-        // outside!
-        *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
-        class BoundaryLineSet *OldLines[3], *NewLines[3];
-        class BoundaryPointSet *OldPoints[3], *NewPoint;
-        // store the three old lines and old points
-        for (int i=0;i<3;i++) {
-          OldLines[i] = BTS->lines[i];
-          OldPoints[i] = BTS->endpoints[i];
+        }
-        // add Walker to boundary points
-        AddPoint(Walker);
-        if (BPS[0] == NULL)
-          NewPoint = BPS[0];
-        else
-          continue;
-        // remove triangle
-        TrianglesOnBoundary.erase(BTS->Nr);
-        // create three new boundary lines
-        for (int i=0;i<3;i++) {
-          BPS[0] = NewPoint;
-          BPS[1] = OldPoints[i];
-          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
-          LinesOnBoundaryCount++;
+        }
-        // create three new triangle with new point
-        for (int i=0;i<3;i++) { // find all baselines
-          BLS[0] = OldLines[i];
-          int n = 1;
-          for (int j=0;j<3;j++) {
-            if (NewLines[j]->IsConnectedTo(BLS[0])) {
-              if (n>2) {
-                *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
-                return false;
-              } else
-                BLS[n++] = NewLines[j];
+            }
+          }
-          // create the triangle
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-          TrianglesOnBoundaryCount++;
+        }
+      }
-    } else { // something is wrong with FindClosestTriangleToPoint!
-      *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
-      return false;
+    }
+  }
-  // exit
-  delete(MolCenter);
-  return true;
-};
-/** Goes over all baselines and checks whether adjacent triangles and convex to each other.
- * \param *out output stream for debugging
- * \return true - all baselines were corrected, false - there are still concave pieces
- */
-bool Tesselation::CorrectConcaveBaselines(ofstream *out)
+{
-  class BoundaryTriangleSet *triangle[2];
-  class BoundaryLineSet *OldLines[4], *NewLine;
-  class BoundaryPointSet *OldPoints[2];
-  Vector BaseLineNormal;
-  class BoundaryLineSet *Base = NULL;
-  int OldTriangles[2], OldBaseLine;
-  int i;
-  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
-    Base = baseline->second;
-    // check convexity
-    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
-      *out << Verbose(3) << Base << " has two convex triangles." << endl;
-    } else { // not convex!
-      // get the two triangles
-      i=0;
-      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-        triangle[i++] = runner->second;
-      // gather four endpoints and four lines
-      for (int j=0;j<4;j++)
-        OldLines[j] = NULL;
-      for (int j=0;j<2;j++)
-        OldPoints[j] = NULL;
-      i=0;
-      for (int m=0;m<2;m++) { // go over both triangles
-        for (int j=0;j<3;j++) { // all of their endpoints and baselines
-          if (triangle[m]->lines[j] != Base) // pick not the central baseline
-            OldLines[i++] = triangle[m]->lines[j];
-          if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
-            OldPoints[m] = triangle[m]->endpoints[j];
+        }
+      }
-      if (i<4) {
-        *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-        return false;
+      }
-      for (int j=0;j<4;j++)
-        if (OldLines[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-          return false;
+        }
-      for (int j=0;j<2;j++)
-        if (OldPoints[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
-          return false;
+        }
-      // remove triangles
-      for (int j=0;j<2;j++) {
-        OldTriangles[j] = triangle[j]->Nr;
-        TrianglesOnBoundary.erase(OldTriangles[j]);
-        triangle[j] = NULL;
+      }
-      // remove baseline
-      OldBaseLine = Base->Nr;
-      LinesOnBoundary.erase(OldBaseLine);
-      Base = NULL;
-      // construct new baseline (with same number as old one)
-      BPS[0] = OldPoints[0];
-      BPS[1] = OldPoints[1];
-      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
-      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
-      // construct new triangles with flipped baseline
-      i=-1;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
-        i=2;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
-        i=3;
-      if (i!=-1) {
-        BLS[0] = OldLines[0];
-        BLS[1] = OldLines[i];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
-        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
-        BLS[1] = OldLines[1];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
-      } else {
-        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
-        return false;
+      }
+    }
+  }
-  return true;
-};
-/** States whether point is in- or outside of a tesselated surface.
- * \param *pointer point to be checked
- * \return true - is inside, false - is outside
- */
-bool Tesselation::IsInside(Vector *pointer)
+{
-  // hier kommt dann Saskias Routine hin...
-  return true;
-};
-/** Finds the closest triangle to a given point.
- * \param *out output stream for debugging
- * \param *x second endpoint of line
- * \return pointer triangle that is closest, NULL if none was found
- */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
+{
-  class BoundaryTriangleSet *triangle = NULL;
-  // hier kommt dann Saskias Routine hin...
-  return triangle;
-};
-/** Adds an atom to the tesselation::PointsOnBoundary list.
- * \param *Walker atom to add
- */
-void
-Tesselation::AddPoint(atom *Walker)
+{
-  PointTestPair InsertUnique;
-  BPS[0] = new class BoundaryPointSet(Walker);
-  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
-  if (InsertUnique.second) // if new point was not present before, increase counter
-    PointsOnBoundaryCount++;
-  else {
-    delete(BPS[0]);
-    BPS[0] = NULL;
+  }
+}
+;
-/** Adds point to Tesselation::PointsOnBoundary if not yet present.
- * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
- * @param Candidate point to add
- * @param n index for this point in Tesselation::TPS array
- */
-void
-Tesselation::AddTrianglePoint(atom* Candidate, int n)
+{
-  PointTestPair InsertUnique;
-  TPS[n] = new class BoundaryPointSet(Candidate);
-  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
-  if (InsertUnique.second) { // if new point was not present before, increase counter
-    PointsOnBoundaryCount++;
-  } else {
-    delete TPS[n];
-    cout << Verbose(3) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
-    TPS[n] = (InsertUnique.first)->second;
+  }
+}
+;
-/** Function tries to add line from current Points in BPS to BoundaryLineSet.
- * If successful it raises the line count and inserts the new line into the BLS,
- * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
- * @param *a first endpoint
- * @param *b second endpoint
- * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
-  bool insertNewLine = true;
-  if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine;
-    pair<LineMap::iterator,LineMap::iterator> FindPair;
-    FindPair = a->lines.equal_range(b->node->nr);
-    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-      // If there is a line with less than two attached triangles, we don't need a new line.
-      if (FindLine->second->TrianglesCount < 2) {
-        insertNewLine = false;
-        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
-        BPS[0] = FindLine->second->endpoints[0];
-        BPS[1] = FindLine->second->endpoints[1];
-        BLS[n] = FindLine->second;
-        break;
+      }
+    }
+  }
-  if (insertNewLine) {
-    AlwaysAddTriangleLine(a, b, n);
+  }
+}
+;
-/**
- * Adds lines from each of the current points in the BPS to BoundaryLineSet.
- * Raises the line count and inserts the new line into the BLS.
+ *
- * @param *a first endpoint
- * @param *b second endpoint
- * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
+{
-  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
-  BPS[0] = a;
-  BPS[1] = b;
-  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
-  // add line to global map
-  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
-  // increase counter
-  LinesOnBoundaryCount++;
-};
-/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
- * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
- */
-void
-Tesselation::AddTriangle()
+{
-  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
-  // add triangle to global map
-  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-  TrianglesOnBoundaryCount++;
-  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
+}
+;
-double det_get(gsl_matrix *A, int inPlace) {
-  /*
-  inPlace = 1 => A is replaced with the LU decomposed copy.
-  inPlace = 0 => A is retained, and a copy is used for LU.
-  */
-  double det;
-  int signum;
-  gsl_permutation *p = gsl_permutation_alloc(A->size1);
-  gsl_matrix *tmpA;
-  if (inPlace)
-  tmpA = A;
-  else {
-  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
-  gsl_matrix_memcpy(tmpA , A);
+  }
-  gsl_linalg_LU_decomp(tmpA , p , &signum);
-  det = gsl_linalg_LU_det(tmpA , signum);
-  gsl_permutation_free(p);
-  if (! inPlace)
-  gsl_matrix_free(tmpA);
-  return det;
-};
-void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
+{
-  gsl_matrix *A = gsl_matrix_calloc(3,3);
-  double m11, m12, m13, m14;
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
-  m11 = det_get(A, 1);
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
-  m12 = det_get(A, 1);
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
-  m13 = det_get(A, 1);
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, b.x[i]);
+  }
-  m14 = det_get(A, 1);
-  if (fabs(m11) < MYEPSILON)
-    cerr << "ERROR: three points are colinear." << endl;
-  center->x[0] =  0.5 * m12/ m11;
-  center->x[1] = -0.5 * m13/ m11;
-  center->x[2] =  0.5 * m14/ m11;
-  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
-    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
-  gsl_matrix_free(A);
-};
-/**
- * Function returns center of sphere with RADIUS, which rests on points a, b, c
- * @param Center this vector will be used for return
- * @param a vector first point of triangle
- * @param b vector second point of triangle
- * @param c vector third point of triangle
- * @param *Umkreismittelpunkt new cneter point of circumference
- * @param Direction vector indicates up/down
- * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
- * @param Halfplaneindicator double indicates whether Direction is up or down
- * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
- * @param alpha double angle at a
- * @param beta double, angle at b
- * @param gamma, double, angle at c
- * @param Radius, double
- * @param Umkreisradius double radius of circumscribing circle
- */
-void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
-    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+{
-  Vector TempNormal, helper;
-  double Restradius;
-  Vector OtherCenter;
-  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
-  Center->Zero();
-  helper.CopyVector(&a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(&b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(&c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
-  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
-  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-  NewUmkreismittelpunkt->CopyVector(Center);
-  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
-  // Here we calculated center of circumscribing circle, using barycentric coordinates
-  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
-  TempNormal.CopyVector(&a);
-  TempNormal.SubtractVector(&b);
-  helper.CopyVector(&a);
-  helper.SubtractVector(&c);
-  TempNormal.VectorProduct(&helper);
-  if (fabs(HalfplaneIndicator) < MYEPSILON)
+    {
-      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+        {
-          TempNormal.Scale(-1);
+        }
+    }
-  else
+    {
-      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+        {
-          TempNormal.Scale(-1);
+        }
+    }
-  TempNormal.Normalize();
-  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
-  TempNormal.Scale(Restradius);
-  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
-  Center->AddVector(&TempNormal);
-  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
-  get_sphere(&OtherCenter, a, b, c, RADIUS);
-  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
-  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
-};
-/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
- * \param *Center new center on return
- * \param *a first point
- * \param *b second point
- * \param *c third point
- */
-void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
+{
-  Vector helper;
-  double alpha, beta, gamma;
-  Vector SideA, SideB, SideC;
-  SideA.CopyVector(b);
-  SideA.SubtractVector(c);
-  SideB.CopyVector(c);
-  SideB.SubtractVector(a);
-  SideC.CopyVector(a);
-  SideC.SubtractVector(b);
-  alpha = M_PI - SideB.Angle(&SideC);
-  beta = M_PI - SideC.Angle(&SideA);
-  gamma = M_PI - SideA.Angle(&SideB);
-  //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
-  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
-    cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
-  Center->Zero();
-  helper.CopyVector(a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
-  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-};
-/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
- * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
- * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
- * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
- * \param CircleCenter Center of the parameter circle
- * \param CirclePlaneNormal normal vector to plane of the parameter circle
- * \param CircleRadius radius of the parameter circle
- * \param NewSphereCenter new center of a circumcircle
- * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
- * \param NormalVector normal vector
- * \param SearchDirection search direction to make angle unique on return.
- * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
- */
-double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
+{
-  Vector helper;
-  double radius, alpha;
-  helper.CopyVector(&NewSphereCenter);
-  // test whether new center is on the parameter circle's plane
-  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
-    helper.ProjectOntoPlane(&CirclePlaneNormal);
+  }
-  radius = helper.ScalarProduct(&helper);
-  // test whether the new center vector has length of CircleRadius
-  if (fabs(radius - CircleRadius) > HULLEPSILON)
-    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-  alpha = helper.Angle(&OldSphereCenter);
-  // make the angle unique by checking the halfplanes/search direction
-  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
-    alpha = 2.*M_PI - alpha;
-  //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
-  radius = helper.Distance(&OldSphereCenter);
-  helper.ProjectOntoPlane(&NormalVector);
-  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
-  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
-    //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
-    return alpha;
-  } else {
-    //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
-    return 2.*M_PI;
+  }
-};
-/** Checks whether the triangle consisting of the three atoms is already present.
- * Searches for the points in Tesselation::PointsOnBoundary and checks their
- * lines. If any of the three edges already has two triangles attached, false is
- * returned.
- * \param *out output stream for debugging
- * \param *Candidates endpoints of the triangle candidate
- * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
- *                 triangles exist which is the maximum for three points
- */
-int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
-// TODO: use findTriangles here and return result.size();
-  LineMap::iterator FindLine;
-  PointMap::iterator FindPoint;
-  TriangleMap::iterator FindTriangle;
-  int adjacentTriangleCount = 0;
-  class BoundaryPointSet *Points[3];
-  //*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
-  // builds a triangle point set (Points) of the end points
-  for (int i = 0; i < 3; i++) {
-    FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      Points[i] = FindPoint->second;
-    } else {
-      Points[i] = NULL;
+    }
+  }
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (Points[i] != NULL) {
-      for (int j = i; j < 3; j++) {
-        if (Points[j] != NULL) {
-          FindLine = Points[i]->lines.find(Points[j]->node->nr);
-          if (FindLine != Points[i]->lines.end()) {
-            for (; FindLine->first == Points[j]->node->nr; FindLine++) {
-              FindTriangle = FindLine->second->triangles.begin();
-              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if ((
-                  (FindTriangle->second->endpoints[0] == Points[0])
-                    || (FindTriangle->second->endpoints[0] == Points[1])
-                    || (FindTriangle->second->endpoints[0] == Points[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[1] == Points[0])
-                    || (FindTriangle->second->endpoints[1] == Points[1])
-                    || (FindTriangle->second->endpoints[1] == Points[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[2] == Points[0])
-                    || (FindTriangle->second->endpoints[2] == Points[1])
-                    || (FindTriangle->second->endpoints[2] == Points[2])
+                  )
-                ) {
-                  adjacentTriangleCount++;
+                }
+              }
+            }
-            // Only one of the triangle lines must be considered for the triangle count.
-            *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-            return adjacentTriangleCount;
+          }
+        }
+      }
+    }
+  }
-  *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-  return adjacentTriangleCount;
-};
-/** This recursive function finds a third point, to form a triangle with two given ones.
- * Note that this function is for the starting triangle.
- * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
- * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
- * the center of the sphere is still fixed up to a single parameter. The band of possible values
- * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
- * us the "null" on this circle, the new center of the candidate point will be some way along this
- * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
- * by the normal vector of the base triangle that always points outwards by construction.
- * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
- * We construct the normal vector that defines the plane this circle lies in, it is just in the
- * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
- * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
- * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
- * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
- * both.
- * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
- * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
- * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
- * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
- * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
- * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
- * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
- * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
- * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
- * @param BaseLine BoundaryLineSet with the current base line
- * @param ThirdNode third atom to avoid in search
- * @param candidates list of equally good candidates to return
- * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
- * @param RADIUS radius of sphere
- * @param *LC LinkedCell structure with neighbouring atoms
- */
-void Find_third_point_for_Tesselation(
-  Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter,
-  class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates,
-  double *ShortestAngle, const double RADIUS, LinkedCell *LC
-) {
-  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
-  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
-  Vector SphereCenter;
-  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
-  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
-  Vector NewNormalVector;   // normal vector of the Candidate's triangle
-  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
-  LinkedAtoms *List = NULL;
-  double CircleRadius; // radius of this circle
-  double radius;
-  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-  atom *Candidate = NULL;
-  CandidateForTesselation *optCandidate = NULL;
-  cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
-  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
-  // construct center of circle
-  CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
-  CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
-  CircleCenter.Scale(0.5);
-  // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
-  CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
-  // calculate squared radius atom *ThirdNode,f circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
-    CirclePlaneNormal.Normalize();
-    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
-    // test whether old center is on the band's plane
-    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
-      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
-    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
-    if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      // check SearchDirection
-      //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
-        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+      }
-      // get cell for the starting atom
-      if (LC->SetIndexToVector(&CircleCenter)) {
-        for(int i=0;i<NDIM;i++) // store indices of this cell
-        N[i] = LC->n[i];
-        //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-      } else {
-        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
-        return;
+      }
-      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-      //cout << Verbose(2) << "LC Intervals:";
-      for (int i=0;i<NDIM;i++) {
-        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-        //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+      }
-      //cout << endl;
-      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-            List = LC->GetCurrentCell();
-            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-            if (List != NULL) {
-              for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-                Candidate = (*Runner);
-                // check for three unique points
-                //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
-                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){
-                  // construct both new centers
-                  GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
-                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
-                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
-                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
-                  ) {
-                    helper.CopyVector(&NewNormalVector);
-                    //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
-                    if (radius < RADIUS*RADIUS) {
-                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
-                      //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
-                      NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
-                      //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-                      // OtherNewSphereCenter is created by the same vector just in the other direction
-                      helper.Scale(-1.);
-                      OtherNewSphereCenter.AddVector(&helper);
-                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
-                      //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
-                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                      alpha = min(alpha, Otheralpha);
-                      // if there is a better candidate, drop the current list and add the new candidate
-                      // otherwise ignore the new candidate and keep the list
-                      if (*ShortestAngle > (alpha - HULLEPSILON)) {
-                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
-                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
-                        } else {
-                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
+                        }
-                        // if there is an equal candidate, add it to the list without clearing the list
-                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
-                          candidates->push_back(optCandidate);
-                          cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
-                        } else {
-                          // remove all candidates from the list and then the list itself
-                          class CandidateForTesselation *remover = NULL;
-                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
-                            remover = *it;
-                            delete(remover);
+                          }
-                          candidates->clear();
-                          candidates->push_back(optCandidate);
-                          cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
+                        }
-                        *ShortestAngle = alpha;
-                        //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
-                      } else {
-                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
-                          //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
-                        } else {
-                          //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                        }
+                      }
-                    } else {
-                      //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+                    }
-                  } else {
-                    //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+                  }
-                } else {
-                  if (ThirdNode != NULL) {
-                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
-                  } else {
-                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                  }
+                }
+              }
+            }
+          }
-    } else {
-      cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+    }
-  } else {
-    if (ThirdNode != NULL)
-      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
-    else
-      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
+  }
-  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
-  if (candidates->size() > 1) {
-    candidates->unique();
-    candidates->sort(sortCandidates);
+  }
-  cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
-};
-struct Intersection {
-  Vector x1;
-  Vector x2;
-  Vector x3;
-  Vector x4;
-};
-/**
- * Intersection calculation function.
+ *
- * @param x to find the result for
- * @param function parameter
- */
-double MinIntersectDistance(const gsl_vector * x, void *params) {
-  double retval = 0;
-  struct Intersection *I = (struct Intersection *)params;
-  Vector intersection;
-  Vector SideA,SideB,HeightA, HeightB;
-  for (int i=0;i<NDIM;i++)
-    intersection.x[i] = gsl_vector_get(x, i);
-  SideA.CopyVector(&(I->x1));
-  SideA.SubtractVector(&I->x2);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&I->x1);
-  HeightA.ProjectOntoPlane(&SideA);
-  SideB.CopyVector(&I->x3);
-  SideB.SubtractVector(&I->x4);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&I->x3);
-  HeightB.ProjectOntoPlane(&SideB);
-  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
-  //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
-  return retval;
-};
-/**
- * Calculates whether there is an intersection between two lines. The first line
- * always goes through point 1 and point 2 and the second line is given by the
- * connection between point 4 and point 5.
+ *
- * @param point 1 of line 1
- * @param point 2 of line 1
- * @param point 1 of line 2
- * @param point 2 of line 2
+ *
- * @return true if there is an intersection between the given lines, false otherwise
- */
-bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) {
-  bool result;
-  struct Intersection par;
-    par.x1.CopyVector(&point1);
-    par.x2.CopyVector(&point2);
-    par.x3.CopyVector(&point3);
-    par.x4.CopyVector(&point4);
-    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
-    gsl_multimin_fminimizer *s = NULL;
-    gsl_vector *ss, *x;
-    gsl_multimin_function minex_func;
-    size_t iter = 0;
-    int status;
-    double size;
-    /* Starting point */
-    x = gsl_vector_alloc(NDIM);
-    gsl_vector_set(x, 0, point1.x[0]);
-    gsl_vector_set(x, 1, point1.x[1]);
-    gsl_vector_set(x, 2, point1.x[2]);
-    /* Set initial step sizes to 1 */
-    ss = gsl_vector_alloc(NDIM);
-    gsl_vector_set_all(ss, 1.0);
-    /* Initialize method and iterate */
-    minex_func.n = NDIM;
-    minex_func.f = &MinIntersectDistance;
-    minex_func.params = (void *)&par;
-    s = gsl_multimin_fminimizer_alloc(T, NDIM);
-    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
-    do {
-        iter++;
-        status = gsl_multimin_fminimizer_iterate(s);
-        if (status) {
-          break;
+        }
-        size = gsl_multimin_fminimizer_size(s);
-        status = gsl_multimin_test_size(size, 1e-2);
-        if (status == GSL_SUCCESS) {
-          cout << Verbose(2) << "converged to minimum" <<  endl;
+        }
-    } while (status == GSL_CONTINUE && iter < 100);
-    // check whether intersection is in between or not
-  Vector intersection, SideA, SideB, HeightA, HeightB;
-  double t1, t2;
-  for (int i = 0; i < NDIM; i++) {
-    intersection.x[i] = gsl_vector_get(s->x, i);
+  }
-  SideA.CopyVector(&par.x2);
-  SideA.SubtractVector(&par.x1);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&par.x1);
-  t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);
-  SideB.CopyVector(&par.x4);
-  SideB.SubtractVector(&par.x3);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&par.x3);
-  t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);
-  cout << Verbose(2) << "Intersection " << intersection << " is at "
-    << t1 << " for (" << point1 << "," << point2 << ") and at "
-    << t2 << " for (" << point3 << "," << point4 << "): ";
-  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
-    cout << "true intersection." << endl;
-    result = true;
-  } else {
-    cout << "intersection out of region of interest." << endl;
-    result = false;
+  }
-  // free minimizer stuff
-    gsl_vector_free(x);
-    gsl_vector_free(ss);
-    gsl_multimin_fminimizer_free(s);
-  return result;
+}
-/** Finds the second point of starting triangle.
- * \param *a first atom
- * \param *Candidate pointer to candidate atom on return
- * \param Oben vector indicating the outside
- * \param Opt_Candidate reference to recommended candidate on return
- * \param Storage[3] array storing angles and other candidate information
- * \param RADIUS radius of virtual sphere
- * \param *LC LinkedCell structure with neighbouring atoms
- */
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
+{
-  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
-  Vector AngleCheck;
-  double norm = -1., angle;
-  LinkedAtoms *List = NULL;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-  if (LC->SetIndexToAtom(a)) {  // get cell for the starting atom
-    for(int i=0;i<NDIM;i++) // store indices of this cell
-      N[i] = LC->n[i];
-  } else {
-    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
-    return;
+  }
-  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-  cout << Verbose(3) << "LC Intervals from [";
-  for (int i=0;i<NDIM;i++) {
-  cout << " " << N[i] << "<->" << LC->N[i];
+  }
-  cout << "] :";
-  for (int i=0;i<NDIM;i++) {
-    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+  }
-  cout << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            Candidate = (*Runner);
-            // check if we only have one unique point yet ...
-            if (a != Candidate) {
-              // Calculate center of the circle with radius RADIUS through points a and Candidate
-              Vector OrthogonalizedOben, a_Candidate, Center;
-              double distance, scaleFactor;
-              OrthogonalizedOben.CopyVector(&Oben);
-              a_Candidate.CopyVector(&(a->x));
-              a_Candidate.SubtractVector(&(Candidate->x));
-              OrthogonalizedOben.ProjectOntoPlane(&a_Candidate);
-              OrthogonalizedOben.Normalize();
-              distance = 0.5 * a_Candidate.Norm();
-              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
-              OrthogonalizedOben.Scale(scaleFactor);
-              Center.CopyVector(&(Candidate->x));
-              Center.AddVector(&(a->x));
-              Center.Scale(0.5);
-              Center.AddVector(&OrthogonalizedOben);
-              AngleCheck.CopyVector(&Center);
-              AngleCheck.SubtractVector(&(a->x));
-              norm = a_Candidate.Norm();
-              // second point shall have smallest angle with respect to Oben vector
-              if (norm < RADIUS*2.) {
-                angle = AngleCheck.Angle(&Oben);
-                if (angle < Storage[0]) {
-                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
-                  Opt_Candidate = Candidate;
-                  Storage[0] = angle;
-                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
-                } else {
-                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+                }
-              } else {
-                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
+              }
-            } else {
-              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
+            }
+          }
-        } else {
-          cout << Verbose(3) << "Linked cell list is empty." << endl;
+        }
+      }
-  cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
-};
-/** Finds the starting triangle for find_non_convex_border().
- * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
- * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
- * point are called.
- * \param RADIUS radius of virtual rolling sphere
- * \param *LC LinkedCell structure with neighbouring atoms
- */
-void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
+{
-  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
-  int i = 0;
-  LinkedAtoms *List = NULL;
-  atom* FirstPoint = NULL;
-  atom* SecondPoint = NULL;
-  atom* MaxAtom[NDIM];
-  double max_coordinate[NDIM];
-  Vector Oben;
-  Vector helper;
-  Vector Chord;
-  Vector SearchDirection;
-  Oben.Zero();
-  for (i = 0; i < 3; i++) {
-    MaxAtom[i] = NULL;
-    max_coordinate[i] = -1;
+  }
-  // 1. searching topmost atom with respect to each axis
-  for (int i=0;i<NDIM;i++) { // each axis
-    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
-    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
-      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
-        List = LC->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
-            if ((*Runner)->x.x[i] > max_coordinate[i]) {
-              cout << Verbose(2) << "New maximal for axis " << i << " atom is " << *(*Runner) << " at " << (*Runner)->x << "." << endl;
-              max_coordinate[i] = (*Runner)->x.x[i];
-              MaxAtom[i] = (*Runner);
+            }
+          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  }
-  cout << Verbose(2) << "Found maximum coordinates: ";
-  for (int i=0;i<NDIM;i++)
-    cout << i << ": " << *MaxAtom[i] << "\t";
-  cout << endl;
-  BTS = NULL;
-  CandidateList *Opt_Candidates = new CandidateList();
-  for (int k=0;k<NDIM;k++) {
-    Oben.x[k] = 1.;
-    FirstPoint = MaxAtom[k];
-    cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
-    double ShortestAngle;
-    atom* Opt_Candidate = NULL;
-    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.
-    Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
-    SecondPoint = Opt_Candidate;
-    if (SecondPoint == NULL)  // have we found a second point?
-      continue;
-    else
-      cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
-    helper.CopyVector(&(FirstPoint->x));
-    helper.SubtractVector(&(SecondPoint->x));
-    helper.Normalize();
-    Oben.ProjectOntoPlane(&helper);
-    Oben.Normalize();
-    helper.VectorProduct(&Oben);
-    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-    Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-    Chord.SubtractVector(&(SecondPoint->x));
-    double radius = Chord.ScalarProduct(&Chord);
-    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
-    helper.CopyVector(&Oben);
-    helper.Scale(CircleRadius);
-    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
-    // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
-    // adding point 1 and point 2 and the line between them
-    AddTrianglePoint(FirstPoint, 0);
-    AddTrianglePoint(SecondPoint, 1);
-    AddTriangleLine(TPS[0], TPS[1], 0);
-    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    Find_third_point_for_Tesselation(
-      Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
-    );
-    cout << Verbose(1) << "List of third Points is ";
-    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-        cout << " " << *(*it)->point;
+    }
-    cout << endl;
-    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-      // add third triangle point
-      AddTrianglePoint((*it)->point, 2);
-      // add the second and third line
-      AddTriangleLine(TPS[1], TPS[2], 1);
-      AddTriangleLine(TPS[0], TPS[2], 2);
-      // ... and triangles to the Maps of the Tesselation class
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangle();
-      // ... and calculate its normal vector (with correct orientation)
-      (*it)->OptCenter.Scale(-1.);
-      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
-      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
-      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
-      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
-      // if we do not reach the end with the next step of iteration, we need to setup a new first line
-      if (it != Opt_Candidates->end()--) {
-        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
-        SecondPoint = (*it)->point;
-        // adding point 1 and point 2 and the line between them
-        AddTrianglePoint(FirstPoint, 0);
-        AddTrianglePoint(SecondPoint, 1);
-        AddTriangleLine(TPS[0], TPS[1], 0);
+      }
-      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
+    }
-    if (BTS != NULL) // we have created one starting triangle
-      break;
-    else {
-      // remove all candidates from the list and then the list itself
-      class CandidateForTesselation *remover = NULL;
-      for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-        remover = *it;
-        delete(remover);
+      }
-      Opt_Candidates->clear();
+    }
+  }
-  // remove all candidates from the list and then the list itself
-  class CandidateForTesselation *remover = NULL;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    remover = *it;
-    delete(remover);
+  }
-  delete(Opt_Candidates);
-  cout << Verbose(1) << "End of Find_starting_triangle\n";
-};
-/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
- * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
- * make it bigger (i.e. closing one (the baseline) and opening two new ones).
- * \param TPS[3] nodes of the triangle
- * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
- */
-bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3])
+{
-  bool result = false;
-  int counter = 0;
-  // check all three points
-  for (int i=0;i<3;i++)
-    for (int j=i+1; j<3; j++) {
-      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
-        LineMap::iterator FindLine;
-        pair<LineMap::iterator,LineMap::iterator> FindPair;
-        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
-        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-          // If there is a line with less than two attached triangles, we don't need a new line.
-          if (FindLine->second->TrianglesCount < 2) {
-            counter++;
-            break;  // increase counter only once per edge
+          }
+        }
-      } else { // no line
-        cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
-        result = true;
+      }
+    }
-  if (counter > 1) {
-    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
-    result = true;
+  }
-  return result;
-};
-/** This function finds a triangle to a line, adjacent to an existing one.
- * @param out output stream for debugging
- * @param *mol molecule with Atom's and Bond's
- * @param Line current baseline to search from
- * @param T current triangle which \a Line is edge of
- * @param RADIUS radius of the rolling ball
- * @param N number of found triangles
- * @param *filename filename base for intermediate envelopes
- * @param *LC LinkedCell structure with neighbouring atoms
- */
-bool Tesselation::Find_next_suitable_triangle(ofstream *out,
-    molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
-    const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
+{
-  cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n";
-  ofstream *tempstream = NULL;
-  char NumberName[255];
-  bool result = true;
-  CandidateList *Opt_Candidates = new CandidateList();
-  Vector CircleCenter;
-  Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
-  Vector SearchDirection;
-  Vector helper;
-  atom *ThirdNode = NULL;
-  LineMap::iterator testline;
-  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-  double radius, CircleRadius;
-  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-  for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
-      ThirdNode = T.endpoints[i]->node;
-  // construct center of circle
-  CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
-  CircleCenter.AddVector(&Line.endpoints[1]->node->x);
-  CircleCenter.Scale(0.5);
-  // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
-  CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
-  // calculate squared radius of circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
-    CirclePlaneNormal.Normalize();
-    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
-    // construct old center
-    GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
-    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
-    radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
-    helper.Scale(sqrt(RADIUS*RADIUS - radius));
-    OldSphereCenter.AddVector(&helper);
-    OldSphereCenter.SubtractVector(&CircleCenter);
-    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-    // construct SearchDirection
-    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(&Line.endpoints[0]->node->x);
-    helper.SubtractVector(&ThirdNode->x);
-    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
-      SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
-    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
-      // rotated the wrong way!
-      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+    }
-    // add third point
-    Find_third_point_for_Tesselation(
-      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates,
-      &ShortestAngle, RADIUS, LC
-    );
-  } else {
-    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+  }
-  if (Opt_Candidates->begin() == Opt_Candidates->end()) {
-    cerr << "WARNING: Could not find a suitable candidate." << endl;
-    return false;
+  }
-  cout << Verbose(1) << "Third Points are ";
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    cout << " " << *(*it)->point;
+  }
-  cout << endl;
-  BoundaryLineSet *BaseRay = &Line;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    cout << Verbose(1) << " Third point candidate is " << *(*it)->point
-    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
-    cout << Verbose(1) << " Baseline is " << *BaseRay << endl;
-    // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-    atom *AtomCandidates[3];
-    AtomCandidates[0] = (*it)->point;
-    AtomCandidates[1] = BaseRay->endpoints[0]->node;
-    AtomCandidates[2] = BaseRay->endpoints[1]->node;
-    int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);
-    BTS = NULL;
-    // If there is no triangle, add it regularly.
-    if (existentTrianglesCount == 0) {
-      AddTrianglePoint((*it)->point, 0);
-      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-      AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
-      AddTriangleLine(TPS[0], TPS[1], 0);
-      AddTriangleLine(TPS[0], TPS[2], 1);
-      AddTriangleLine(TPS[1], TPS[2], 2);
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangle();
-      (*it)->OptCenter.Scale(-1.);
-      BTS->GetNormalVector((*it)->OptCenter);
-      (*it)->OptCenter.Scale(-1.);
-      cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-        << " for this triangle ... " << endl;
-      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
-    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
-        AddTrianglePoint((*it)->point, 0);
-        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-        AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
-        // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
-        // i.e. at least one of the three lines must be present with TriangleCount <= 1
-        if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
-          AddTriangleLine(TPS[0], TPS[1], 0);
-          AddTriangleLine(TPS[0], TPS[2], 1);
-          AddTriangleLine(TPS[1], TPS[2], 2);
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          AddTriangle();
-          (*it)->OtherOptCenter.Scale(-1.);
-          BTS->GetNormalVector((*it)->OtherOptCenter);
-          (*it)->OtherOptCenter.Scale(-1.);
-          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-          << " for this triangle ... " << endl;
-          cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
-        } else {
-          cout << Verbose(1) << "WARNING: This triangle consisting of ";
-          cout << *(*it)->point << ", ";
-          cout << *BaseRay->endpoints[0]->node << " and ";
-          cout << *BaseRay->endpoints[1]->node << " ";
-          cout << "exists and is not added, as it does not seem helpful!" << endl;
-          result = false;
+        }
-    } else {
-      cout << Verbose(1) << "This triangle consisting of ";
-      cout << *(*it)->point << ", ";
-      cout << *BaseRay->endpoints[0]->node << " and ";
-      cout << *BaseRay->endpoints[1]->node << " ";
-      cout << "is invalid!" << endl;
-      result = false;
+    }
-    if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
-      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
-      if (DoTecplotOutput) {
-        string NameofTempFile(tempbasename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(TecplotSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
+      }
-      if (DoRaster3DOutput) {
-        string NameofTempFile(tempbasename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(Raster3DSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_raster3d_file(out, tempstream, this, mol);
-        // include the current position of the virtual sphere in the temporary raster3d file
-        // make the circumsphere's center absolute again
-        helper.CopyVector(&BaseRay->endpoints[0]->node->x);
-        helper.AddVector(&BaseRay->endpoints[1]->node->x);
-        helper.Scale(0.5);
-        (*it)->OptCenter.AddVector(&helper);
-        Vector *center = mol->DetermineCenterOfAll(out);
-        (*it)->OptCenter.SubtractVector(center);
-        delete(center);
-        // and add to file plus translucency object
-        *tempstream << "# current virtual sphere\n";
-        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
-          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
-          << "\t" << RADIUS << "\t1 0 0\n";
-        *tempstream << "9\n  terminating special property\n";
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
+      }
-      if (DoTecplotOutput || DoRaster3DOutput)
-        TriangleFilesWritten++;
+    }
-    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
-    BaseRay = BLS[0];
+  }
-  // remove all candidates from the list and then the list itself
-  class CandidateForTesselation *remover = NULL;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    remover = *it;
-    delete(remover);
+  }
-  delete(Opt_Candidates);
-  cout << Verbose(0) << "End of Find_next_suitable_triangle\n";
-  return result;
-};
-/**
- * Sort function for the candidate list.
- */
-bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
-  // TODO: use get angle and remove duplicate code
-  Vector BaseLineVector, OrthogonalVector, helper;
-        if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
-          cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
-          //return false;
-          exit(1);
+        }
-        // create baseline vector
-        BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x));
-        BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-        BaseLineVector.Normalize();
-  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
-  helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  helper.SubtractVector(&(candidate1->point->x));
-  OrthogonalVector.CopyVector(&helper);
-  helper.VectorProduct(&BaseLineVector);
-  OrthogonalVector.SubtractVector(&helper);
-  OrthogonalVector.Normalize();
-  // calculate both angles and correct with in-plane vector
-  helper.CopyVector(&(candidate1->point->x));
-  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  double phi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    phi = 2.*M_PI - phi;
+  }
-  helper.CopyVector(&(candidate2->point->x));
-  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  double psi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    psi = 2.*M_PI - psi;
+  }
-  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
-  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
-  // return comparison
-  return phi < psi;
+}
-/**
- * Checks whether the provided atom is within the tesselation structure.
+ *
- * @param atom of which to check the position
- * @param tesselation structure
+ *
- * @return true if the atom is inside the tesselation structure, false otherwise
- */
-bool IsInnerAtom(atom *Atom, class Tesselation *Tess, LinkedCell* LC) {
-  if (Tess->LinesOnBoundary.begin() == Tess->LinesOnBoundary.end()) {
-    cout << Verbose(0) << "Error: There is no tesselation structure to compare the atom with, "
-        << "please create one first.";
-    exit(1);
+  }
-  class atom *trianglePoints[3];
-  trianglePoints[0] = findClosestAtom(Atom, LC);
-  // check whether closest atom is "too close" :), then it's inside
-  if (trianglePoints[0]->x.DistanceSquared(&Atom->x) < MYEPSILON)
-    return true;
-  list<atom*> *connectedClosestAtoms = Tess->getClosestConnectedAtoms(trianglePoints[0], Atom);
-  trianglePoints[1] = connectedClosestAtoms->front();
-  trianglePoints[2] = connectedClosestAtoms->back();
-  for (int i=0;i<3;i++) {
-    if (trianglePoints[i] == NULL) {
-      cout << Verbose(1) << "IsInnerAtom encounters serious error, point " << i << " not found." << endl;
+    }
-    cout << Verbose(1) << "List of possible atoms:" << endl;
-    cout << *trianglePoints[i] << endl;
-//    for(list<atom*>::iterator runner = connectedClosestAtoms->begin(); runner != connectedClosestAtoms->end(); runner++)
-//      cout << Verbose(2) << *(*runner) << endl;
+  }
-  delete(connectedClosestAtoms);
-  list<BoundaryTriangleSet*> *triangles = Tess->FindTriangles(trianglePoints);
-  if (triangles->empty()) {
-    cout << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
-    exit(1);
+  }
-  Vector helper;
-  helper.CopyVector(&Atom->x);
-  // Only in case of degeneration, there will be two different scalar products.
-  // If at least one scalar product is positive, the point is considered to be outside.
-  bool status = (helper.ScalarProduct(&triangles->front()->NormalVector) < 0)
-      && (helper.ScalarProduct(&triangles->back()->NormalVector) < 0);
-  delete(triangles);
-  return status;
+}
-/**
- * Finds the atom which is closest to the provided one.
+ *
- * @param atom to which to find the closest other atom
- * @param linked cell structure
+ *
- * @return atom which is closest to the provided one
- */
-atom* findClosestAtom(const atom* Atom, LinkedCell* LC) {
-  LinkedAtoms *List = NULL;
-  atom* closestAtom = NULL;
-  double distance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-  LC->SetIndexToVector(&Atom->x); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //cout << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        //cout << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(&Atom->x);
-            helper.SubtractVector(&(*Runner)->x);
-            double currentNorm = helper. Norm();
-            if (currentNorm < distance) {
-              distance = currentNorm;
-              closestAtom = (*Runner);
+            }
+          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
-  return closestAtom;
+}
-/**
- * Gets all atoms connected to the provided atom by triangulation lines.
+ *
- * @param atom of which get all connected atoms
- * @param atom to be checked whether it is an inner atom
+ *
- * @return list of the two atoms linked to the provided one and closest to the atom to be checked,
- */
-list<atom*> * Tesselation::getClosestConnectedAtoms(atom* Atom, atom* AtomToCheck) {
-  list<atom*> connectedAtoms;
-  map<double, atom*> anglesOfAtoms;
-  map<double, atom*>::iterator runner;
-  LineMap::iterator findLines = LinesOnBoundary.begin();
-  list<atom*>::iterator listRunner;
-  Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
-  atom* current;
-  bool takeAtom = false;
-  planeNorm.CopyVector(&Atom->x);
-  planeNorm.SubtractVector(&AtomToCheck->x);
-  planeNorm.Normalize();
-  while (findLines != LinesOnBoundary.end()) {
-    takeAtom = false;
-    if (findLines->second->endpoints[0]->Nr == Atom->nr) {
-      takeAtom = true;
-      current = findLines->second->endpoints[1]->node;
-    } else if (findLines->second->endpoints[1]->Nr == Atom->nr) {
-      takeAtom = true;
-      current = findLines->second->endpoints[0]->node;
+    }
-    if (takeAtom) {
-      connectedAtoms.push_back(current);
-      currentPoint.CopyVector(&current->x);
-      currentPoint.ProjectOntoPlane(&planeNorm);
-      center.AddVector(&currentPoint);
+    }
-    findLines++;
+  }
-  cout << "Summed vectors " << center << "; number of atoms " << connectedAtoms.size()
-    << "; scale factor " << 1.0/connectedAtoms.size();
-  center.Scale(1.0/connectedAtoms.size());
-  listRunner = connectedAtoms.begin();
-  cout << " calculated center " << center <<  endl;
-  // construct one orthogonal vector
-  helper.CopyVector(&AtomToCheck->x);
-  helper.ProjectOntoPlane(&planeNorm);
-  OrthogonalVector.MakeNormalVector(&center, &helper, &(*listRunner)->x);
-  while (listRunner != connectedAtoms.end()) {
-    double angle = getAngle((*listRunner)->x, AtomToCheck->x, center, OrthogonalVector);
-    cout << "Calculated angle " << angle << " for atom " << **listRunner << endl;
-    anglesOfAtoms.insert(pair<double, atom*>(angle, (*listRunner)));
-    listRunner++;
+  }
-  list<atom*> *result = new list<atom*>;
-  runner = anglesOfAtoms.begin();
-  cout << "First value is " << *runner->second << endl;
-  result->push_back(runner->second);
-  runner = anglesOfAtoms.end();
-  runner--;
-  cout << "Second value is " << *runner->second << endl;
-  result->push_back(runner->second);
-  cout << "List of closest atoms has " << result->size() << " elements, which are "
-    << *(result->front()) << " and " << *(result->back()) << endl;
-  return result;
+}
-/**
- * Finds triangles belonging to the three provided atoms.
+ *
- * @param atom list, is expected to contain three atoms
+ *
- * @return triangles which belong to the provided atoms, will be empty if there are none,
- *         will usually be one, in case of degeneration, there will be two
- */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(atom* Points[3]) {
-  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
-  LineMap::iterator FindLine;
-  PointMap::iterator FindPoint;
-  TriangleMap::iterator FindTriangle;
-  class BoundaryPointSet *TrianglePoints[3];
-  for (int i = 0; i < 3; i++) {
-    FindPoint = PointsOnBoundary.find(Points[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      TrianglePoints[i] = FindPoint->second;
-    } else {
-      TrianglePoints[i] = NULL;
+    }
+  }
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (TrianglePoints[i] != NULL) {
-      for (int j = i; j < 3; j++) {
-        if (TrianglePoints[j] != NULL) {
-          FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr);
-          if (FindLine != TrianglePoints[i]->lines.end()) {
-            for (; FindLine->first == TrianglePoints[j]->node->nr; FindLine++) {
-              FindTriangle = FindLine->second->triangles.begin();
-              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if ((
-                  (FindTriangle->second->endpoints[0] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[1] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[2] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
+                  )
-                ) {
-                  result->push_back(FindTriangle->second);
+                }
+              }
+            }
-            // Is it sufficient to consider one of the triangle lines for this.
-            return result;
+          }
+        }
+      }
+    }
+  }
-  return result;
+}
-/**
- * Gets the angle between a point and a reference relative to the provided center.
+ *
- * @param point to calculate the angle for
- * @param reference to which to calculate the angle
- * @param center for which to calculate the angle between the vectors
- * @param OrthogonalVector helps discern between [0,pi] and [pi,2pi] in acos()
+ *
- * @return angle between point and reference
- */
-double getAngle(Vector point, Vector reference, Vector center, Vector OrthogonalVector) {
-  Vector ReferenceVector, helper;
-  cout << Verbose(2) << reference << " is our reference " << endl;
-  cout << Verbose(2) << center << " is our center " << endl;
-  // create baseline vector
-  ReferenceVector.CopyVector(&reference);
-  ReferenceVector.SubtractVector(&center);
-  if (ReferenceVector.IsNull())
-    return M_PI;
-  // calculate both angles and correct with in-plane vector
-  helper.CopyVector(&point);
-  helper.SubtractVector(&center);
-  if (helper.IsNull())
-    return M_PI;
-  double phi = ReferenceVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    phi = 2.*M_PI - phi;
+  }
-  cout << Verbose(2) << point << " has angle " << phi << endl;
-  return phi;
+}
-/**
- * Checks whether the provided point is within the tesselation structure.
+ *
- * This is a wrapper function for IsInnerAtom, so it can be used with a simple
- * vector, too.
+ *
- * @param point of which to check the position
- * @param tesselation structure
+ *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool IsInnerPoint(Vector Point, class Tesselation *Tess, LinkedCell* LC) {
-  atom *temp_atom = new atom;
-  bool status = false;
-  temp_atom->x.CopyVector(&Point);
-  status = IsInnerAtom(temp_atom, Tess, LC);
-  delete(temp_atom);
-  return status;
+}
 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
  * \param *out output stream for debugging
 …
   bool freeTess = false;
   bool freeLC = false;
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  int TriangleFilesWritten = 0;
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
   if (Tess == NULL) {
 …
+  }
   Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
+  Tess->Find_starting_triangle(out, RADIUS, LCList);
   baseline = Tess->LinesOnBoundary.begin();
   while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
     if (baseline->second->TrianglesCount == 1) {
       failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
+      failflag = Tess->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
 …
       flag = false;
+    }
+  }
+    // write temporary envelope
+    if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+      class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
+      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
+      if (DoTecplotOutput) {
+        string NameofTempFile(filename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(TecplotSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+      if (DoRaster3DOutput) {
+        string NameofTempFile(filename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(Raster3DSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_raster3d_file(out, tempstream, Tess, mol);
+//        // include the current position of the virtual sphere in the temporary raster3d file
+//        // make the circumsphere's center absolute again
+//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
+//        helper.AddVector(BaseRay->endpoints[1]->node->node);
+//        helper.Scale(0.5);
+//        (*it)->OptCenter.AddVector(&helper);
+//        Vector *center = mol->DetermineCenterOfAll(out);
+//        (*it)->OptCenter.SubtractVector(center);
+//        delete(center);
+//        // and add to file plus translucency object
+//        *tempstream << "# current virtual sphere\n";
+//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
+//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
+//          << "\t" << RADIUS << "\t1 0 0\n";
+//        *tempstream << "9\n  terminating special property\n";
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+      if (DoTecplotOutput || DoRaster3DOutput)
+        TriangleFilesWritten++;
+    }
+  }
+  // write final envelope
   if (filename != 0) {
     *out << Verbose(1) << "Writing final tecplot file\n";
 …
   cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(x, Tess, LCList)
     << "for vector " << x << "." << endl;
   atom* a = Tess->PointsOnBoundary.begin()->second->node;
+  TesselPoint* a = Tess->PointsOnBoundary.begin()->second->node;
   cout << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
   cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
+  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerPoint(a, Tess, LCList)
     << "for atom " << a << " (on boundary)." << endl;
   LinkedAtoms *List = NULL;
+  LinkedNodes *List = NULL;
   for (int i=0;i<NDIM;i++) { // each axis
     LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
 …
         //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
         if (List != NULL) {
           for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
+          for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
             if (Tess->PointsOnBoundary.find((*Runner)->nr) == Tess->PointsOnBoundary.end()) {
               a = *Runner;
 …
+      }
+  }
   cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
+  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(a, Tess, LCList)
     << "for atom " << a << " (inside)." << endl;

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset ab1932 for src/boundary.cpp

Legend:

src/boundary.cpp

Download in other formats: