Diff [87e2e39b6683b88cf56f1229adaf4e1bc2358143:c111dbf810e0e3f7f92d4afa3b50cd427f3d515d] for / – MoleCuilder

configure.ac

-              r87e2e39
+              rc111db
 # Checks for libraries.
 AC_CHECK_LIB(m, sqrt, ,AC_MSG_ERROR([compatible libc math library not found]))
+# Boost libraries
+AX_BOOST_BASE([1.33.1])
+AX_BOOST_PROGRAM_OPTIONS
+#AX_BOOST_FOREACH
+#AX_BOOST_FILESYSTEM
+#AX_BOOST_THREAD
+#AX_BOOST_PROGRAM_OPTIONS
+#AX_BOOST_SERIALIZATION
 # Checks for header files.

src/Makefile.am

-              r87e2e39
+              rc111db
 SOURCE = atom.cpp bond.cpp boundary.cpp config.cpp          element.cpp ellipsoid.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp memoryusageobserver.cpp molecules.cpp moleculelist.cpp parser.cpp periodentafel.cpp                tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
 HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp memoryusageobserver.hpp molecules.hpp                  parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
+SOURCE = atom.cpp bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp graph.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
+HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp graph.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
+BOOST_LIB = $(BOOST_LDFLAGS) $(BOOST_MPL_LIB)
 INCLUDES = -I$(top_srcdir)/src/unittests
 …
 libmolecuilder_a_SOURCES = ${SOURCE} ${HEADER}
 molecuilder_DATA = elements.db valence.db orbitals.db Hbonddistance.db Hbondangle.db
+molecuilder_LDFLAGS = $(BOOST_LIB)
 molecuilder_SOURCES = builder.cpp
 molecuilder_LDADD = libmolecuilder.a

src/atom.cpp

-              r87e2e39
+              rc111db
 #include "atom.hpp"
+#include "bond.hpp"
+#include "element.hpp"
 #include "memoryallocator.hpp"
+#include "vector.hpp"
 /************************************* Functions for class atom *************************************/
 …
 atom::~atom()
+{
+  Free(&ComponentNr);
+  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
+  Free<char>(&Name, "atom::~atom: *Name");
+  Trajectory.R.clear();
+  Trajectory.U.clear();
+  Trajectory.F.clear();
 };
 …
 };
+/** Sets father to itself or its father in case of copying a molecule.
+ */
+void atom::CorrectFather()
+{
+  if (father->father == father)   // same atom in copy's father points to itself
+    father = this;  // set father to itself (copy of a whole molecule)
+  else
+   father = father->father;  // set father to original's father
+};
+/** Check whether father is equal to given atom.
+ * \param *ptr atom to compare father to
+ * \param **res return value (only set if atom::father is equal to \a *ptr)
+ */
+void atom::EqualsFather ( atom *ptr, atom **res )
+{
+  if ( ptr == father )
+    *res = this;
+};
+/** Checks whether atom is within the given box.
+ * \param offset offset to box origin
+ * \param *parallelepiped box matrix
+ * \return true - is inside, false - is not
+ */
+bool atom::IsInParallelepiped(Vector offset, double *parallelepiped)
+{
+  return (node->IsInParallelepiped(offset, parallelepiped));
+};
 /** Output of a single atom.
  * \param ElementNo cardinal number of the element
 …
  * \param *out stream to output to
  * \param *comment commentary after '#' sign
+ */
+bool atom::Output(int ElementNo, int AtomNo, ofstream *out, const char *comment) const
+  * \return true - \a *out present, false - \a *out is NULL
+ */
+bool atom::Output(ofstream *out, int ElementNo, int AtomNo, const char *comment) const
+{
   if (out != NULL) {
 …
     return false;
 };
+bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment)
+{
+  AtomNo[type->Z]++;  // increment number
+  if (out != NULL) {
+    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
+    *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
+    *out << "\t" << FixedIon;
+    if (v.Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
+    if (comment != NULL)
+      *out << " # " << comment << endl;
+    else
+      *out << " # molecule nr " << nr << endl;
+    return true;
+  } else
+    return false;
+};
 /** Output of a single atom as one lin in xyz file.
  * \param *out stream to output to
+  * \return true - \a *out present, false - \a *out is NULL
  */
 bool atom::OutputXYZLine(ofstream *out) const
 …
   if (out != NULL) {
     *out << type->symbol << "\t" << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2] << "\t" << endl;
+    return true;
+  } else
+    return false;
+};
+/** Output of a single atom as one lin in xyz file.
+ * \param *out stream to output to
+ * \param *ElementNo array with ion type number in the config file this atom's element shall have
+ * \param *AtomNo array with atom number in the config file this atom shall have, is increase by one automatically
+ * \param step Trajectory time step to output
+  * \return true - \a *out present, false - \a *out is NULL
+ */
+bool atom::OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const
+{
+  AtomNo[type->Z]++;
+  if (out != NULL) {
+    *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
+    *out << Trajectory.R.at(step).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2];
+    *out << "\t" << FixedIon;
+    if (Trajectory.U.at(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step).x[0] << "\t" << Trajectory.U.at(step).x[1] << "\t" << Trajectory.U.at(step).x[2] << "\t";
+    if (Trajectory.F.at(step).Norm() > MYEPSILON)
+      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step).x[0] << "\t" << Trajectory.F.at(step).x[1] << "\t" << Trajectory.F.at(step).x[2] << "\t";
+    *out << "\t# Number in molecule " << nr << endl;
+    return true;
+  } else
+    return false;
+};
+/** Output of a single atom as one lin in xyz file.
+ * \param *out stream to output to
+ * \param step Trajectory time step to output
+ * \return true - \a *out present, false - \a *out is NULL
+ */
+bool atom::OutputTrajectoryXYZ(ofstream *out, int step) const
+{
+  if (out != NULL) {
+    *out << type->symbol << "\t";
+    *out << Trajectory.R.at(step).x[0] << "\t";
+    *out << Trajectory.R.at(step).x[1] << "\t";
+    *out << Trajectory.R.at(step).x[2] << endl;
+    return true;
+  } else
+    return false;
+};
+/** Prints all bonds of this atom from given global lists.
+ * \param *out stream to output to
+ * \param *NumberOfBondsPerAtom array with number of bonds per atomic index
+ * \param ***ListOfBondsPerAtom array per atomic index of array with pointer to bond
+ * \return true - \a *out present, false - \a *out is NULL
+ */
+bool atom::OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const
+{
+  if (out != NULL) {
+#ifdef ADDHYDROGEN
+    if (type->Z != 1) {   // regard only non-hydrogen
+#endif
+      *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << NumberOfBondsPerAtom[nr] << " bonds: ";
+      int TotalDegree = 0;
+      for (int j=0;j<NumberOfBondsPerAtom[nr];j++) {
+        *out << *ListOfBondsPerAtom[nr][j] << "\t";
+        TotalDegree += ListOfBondsPerAtom[nr][j]->BondDegree;
+      }
+      *out << " -- TotalDegree: " << TotalDegree << endl;
+#ifdef ADDHYDROGEN
+    }
+#endif
     return true;
   } else
 …
 };
+/** Returns squared distance to a given vector.
+ * \param origin vector to calculate distance to
+ * \return distance squared
+ */
+double atom::DistanceSquaredToVector(Vector &origin)
+{
+  return origin.DistanceSquared(&x);
+};
+/** Adds kinetic energy of this atom to given temperature value.
+ * \param *temperature add on this value
+ * \param step given step of trajectory to add
+ */
+void atom::AddKineticToTemperature(double *temperature, int step) const
+{
+  for (int i=NDIM;i--;)
+    *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
+};
+/** Returns distance to a given vector.
+ * \param origin vector to calculate distance to
+ * \return distance
+ */
+double atom::DistanceToVector(Vector &origin)
+{
+  return origin.Distance(&x);
+};
 bool operator < (atom &a, atom &b)
+{

src/atom.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #endif
+#include <iostream>
+#include <vector>
 #include <iostream>
+#include "tesselation.hpp"
+#include "element.hpp"
+#include "tesselation.hpp"
+#include "vector.hpp"
+/****************************************** forward declarations *****************************/
+class bond;
+class element;
+class Vector;
+/********************************************** declarations *******************************/
 /** Single atom.
 …
 class atom : public TesselPoint {
   public:
+    Vector x;       //!< coordinate array of atom, giving position within cell
+    Vector v;       //!< velocity array of atom
+    struct
+    {
+      vector<Vector> R;  //!< position vector
+      vector<Vector> U;  //!< velocity vector
+      vector<Vector> F;  //!< last force vector
+    } Trajectory;
+    Vector x;       //!< coordinate vector of atom, giving last position within cell
+    Vector v;       //!< velocity vector of atom, giving last velocity within cell
+    Vector F;       //!< Force vector of atom, giving last force within cell
     element *type;  //!< pointing to element
     atom *previous; //!< previous atom in molecule list
 …
   virtual ~atom();
+  bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
+  bool Output(ofstream *out, int ElementNo, int AtomNo, const char *comment = NULL) const;
+  bool Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment = NULL);
   bool OutputXYZLine(ofstream *out) const;
+  bool OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const;
+  bool OutputTrajectoryXYZ(ofstream *out, int step) const;
+  bool OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const;
+  void EqualsFather ( atom *ptr, atom **res );
+  void CorrectFather();
   atom *GetTrueFather();
   bool Compare(const atom &ptr);
+  double DistanceToVector(Vector &origin);
+  double DistanceSquaredToVector(Vector &origin);
+  void AddKineticToTemperature(double *temperature, int step) const;
+  bool IsInParallelepiped(Vector offset, double *parallelepiped);
   ostream & operator << (ostream &ost);

src/bond.cpp

-              r87e2e39
+              rc111db
  */
+#include "atom.hpp"
 #include "bond.hpp"
+#include "element.hpp"
+#include "lists.hpp"
 /***************************************** Functions for class bond ********************************/

src/bond.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #endif
+#include "atom.hpp"
+/****************************************** forward declarations *****************************/
+class atom;
+/********************************************** declarations *******************************/
 /** Bonds between atoms.

src/boundary.cpp

-              r87e2e39
+              rc111db
 /** \file boundary.hpp
+/** \file boundary.cpp
+ *
  * Implementations and super-function for envelopes
  */
+#include "atom.hpp"
+#include "bond.hpp"
 #include "boundary.hpp"
+#include "config.hpp"
+#include "element.hpp"
+#include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "memoryallocator.hpp"
+#include "molecule.hpp"
+#include "tesselation.hpp"
+#include "tesselationhelpers.hpp"
 #include<gsl/gsl_poly.h>
 …
  * \return NDIM array of the diameters
  */
+double *
+GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
+    bool IsAngstroem)
+double *GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
+{
   // get points on boundary of NULL was given as parameter
 …
+;
-/** Creates the objects in a VRML file.
- * \param *out output stream for debugging
- * \param *vrmlfile output stream for tecplot data
- * \param *Tess Tesselation structure with constructed triangles
- * \param *mol molecule structure with atom positions
- */
-void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
+{
-  atom *Walker = mol->start;
-  bond *Binder = mol->first;
-  int i;
-  Vector *center = mol->DetermineCenterOfAll(out);
-  if (vrmlfile != NULL) {
-    //cout << Verbose(1) << "Writing Raster3D file ... ";
-    *vrmlfile << "#VRML V2.0 utf8" << endl;
-    *vrmlfile << "#Created by molecuilder" << endl;
-    *vrmlfile << "#All atoms as spheres" << endl;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
-      for (i=0;i<NDIM;i++)
-        *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
-      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
-    *vrmlfile << "# All bonds as vertices" << endl;
-    while (Binder->next != mol->last) {
-      Binder = Binder->next;
-      *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
-      for (i=0;i<NDIM;i++)
-        *vrmlfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
-      *vrmlfile << "\t0.03\t";
-      for (i=0;i<NDIM;i++)
-        *vrmlfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
-      *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
-    *vrmlfile << "# All tesselation triangles" << endl;
-    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
-      *vrmlfile << "1" << endl << "  "; // 1 is triangle type
-      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->node->x[j]-center->x[j] << " ";
-        *vrmlfile << "\t";
+      }
-      *vrmlfile << "1. 0. 0." << endl;  // red as colour
-      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+    }
-  } else {
-    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
+  }
-  delete(center);
-};
-/** Creates the objects in a raster3d file (renderable with a header.r3d).
- * \param *out output stream for debugging
- * \param *rasterfile output stream for tecplot data
- * \param *Tess Tesselation structure with constructed triangles
- * \param *mol molecule structure with atom positions
- */
-void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
+{
-  atom *Walker = mol->start;
-  bond *Binder = mol->first;
-  int i;
-  Vector *center = mol->DetermineCenterOfAll(out);
-  if (rasterfile != NULL) {
-    //cout << Verbose(1) << "Writing Raster3D file ... ";
-    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
-    *rasterfile << "@header.r3d" << endl;
-    *rasterfile << "# All atoms as spheres" << endl;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
-      for (i=0;i<NDIM;i++)
-        *rasterfile << Walker->x.x[i]-center->x[i] << " ";
-      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
-    *rasterfile << "# All bonds as vertices" << endl;
-    while (Binder->next != mol->last) {
-      Binder = Binder->next;
-      *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
-      for (i=0;i<NDIM;i++)
-        *rasterfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
-      *rasterfile << "\t0.03\t";
-      for (i=0;i<NDIM;i++)
-        *rasterfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
-      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
-    *rasterfile << "# All tesselation triangles" << endl;
-    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
-    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
-      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
-      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->node->x[j]-center->x[j] << " ";
-        *rasterfile << "\t";
+      }
-      *rasterfile << "1. 0. 0." << endl;  // red as colour
-      //*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
+    }
-    *rasterfile << "9\n#  terminating special property\n";
-  } else {
-    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
-  delete(center);
-};
-/** This function creates the tecplot file, displaying the tesselation of the hull.
- * \param *out output stream for debugging
- * \param *tecplot output stream for tecplot data
- * \param N arbitrary number to differentiate various zones in the tecplot format
- */
-void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
+{
-  if ((tecplot != NULL) && (TesselStruct != NULL)) {
-    // write header
-    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
-    int *LookupList = new int[mol->AtomCount];
-    for (int i = 0; i < mol->AtomCount; i++)
-      LookupList[i] = -1;
-    // print atom coordinates
-    *out << Verbose(2) << "The following triangles were created:";
-    int Counter = 1;
-    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->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
+    }
-    *tecplot << endl;
-    // print connectivity
-    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
-      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
-      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+    }
-    delete[] (LookupList);
-    *out << endl;
+  }
+}
 /** Determines the boundary points of a cluster.
 …
         // flip the line
         if (!mol->TesselStruct->PickFarthestofTwoBaselines(out, line))
+        if (mol->TesselStruct->PickFarthestofTwoBaselines(out, line) == 0.)
           *out << Verbose(1) << "ERROR: Correction of concave baselines failed!" << endl;
+        else
+        else {
+          mol->TesselStruct->FlipBaseline(out, line);
           *out << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
+        }
+      }
+    }
 …
   cout << Verbose(1) << "End of FindConvexBorder" << endl;
+};
+/** For testing removes one boundary point after another to check for leaks.
+ * \param *out output stream for debugging
+ * \param *TesselStruct Tesselation containing envelope with boundary points
+ * \param *mol molecule
+ * \param *filename name of file
+ * \return true - all removed, false - something went wrong
+ */
+bool RemoveAllBoundaryPoints(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename)
+{
+  int i=0;
+  char number[MAXSTRINGSIZE];
+  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
+    *out << Verbose(2) << "ERROR: TesselStruct is empty." << endl;
+    return false;
+  }
+  PointMap::iterator PointRunner;
+  while (!TesselStruct->PointsOnBoundary.empty()) {
+    *out << Verbose(2) << "Remaining points are: ";
+    for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
+      *out << *(PointSprinter->second) << "\t";
+      *out << endl;
+    PointRunner = TesselStruct->PointsOnBoundary.begin();
+    // remove point
+    TesselStruct->RemovePointFromTesselatedSurface(out, PointRunner->second);
+    // store envelope
+    sprintf(number, "-%04d", i++);
+    StoreTrianglesinFile(out, mol, filename, number);
+  }
+  return true;
 };
 …
   class BoundaryLineSet *line = NULL;
   bool Concavity;
+  char dummy[MAXSTRINGSIZE];
   PointMap::iterator PointRunner, PointAdvance;
   LineMap::iterator LineRunner, LineAdvance;
 …
+  }
-  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
-  StoreTrianglesinFile(out, mol, filename, "-first");
   // First step: RemovePointFromTesselatedSurface
+  int run = 0;
+  double tmp;
   do {
     Concavity = false;
+    sprintf(dummy, "-first-%d", run);
+    //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+    StoreTrianglesinFile(out, mol, filename, dummy);
     PointRunner = TesselStruct->PointsOnBoundary.begin();
     PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
 …
         line = LineRunner->second;
         *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+      }
+      if (!line->CheckConvexityCriterion(out)) {
+        *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
+        // remove the point
+        Concavity = true;
+        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+        if (!line->CheckConvexityCriterion(out)) {
+          // remove the point if needed
+          *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
+          volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
+          sprintf(dummy, "-first-%d", ++run);
+          StoreTrianglesinFile(out, mol, filename, dummy);
+          Concavity = true;
+          break;
+        }
+      }
       PointRunner = PointAdvance;
+    }
+    sprintf(dummy, "-second-%d", run);
     //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
     //StoreTrianglesinFile(out, mol, filename, "-second");
+    StoreTrianglesinFile(out, mol, filename, dummy);
     // second step: PickFarthestofTwoBaselines
 …
       // take highest of both lines
       if (TesselStruct->IsConvexRectangle(out, line) == NULL) {
+        TesselStruct->PickFarthestofTwoBaselines(out, line);
+        Concavity = true;
+        tmp = TesselStruct->PickFarthestofTwoBaselines(out, line);
+        volume += tmp;
+        if (tmp != 0) {
+          mol->TesselStruct->FlipBaseline(out, line);
+          Concavity = true;
+        }
+      }
       LineRunner = LineAdvance;
+    }
+    run++;
   } while (Concavity);
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  //StoreTrianglesinFile(out, mol, filename, "-third");
+  // third step: IsConvexRectangle
+//  LineRunner = TesselStruct->LinesOnBoundary.begin();
+//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+//    LineAdvance++;
+//    line = LineRunner->second;
+//    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
+//    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+//      //*out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
+//    if (!line->CheckConvexityCriterion(out)) {
+//      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
+//
+//      // take highest of both lines
+//      point = TesselStruct->IsConvexRectangle(out, line);
+//      if (point != NULL)
+//        volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
+//    }
+//    LineRunner = LineAdvance;
+//  }
   CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-third");
+  // third step: IsConvexRectangle
+  LineRunner = TesselStruct->LinesOnBoundary.begin();
+  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+    LineAdvance++;
+    line = LineRunner->second;
+    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
+    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+      //*out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
+    if (!line->CheckConvexityCriterion(out)) {
+      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
+      // take highest of both lines
+      point = TesselStruct->IsConvexRectangle(out, line);
+      if (point != NULL)
+        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+    }
+    LineRunner = LineAdvance;
+  }
+  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-fourth");
+  StoreTrianglesinFile(out, mol, filename, "");
   // end
+  *out << Verbose(1) << "Volume is " << volume << "." << endl;
   *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
   return volume;
 };
-/** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
- * Sets BoundaryPointSet::value equal to the number of connected lines that are not convex.
- * \param *out output stream for debugging
- * \param *TesselStruct pointer to Tesselation structure
- */
-void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct)
+{
-  class BoundaryPointSet *point = NULL;
-  class BoundaryLineSet *line = NULL;
-  // calculate remaining concavity
-  for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
-    point = PointRunner->second;
-    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
-    point->value = 0;
-    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
-      line = LineRunner->second;
-      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
-      if (!line->CheckConvexityCriterion(out))
-        point->value += 1;
+    }
+  }
-};
-/** Stores triangles to file.
- * \param *out output stream for debugging
- * \param *mol molecule with atoms and bonds
- * \param *filename prefix of filename
- * \param *extraSuffix intermediate suffix
- */
-void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix)
+{
-  // 4. Store triangles in tecplot file
-  if (filename != NULL) {
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append(extraSuffix);
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
-      tecplot->close();
-      delete(tecplot);
+    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append(extraSuffix);
-      OutputName.append(Raster3DSuffix);
-      ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
-      rasterplot->close();
-      delete(rasterplot);
+    }
+  }
-};
 /** Determines the volume of a cluster.
 …
   return volume;
+}
+;
+};
+/** Stores triangles to file.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms and bonds
+ * \param *filename prefix of filename
+ * \param *extraSuffix intermediate suffix
+ */
+void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix)
+{
+  // 4. Store triangles in tecplot file
+  if (filename != NULL) {
+    if (DoTecplotOutput) {
+      string OutputName(filename);
+      OutputName.append(extraSuffix);
+      OutputName.append(TecplotSuffix);
+      ofstream *tecplot = new ofstream(OutputName.c_str());
+      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
+      tecplot->close();
+      delete(tecplot);
+    }
+    if (DoRaster3DOutput) {
+      string OutputName(filename);
+      OutputName.append(extraSuffix);
+      OutputName.append(Raster3DSuffix);
+      ofstream *rasterplot = new ofstream(OutputName.c_str());
+      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
+      rasterplot->close();
+      delete(rasterplot);
+    }
+  }
+};
 /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
 …
   int N[NDIM];
   int n[NDIM];
   double *M =  filler->ReturnFullMatrixforSymmetric(filler->cell_size);
+  double *M =  ReturnFullMatrixforSymmetric(filler->cell_size);
   double Rotations[NDIM*NDIM];
   Vector AtomTranslations;
 …
     if ((*ListRunner)->TesselStruct == NULL) {
       *out << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
       FindNonConvexBorder((ofstream *)&cout, (*ListRunner), LCList[i], NULL, 5.);
+      FindNonConvexBorder((ofstream *)&cout, (*ListRunner), LCList[i], 5., NULL);
+    }
     i++;
 …
  * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
  * \param *LCList atoms in LinkedCell list
+ * \param RADIUS radius of the virtual sphere
  * \param *filename filename prefix for output of vertex data
+ * \para RADIUS radius of the virtual sphere
+ */
+void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename, const double RADIUS)
+ */
+void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const double RADIUS, const char *filename = NULL)
+{
-  int N = 0;
   bool freeLC = false;
-  ofstream *tempstream = NULL;
-  char NumberName[255];
-  int TriangleFilesWritten = 0;
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
 …
   LineMap::iterator testline;
   *out << Verbose(0) << "Begin of FindNonConvexBorder\n";
+  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
+  bool failflag = false;
+  bool OneLoopWithoutSuccessFlag = false;  // marks whether we went once through all baselines without finding any without two triangles
+  bool TesselationFailFlag = false;
+  // initialise Linked Cell
   if (LCList == NULL) {
     LCList = new LinkedCell(mol, 2.*RADIUS);
 …
+  }
+  // 1. get starting triangle
   mol->TesselStruct->FindStartingTriangle(out, RADIUS, LCList);
+  // 2. expand from there
   baseline = mol->TesselStruct->LinesOnBoundary.begin();
+  // the outward most line is dangerous, as we may end up with wrapping up the starting triangle, hence
+  // terminating the algorithm too early.
+  if (baseline != mol->TesselStruct->LinesOnBoundary.end()) // skip first line as it its the outwardmost!
+        baseline++;
+  while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (flag)) {
+  baseline++; // skip first line
+  while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (OneLoopWithoutSuccessFlag)) {
     if (baseline->second->triangles.size() == 1) {
+      failflag = mol->TesselStruct->FindNextSuitableTriangle(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)
+      // 3. find next triangle
+      TesselationFailFlag = mol->TesselStruct->FindNextSuitableTriangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one.
+      OneLoopWithoutSuccessFlag = OneLoopWithoutSuccessFlag || TesselationFailFlag;
+      if (!TesselationFailFlag)
         cerr << "WARNING: FindNextSuitableTriangle failed." << endl;
       // write temporary envelope
+      if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % SingleStepWidth == 0))) { // if we have a new triangle and want to output each new triangle configuration
+        TriangleMap::iterator runner = mol->TesselStruct->TrianglesOnBoundary.end();
+        runner--;
+        class BoundaryTriangleSet *triangle = runner->second;
+        if (triangle != NULL) {
+          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);
+            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+            WriteTecplotFile(out, tempstream, mol->TesselStruct, 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);
+            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+            WriteRaster3dFile(out, tempstream, mol->TesselStruct, 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++;
+      if (filename != NULL) {
+        if ((DoSingleStepOutput && ((mol->TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
+          mol->TesselStruct->Output(out, filename, mol);
+        }
+      }
+      baseline = mol->TesselStruct->LinesOnBoundary.end();
+      *out << Verbose(2) << "Baseline set to end." << endl;
     } else {
       //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->triangles.size() << " triangles adjacent" << endl;
 …
+    }
+    N++;
+    if ((baseline == mol->TesselStruct->LinesOnBoundary.end()) && (OneLoopWithoutSuccessFlag)) {
+      baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+      OneLoopWithoutSuccessFlag = false;
+    }
     baseline++;
+    if ((baseline == mol->TesselStruct->LinesOnBoundary.end()) && (flag)) {
+      baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+      flag = false;
+    }
+  }
+  }
+  // check envelope for consistency
+  CheckListOfBaselines(out, mol->TesselStruct);
+  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
+  //mol->TesselStruct->InsertStraddlingPoints(out, mol, LCList);
+//  mol->GoToFirst();
+//  class TesselPoint *Runner = NULL;
+//  while (!mol->IsEnd()) {
+//    Runner = mol->GetPoint();
+//    *out << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
+//    if (!mol->TesselStruct->IsInnerPoint(out, Runner, LCList)) {
+//      *out << Verbose(2) << Runner->Name << " is outside of envelope, adding via degenerated triangles." << endl;
+//      mol->TesselStruct->AddBoundaryPointByDegeneratedTriangle(out, Runner, LCList);
+//    } else {
+//      *out << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
+//    }
+//    mol->GoToNext();
+//  }
   // Purges surplus triangles.
   mol->TesselStruct->RemoveDegeneratedTriangles();
+  // check envelope for consistency
+  CheckListOfBaselines(out, mol->TesselStruct);
   // write final envelope
+  if (filename != 0) {
+    *out << Verbose(1) << "Writing final tecplot file\n";
+    if (DoTecplotOutput) {
+      string OutputName(filename);
+      OutputName.append(TecplotSuffix);
+      ofstream *tecplot = new ofstream(OutputName.c_str());
+      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, -1);
+      tecplot->close();
+      delete(tecplot);
+    }
+    if (DoRaster3DOutput) {
+      string OutputName(filename);
+      OutputName.append(Raster3DSuffix);
+      ofstream *raster = new ofstream(OutputName.c_str());
+      WriteRaster3dFile(out, raster, mol->TesselStruct, mol);
+      raster->close();
+      delete(raster);
+    }
+  } else {
+    cerr << "ERROR: Could definitively not find all necessary triangles!" << endl;
+  }
+  cout << Verbose(2) << "Check: List of Baselines with not two connected triangles:" << endl;
+  int counter = 0;
+  for (testline = mol->TesselStruct->LinesOnBoundary.begin(); testline != mol->TesselStruct->LinesOnBoundary.end(); testline++) {
+    if (testline->second->triangles.size() != 2) {
+      cout << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl;
+      counter++;
+    }
+  }
+  if (counter == 0)
+    *out << Verbose(2) << "None." << endl;
+//  // Tests the IsInnerAtom() function.
+//  Vector x (0, 0, 0);
+//  *out << Verbose(0) << "Point to check: " << x << endl;
+//  *out << Verbose(0) << "Check: IsInnerPoint() returns " << mol->TesselStruct->IsInnerPoint(out, x, LCList)
+//    << "for vector " << x << "." << endl;
+//  TesselPoint* a = mol->TesselStruct->PointsOnBoundary.begin()->second->node;
+//  *out << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
+//  *out << Verbose(0) << "Check: IsInnerAtom() returns " << mol->TesselStruct->IsInnerPoint(out, a, LCList)
+//    << "for atom " << a << " (on boundary)." << endl;
+//  LinkedNodes *List = NULL;
+//  for (int i=0;i<NDIM;i++) { // each axis
+//    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
+//    for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
+//      for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
+//        List = LCList->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 (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
+//            if (mol->TesselStruct->PointsOnBoundary.find((*Runner)->nr) == mol->TesselStruct->PointsOnBoundary.end()) {
+//              a = *Runner;
+//              i=3;
+//              for (int j=0;j<NDIM;j++)
+//                LCList->n[j] = LCList->N[j];
+//              break;
+//            }
+//          }
+//        }
+//      }
+//  }
+//  *out << Verbose(0) << "Check: IsInnerPoint() returns " << mol->TesselStruct->IsInnerPoint(out, a, LCList)
+//    << "for atom " << a << " (inside)." << endl;
+  CalculateConcavityPerBoundaryPoint(out, mol->TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "");
   if (freeLC)
 …
   *out << Verbose(0) << "End of FindNonConvexBorder\n";
 };
 /** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.

src/boundary.hpp

-              r87e2e39
+              rc111db
 #ifndef BOUNDARY_HPP_
 #define BOUNDARY_HPP_
+using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 …
 #endif
+#include <fstream>
+#include <iostream>
 // STL headers
 #include <map>
+#include "config.hpp"
+#include "linkedcell.hpp"
+#include "molecules.hpp"
+#include "tesselation.hpp"
+#include "defs.hpp"
+/****************************************** forward declarations *****************************/
+class atom;
+class BoundaryPointSet;
+class BoundaryLineSet;
+class BoundaryTriangleSet;
+class config;
+class LinkedCell;
+class molecule;
+class MoleculeListClass;
+class Tesselation;
+class Vector;
+/********************************************** definitions *********************************/
 #define DEBUG 1
 #define DoSingleStepOutput 0
 #define SingleStepWidth 1
-#define DoTecplotOutput 1
-#define DoRaster3DOutput 1
-#define DoVRMLOutput 1
-#define TecplotSuffix ".dat"
-#define Raster3DSuffix ".r3d"
-#define VRMLSUffix ".wrl"
 #define DistancePair pair < double, atom* >
 …
 #define BoundariesTestPair pair< Boundaries::iterator, bool>
+/********************************************** declarations *******************************/
 double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration);
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
 …
 molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
 void FindConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename);
 void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
+void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LC, const double RADIUS, const char *tempbasename);
 double ConvexizeNonconvexEnvelope(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename);
 void FindNextSuitablePoint(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
 Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);
-void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct);
 void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix);
+bool RemoveAllBoundaryPoints(ofstream *out, class Tesselation *TesselStruct, molecule *mol, char *filename);

src/builder.cpp

-              r87e2e39
+              rc111db
 using namespace std;
+#include "atom.hpp"
+#include "bond.hpp"
 #include "boundary.hpp"
+#include "config.hpp"
+#include "element.hpp"
 #include "ellipsoid.hpp"
 #include "helpers.hpp"
+#include "leastsquaremin.hpp"
+#include "linkedcell.hpp"
 #include "memoryusageobserverunittest.hpp"
+#include "molecules.hpp"
+#include "molecule.hpp"
+#include "periodentafel.hpp"
 /********************************************* Subsubmenu routine ************************************/
 …
       cin >> nr;
       count = 1;
       for( MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+      for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
         if (nr == (*ListRunner)->IndexNr) {
           mol = *ListRunner;
           molecules->ListOfMolecules.erase(ListRunner);
           delete(mol);
+          break;
+        }
       break;
 …
     case 'e':
+      cout << Verbose(0) << "Not implemented yet." << endl;
+      {
+        int src, dest;
+        molecule *srcmol = NULL, *destmol = NULL;
+        do {
+          cout << Verbose(0) << "Enter index of matrix molecule (the variable one): ";
+          cin >> src;
+          srcmol = molecules->ReturnIndex(src);
+        } while ((srcmol == NULL) && (src != -1));
+        do {
+          cout << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ";
+          cin >> dest;
+          destmol = molecules->ReturnIndex(dest);
+        } while ((destmol == NULL) && (dest != -1));
+        if ((src != -1) && (dest != -1))
+          molecules->EmbedMerge(destmol, srcmol);
+      }
       break;
 …
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
             cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-R\t\tRemove all atoms out of sphere around a given one." << endl;
+            cout << "\t-r <id>\t\tRemove an atom with given id." << endl;
+            cout << "\t-R <id> <radius>\t\tRemove all atoms out of sphere around a given one." << endl;
             cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
             cout << "\t-S <file> Store temperatures from the config file in <file>." << endl;
 …
           switch(argv[argptr-1][1]) {
             case 'p':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
 …
               break;
             case 'a':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {
                 ExitFlag = 255;
 …
               break;
             case 'D':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
+              {
                 cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
 …
               break;
             case 'E':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
                 ExitFlag = 255;
 …
               break;
             case 'F':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if (argptr+5 >=argc) {
                 ExitFlag = 255;
 …
               break;
             case 'A':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag =255;
 …
               break;
             case 'N':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
                 ExitFlag = 255;
 …
                 cout << Verbose(0) << "Evaluating non-convex envelope.";
                 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
                 start = clock();
+                start = clock();
                 LinkedCell LCList(mol, atof(argv[argptr])*2.);
                 FindNonConvexBorder((ofstream *)&cout, mol, &LCList, argv[argptr+1], atof(argv[argptr]));
+                FindNonConvexBorder((ofstream *)&cout, mol, &LCList, atof(argv[argptr]), argv[argptr+1]);
                 //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
                 end = clock();
                 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+                end = clock();
+                cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
                 argptr+=2;
+              }
               break;
             case 'S':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
 …
               break;
             case 'L':
+              ExitFlag = 1;
+              SaveFlag = true;
+              cout << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
+              if (!mol->LinearInterpolationBetweenConfiguration((ofstream *)&cout, atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration))
+                cout << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
+              else
+                cout << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
+              argptr+=3;
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl;
+              } else {
+                SaveFlag = true;
+                cout << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
+                if (atoi(argv[argptr+3]) == 1)
+                  cout << Verbose(1) << "Using Identity for the permutation map." << endl;
+                if (!mol->LinearInterpolationBetweenConfiguration((ofstream *)&cout, atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration, atoi(argv[argptr+3])) == 1 ? true : false)
+                  cout << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
+                else
+                  cout << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
+                argptr+=4;
+              }
               break;
             case 'P':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
 …
               break;
             case 'R':
               ExitFlag = 1;
               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])))  {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])))  {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl;
 …
               break;
             case 't':
               ExitFlag = 1;
               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
               } else {
                 ExitFlag = 1;
+                if (ExitFlag == 0) ExitFlag = 1;
                 SaveFlag = true;
                 cout << Verbose(1) << "Translating all ions to new origin." << endl;
+                cout << Verbose(1) << "Translating all ions by given vector." << endl;
                 for (int i=NDIM;i--;)
                   x.x[i] = atof(argv[argptr+i]);
 …
                 argptr+=3;
+              }
+              break;
             case 'T':
               ExitFlag = 1;
               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl;
               } else {
                 ExitFlag = 1;
+                if (ExitFlag == 0) ExitFlag = 1;
                 SaveFlag = true;
                 cout << Verbose(1) << "Translating all ions periodically to new origin." << endl;
+                cout << Verbose(1) << "Translating all ions periodically by given vector." << endl;
                 for (int i=NDIM;i--;)
                   x.x[i] = atof(argv[argptr+i]);
 …
               break;
             case 's':
               ExitFlag = 1;
               if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;
 …
               break;
             case 'b':
               ExitFlag = 1;
               if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
 …
               break;
             case 'B':
               ExitFlag = 1;
               if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
 …
               break;
             case 'c':
               ExitFlag = 1;
               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
 …
               break;
             case 'O':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               SaveFlag = true;
               cout << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl;
 …
               break;
             case 'r':
+              ExitFlag = 1;
+              SaveFlag = true;
+              cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])))  {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl;
+              } else {
+                SaveFlag = true;
+                cout << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl;
+                atom *first = mol->FindAtom(atoi(argv[argptr]));
+                mol->RemoveAtom(first);
+                argptr+=1;
+              }
               break;
             case 'f':
               ExitFlag = 1;
               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
 …
               break;
             case 'm':
               ExitFlag = 1;
+              if (ExitFlag == 0) ExitFlag = 1;
               j = atoi(argv[argptr++]);
               if ((j<0) || (j>1)) {
 …
               break;
             case 'o':
               ExitFlag = 1;
               if ((argptr >= argc) || (argv[argptr][0] == '-')){
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
+                cerr << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl;
               } else {
                 cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+                cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
+                cout << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;
+                cout << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl;
                 LinkedCell LCList(mol, 10.);
                 //FindConvexBorder((ofstream *)&cout, mol, &LCList, argv[argptr]);
                 FindNonConvexBorder((ofstream *)&cout, mol, &LCList, argv[argptr], 10.);
+                FindNonConvexBorder((ofstream *)&cout, mol, &LCList, 5., argv[argptr+1]);
+//                RemoveAllBoundaryPoints((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
                 double volumedifference = ConvexizeNonconvexEnvelope((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
                 double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, mol->TesselStruct, &configuration);
                 cout << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                 cout << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
                 argptr+=1;
+                argptr+=2;
+              }
               break;
             case 'U':
               ExitFlag = 1;
               if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
 …
+              }
             case 'u':
               ExitFlag = 1;
               if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
                 if (volume != -1)
                   ExitFlag = 255;
 …
               break;
             case 'd':
               ExitFlag = 1;
               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+              if (ExitFlag == 0) ExitFlag = 1;
+              if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
                 ExitFlag = 255;
                 cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
 …
       return 0;
       break;
+    case 2:  // just for -f option
+      delete(molecules); // also free's all molecules contained
+      delete(periode);
+      cout << Verbose(0) <<  "Maximum of allocated memory: "
+        << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
+      cout << Verbose(0) <<  "Remaining non-freed memory: "
+        << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
+      return 2;
+      break;
     default:
       break;

src/config.cpp

-              r87e2e39
+              rc111db
  */
+#include "atom.hpp"
 #include "config.hpp"
+#include "element.hpp"
 #include "memoryallocator.hpp"
+#include "molecule.hpp"
+#include "periodentafel.hpp"
 /******************************** Functions for class ConfigFileBuffer **********************/
 …
             // check size of vectors
             if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
+            if (neues->Trajectory.R.size() <= (unsigned int)(repetition)) {
               //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
               mol->Trajectories[neues].R.resize(repetition+10);
               mol->Trajectories[neues].U.resize(repetition+10);
               mol->Trajectories[neues].F.resize(repetition+10);
+              neues->Trajectory.R.resize(repetition+10);
+              neues->Trajectory.U.resize(repetition+10);
+              neues->Trajectory.F.resize(repetition+10);
+            }
             // put into trajectories list
             for (int d=0;d<NDIM;d++)
               mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
+              neues->Trajectory.R.at(repetition).x[d] = neues->x.x[d];
             // parse velocities if present
 …
               neues->v.x[2] = 0.;
             for (int d=0;d<NDIM;d++)
               mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
+              neues->Trajectory.U.at(repetition).x[d] = neues->v.x[d];
             // parse forces if present
 …
               value[2] = 0.;
             for (int d=0;d<NDIM;d++)
               mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
+              neues->Trajectory.F.at(repetition).x[d] = value[d];
   //            cout << "Parsed position of step " << (repetition) << ": (";
   //            for (int d=0;d<NDIM;d++)
   //              cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";          // next step
+  //              cout << neues->Trajectory.R.at(repetition).x[d] << " ";          // next step
   //            cout << ")\t(";
   //            for (int d=0;d<NDIM;d++)
   //              cout << mol->Trajectories[neues].U.at(repetition).x[d] << " ";          // next step
+  //              cout << neues->Trajectory.U.at(repetition).x[d] << " ";          // next step
   //            cout << ")\t(";
   //            for (int d=0;d<NDIM;d++)
   //              cout << mol->Trajectories[neues].F.at(repetition).x[d] << " ";          // next step
+  //              cout << neues->Trajectory.F.at(repetition).x[d] << " ";          // next step
   //            cout << ")" << endl;
+          }
 …
       repetition--;
       cout << "Found " << repetition << " trajectory steps." << endl;
+      mol->MDSteps = repetition;
+      if (repetition <= 1)  // if onyl one step, desactivate use of trajectories
+        mol->MDSteps = 0;
+      else
+        mol->MDSteps = repetition;
     } else {
       // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom)

src/config.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #endif
+#include "molecules.hpp"
+#include "periodentafel.hpp"
+#include <string>
+/****************************************** forward declarations *****************************/
+class molecule;
+class periodentafel;
+/********************************************** declarations *******************************/
 class ConfigFileBuffer {

src/datacreator.cpp

-              r87e2e39
+              rc111db
 #include "datacreator.hpp"
+#include "helpers.hpp"
+#include "parser.hpp"
 //=========================== FUNCTIONS============================
 …
 /** Plots an energy vs. order.
  * \param &Fragments EnergyMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
  * \return true if file was written successfully
  */
 bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum)
+bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
   stringstream filename;
 …
   output << "# " << msg << ", created on " << datum;
   output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
         for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
           Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+    }
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
+      for(int j=Fragments.RowCounter[ KeySets.OrderSet[BondOrder][i] ];j--;)
+        for(int k=Fragments.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];k--;)
+          Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+    }
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
     for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
 …
  * \param &Energy EnergyMatrix class containing reference values (in MatrixCounter matrix)
  * \param &Fragments EnergyMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
  * \return true if file was written successfully
  */
 bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum)
+bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
   stringstream filename;
 …
   output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
         for(int k=Fragments.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];k--;)
           Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+    }
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
+      for(int j=Fragments.RowCounter[ KeySets.OrderSet[BondOrder][i] ];j--;)
+        for(int k=Fragments.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];k--;)
+          Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+    }
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
     for (int l=0;l<Fragments.ColumnCounter[Energy.MatrixCounter];l++)
       if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
 …
 /** Plot forces vs. order.
  * \param &Fragments ForceMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
 …
  * \return true if file was written successfully
  */
 bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix,const  char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
+bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix,const  char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
+{
   stringstream filename;
 …
   output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(0.,0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, 1.);
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
     CreateForce(Fragments, Fragments.MatrixCounter);
     for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
 …
  * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
  * \param &Fragments ForceMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
 …
  * \return true if file was written successfully
  */
 bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
+bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int))
+{
   stringstream filename;
 …
   output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, -1.);
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
     CreateForce(Fragments, Fragments.MatrixCounter);
     for (int l=0;l<Fragments.ColumnCounter[Fragments.MatrixCounter];l++)
 …
  * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
  * \param &Fragments ForceMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
 …
  * \return true if file was written successfully
  */
 bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum)
+bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
   stringstream filename;
 …
   output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
     //cout << "Current order is " << BondOrder << "." << endl;
     Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, -1.);
     // errors per atom
     output << endl << "#Order\t" << BondOrder+1 << endl;
 …
 /** Plot forces error vs. vs atom vs. order.
  * \param &Fragments ForceMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
 …
  * \return true if file was written successfully
  */
 bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum)
+bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
   stringstream filename;
 …
   output << "# " << msg << ", created on " << datum;
   output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
     //cout << "Current order is " << BondOrder << "." << endl;
     Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
+    Fragments.SumSubForces(Fragments, KeySets, BondOrder, 1.);
     // errors per atom
     output << endl << "#Order\t" << BondOrder+1 << endl;
 …
  * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
  * \param &Fragments HessianMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
 …
  * \return true if file was written successfully
  */
 bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum)
+bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
   stringstream filename;
 …
   output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
   Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
     //cout << "Current order is " << BondOrder << "." << endl;
     Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.);
+    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, -1.);
     // errors per atom
     output << endl << "#Order\t" << BondOrder+1 << endl;
 …
  * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
  * \param &Fragments HessianMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
 …
  * \return true if file was written successfully
  */
 bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum)
+bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
   stringstream filename;
 …
   output << "# AtomNo\t";
   Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
     output << "Order" << BondOrder+1 << "\t";
+  }
   output << endl;
   output << Fragments.RowCounter[ Fragments.MatrixCounter ] << "\t";
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
     //cout << "Current order is " << BondOrder << "." << endl;
     Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.);
+    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, -1.);
     // frobenius norm of errors per atom
     norm = 0.;
 …
 /** Plot hessian error vs. vs atom vs. order.
  * \param &Fragments HessianMatrix class containing matrix values
  * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
+ * \param KeySets KeySetContainer class holding bond KeySetContainer::Order
  * \param *prefix prefix in filename (without ending)
  * \param *msg message to be place in first line as a comment
 …
  * \return true if file was written successfully
  */
 bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum)
+bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
+{
   stringstream filename;
 …
   output << "# AtomNo\t" << Fragments.Header[ Fragments.MatrixCounter ] << endl;
   Fragments.SetLastMatrix(0., 0);
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
     //cout << "Current order is " << BondOrder << "." << endl;
     Fragments.SumSubHessians(Fragments, KeySet, BondOrder, 1.);
+    Fragments.SumSubHessians(Fragments, KeySets, BondOrder, 1.);
     // errors per atom
     output << endl << "#Order\t" << BondOrder+1 << endl;
 …
 /** Plot matrix vs. fragment.
  */
 bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragment)(class MatrixContainer &, int))
+bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragment)(class MatrixContainer &, int))
+{
   stringstream filename;
 …
   output << "# " << msg << ", created on " << datum << endl;
   output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
       output << BondOrder+1 << "\t" << KeySet.OrderSet[BondOrder][i]+1;
       CreateFragment(Fragment, KeySet.OrderSet[BondOrder][i]);
       for (int l=0;l<Fragment.ColumnCounter[ KeySet.OrderSet[BondOrder][i] ];l++)
         output << scientific << "\t" << Fragment.Matrix[ KeySet.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySet.OrderSet[BondOrder][i] ] ][l];
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
+    for(int i=0;i<KeySets.FragmentsPerOrder[BondOrder];i++) {
+      output << BondOrder+1 << "\t" << KeySets.OrderSet[BondOrder][i]+1;
+      CreateFragment(Fragment, KeySets.OrderSet[BondOrder][i]);
+      for (int l=0;l<Fragment.ColumnCounter[ KeySets.OrderSet[BondOrder][i] ];l++)
+        output << scientific << "\t" << Fragment.Matrix[ KeySets.OrderSet[BondOrder][i] ][ Fragment.RowCounter[ KeySets.OrderSet[BondOrder][i] ] ][l];
       output << endl;
+    }
 …
 /** Copies fragment energy values into last matrix of \a Matrix with greatest total energy.
  * \param &Matrix MatrixContainer with all fragment energy values
  * \param &KeySet KeySetsContainer with associations of each fragment to a bond order
+ * \param &KeySets KeySetsContainer with associations of each fragment to a bond order
  * \param BondOrder current bond order
  */
 void CreateMaxFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, int BondOrder)
+void CreateMaxFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySets, int BondOrder)
+{
   for(int j=Fragments.RowCounter[ Fragments.MatrixCounter ];j--;) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
       if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
+    for(int i=KeySets.FragmentsPerOrder[BondOrder];i--;) {
+      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < fabs(Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][1])) {
         for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
           Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+      }
+    }
 …
 /** Copies fragment energy values into last matrix of \a Matrix with smallest total energy.
  * \param &Matrix MatrixContainer with all fragment energy values
  * \param &KeySet KeySetsContainer with associations of each fragment to a bond order
+ * \param &KeySets KeySetsContainer with associations of each fragment to a bond order
  * \param BondOrder current bond order
  */
 void CreateMinFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySet, int BondOrder)
+void CreateMinFragmentOrder(class MatrixContainer &Fragments, class KeySetsContainer &KeySets, int BondOrder)
+{
   for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
 …
     do {  // first get a minimum value unequal to 0
       for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
         Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+        Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
       i++;
     } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySet.FragmentsPerOrder[BondOrder]));
     for(;i<KeySet.FragmentsPerOrder[BondOrder];i++) { // then find lowest
       if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][1])) {
+    } while ((fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) < MYEPSILON) && (i<KeySets.FragmentsPerOrder[BondOrder]));
+    for(;i<KeySets.FragmentsPerOrder[BondOrder];i++) { // then find lowest
+      if (fabs(Fragments.Matrix[ Fragments.MatrixCounter ][j][1]) > fabs(Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][1])) {
         for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
           Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+          Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
+      }
+    }
 …
 /** Plot matrix vs. fragment.
  */
 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))
+bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))
+{
   stringstream filename;
 …
   output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
   // max
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
+  for (int BondOrder=0;BondOrder<KeySets.Order;BondOrder++) {
     Fragment.SetLastMatrix(0.,0);
     CreateFragmentOrder(Fragment, KeySet, BondOrder);
     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+    CreateFragmentOrder(Fragment, KeySets, BondOrder);
+    output << BondOrder+1 << "\t" << KeySets.FragmentsPerOrder[BondOrder];
     for (int l=0;l<Fragment.ColumnCounter[ Fragment.MatrixCounter ];l++)
       output << scientific << "\t" << Fragment.Matrix[ Fragment.MatrixCounter ][ Fragment.RowCounter[ Fragment.MatrixCounter ]-1 ][l];
 …
 /** Creates the pyxplotfile for energy data.
  * \param Matrix MatrixContainer with matrix values
  * \param KeySet contains bond order
+ * \param KeySets contains bond order
  * \param *dir directory
  * \param *prefix prefix for all filenames (without ending)
 …
  * \return true if file was written successfully
  */
 bool CreatePlotOrder(class MatrixContainer &Matrix, const class KeySetsContainer &KeySet, const char *dir, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel, const char *xrange, const char *yrange, const char *xargument, const char *uses, void (*CreatePlotLines)(ofstream &, class MatrixContainer &, const char *, const char *, const char *))
+bool CreatePlotOrder(class MatrixContainer &Matrix, const class KeySetsContainer &KeySets, const char *dir, const char *prefix, const int keycolumns, const char *key, const char *logscale, const char *extraline, const int mxtics, const int xtics, const char *xlabel, const char *ylabel, const char *xrange, const char *yrange, const char *xargument, const char *uses, void (*CreatePlotLines)(ofstream &, class MatrixContainer &, const char *, const char *, const char *))
+{
   stringstream filename;

src/datacreator.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+//============================ INCLUDES ===========================
+#include <iostream>
+#include "helpers.hpp"
+#include "parser.hpp"
+/****************************************** forward declarations *****************************/
+//=========================== FUNCTIONS============================
+class EnergyMatrix;
+class ForceMatrix;
+class HessianMatrix;
+class KeySetsContainer;
+class MatrixContainer;
+/********************************************** declarations *******************************/
 bool OpenOutputFile(ofstream &output, const char *dir, const char *filename);

src/element.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #endif
 #include <iostream>
 #include "defs.hpp"
+/********************************************** declarations *******************************/
 /** Chemical element.

src/ellipsoid.cpp

-              r87e2e39
+              rc111db
 #include <gsl/gsl_vector.h>
+#include <iomanip>
+#include <set>
 #include "boundary.hpp"
 #include "ellipsoid.hpp"
+#include "linkedcell.hpp"
+#include "tesselation.hpp"
+#include "vector.hpp"
+#include "verbose.hpp"
 /** Determines squared distance for a given point \a x to surface of ellipsoid.

src/ellipsoid.hpp

-              r87e2e39
+              rc111db
 #define ELLIPSOID_HPP_
+using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #endif
+#include "linkedcell.hpp"
+#include "vector.hpp"
+#include <gsl/gsl_vector.h>
+/****************************************** forward declarations *****************************/
+class LinkedCell;
+class Vector;
+/********************************************** declarations *******************************/
 double SquaredDistanceToEllipsoid(Vector &x, Vector &EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle);

src/graph.cpp

-              r87e2e39
+              rc111db
 using namespace std;
+#include "atom.hpp"
+#include "bond.hpp"
+#include "config.hpp"
 #include "graph.hpp"
+#include "molecule.hpp"
 /***************************************** Implementations for graph classes ********************************/
+/** Constructor of class Graph.
+///** Constructor of class Graph.
+// */
+//Graph::Graph()
+//{
+//};
+//
+///** Destructor of class Graph.
+// * Destructor does release memory for nodes and edges contained in its lists as well.
+// */
+//Graph::~Graph()
+//{
+//};
+//
+///** Constructor of class SubGraph.
+// */
+//SubGraph::SubGraph()
+//{
+//};
+//
+///** Destructor of class SubGraph.
+// * Note that destructor does not deallocate either nodes or edges! (this is done by its subgraph!)
+// */
+//SubGraph::~SubGraph()
+//{
+//};
+//
+///** Constructor of class Node.
+// */
+//Node::Node()
+//{
+//};
+//
+///** Destructor of class Node.
+// */
+//Node::~Node()
+//{
+//};
+//
+///** Constructor of class Edge.
+// */
+//Edge::Edge()
+//{
+//};
+//
+///** Destructor of class Edge.
+// */
+//Edge::~Edge()
+//{
+//};
+//bool operator < (KeySet SubgraphA, KeySet SubgraphB)
+//{
+//  return KeyCompare(SubgraphA, SubgraphB);
+//};
+/** Checking whether KeySet is not already present in Graph, if so just adds factor.
+ * \param *out output stream for debugging
+ * \param &set KeySet to insert
+ * \param &graph Graph to insert into
+ * \param *counter pointer to unique fragment count
+ * \param factor energy factor for the fragment
  */
 Graph::Graph()
+void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment)
+{
+  GraphTestPair testGraphInsert;
+  testGraphInsert = Fragment->Leaflet->insert(GraphPair (*Fragment->FragmentSet,pair<int,double>(Fragment->FragmentCounter,Fragment->TEFactor)));  // store fragment number and current factor
+  if (testGraphInsert.second) {
+    *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " successfully inserted." << endl;
+    Fragment->FragmentCounter++;
+  } else {
+    *out << Verbose(2) << "KeySet " << Fragment->FragmentCounter << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
+    ((*(testGraphInsert.first)).second).second += Fragment->TEFactor;  // increase the "created" counter
+    *out << Verbose(2) << "New factor is " << ((*(testGraphInsert.first)).second).second << "." << endl;
+  }
+};
+//void inline InsertIntoGraph(ofstream *out, KeyStack &stack, Graph &graph, int *counter, double factor)
+//{
+//  // copy stack contents to set and call overloaded function again
+//  KeySet set;
+//  for(KeyStack::iterator runner = stack.begin(); runner != stack.begin(); runner++)
+//    set.insert((*runner));
+//  InsertIntoGraph(out, set, graph, counter, factor);
+//};
+/** Inserts each KeySet in \a graph2 into \a graph1.
+ * \param *out output stream for debugging
+ * \param graph1 first (dest) graph
+ * \param graph2 second (source) graph
+ * \param *counter keyset counter that gets increased
+ */
+void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
+{
+  GraphTestPair testGraphInsert;
+  for(Graph::iterator runner = graph2.begin(); runner != graph2.end(); runner++) {
+    testGraphInsert = graph1.insert(GraphPair ((*runner).first,pair<int,double>((*counter)++,((*runner).second).second)));  // store fragment number and current factor
+    if (testGraphInsert.second) {
+      *out << Verbose(2) << "KeySet " << (*counter)-1 << " successfully inserted." << endl;
+    } else {
+      *out << Verbose(2) << "KeySet " << (*counter)-1 << " failed to insert, present fragment is " << ((*(testGraphInsert.first)).second).first << endl;
+      ((*(testGraphInsert.first)).second).second += (*runner).second.second;
+      *out << Verbose(2) << "New factor is " << (*(testGraphInsert.first)).second.second << "." << endl;
+    }
+  }
 };
-/** Destructor of class Graph.
- * Destructor does release memory for nodes and edges contained in its lists as well.
- */
-Graph::~Graph()
+{
-};
-/** Constructor of class SubGraph.
- */
-SubGraph::SubGraph()
+{
-};
-/** Destructor of class SubGraph.
- * Note that destructor does not deallocate either nodes or edges! (this is done by its subgraph!)
- */
-SubGraph::~SubGraph()
+{
-};
-/** Constructor of class Node.
- */
-Node::Node()
+{
-};
-/** Destructor of class Node.
- */
-Node::~Node()
+{
-};
-/** Constructor of class Edge.
- */
-Edge::Edge()
+{
-};
-/** Destructor of class Edge.
- */
-Edge::~Edge()
+{
-};

src/graph.hpp

-              r87e2e39
+              rc111db
 #define GRAPH_HPP_
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 // STL headers
+#include <deque>
 #include <map>
 #include <multimap>
+#include <set>
+#include "molecules.hpp"
+/****************************************** forward declarations *****************************/
+class atom;
+class bond;
+class config;
+class molecule;
 class Graph;
 …
 class Edge;
 /***************************************** Various graph-related STL defines ********************************/
+/********************************************** definitions *********************************/
 #define NodeMap pair < int, class Node* >
 #define EdgeMap multimap < class Node*, class Edge* >
+/***************************************** Definition for classes ********************************/
+#define KeyStack deque<int>
+#define KeySet set<int>
+#define NumberValuePair pair<int, double>
+#define Graph map <KeySet, NumberValuePair, KeyCompare >
+#define GraphPair pair <KeySet, NumberValuePair >
+#define KeySetTestPair pair<KeySet::iterator, bool>
+#define GraphTestPair pair<Graph::iterator, bool>
+/** Graph class containing the graphs behind molecules.
+ */
+class Graph
+/******************************** Some small functions and/or structures **********************************/
+struct KeyCompare
+{
+        NodeMap ListOfNodes;            //!< tree-list of all nodes in this graph
+        EdgeMap ListOfEdges;            //!< tree-multi-list of all nodes, referenced to node id
+  bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
 };
+/** Class describing subgraphs of the Class \a Graph.
+ * SubGraph has its own node and edge lists, however also a pointer to its father graph
+ * and hence access to its list as well.
+//bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)
+void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment); // Insert a KeySet into a Graph
+void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);  // Insert all KeySet's in a Graph into another Graph
+/** Structure containing all values in power set combination generation.
  */
+class SubGraph : class Graph
+{
+        class Graph *FatherGraph;               //!< Graph whose subgraph we are
+struct UniqueFragments {
+  config *configuration;
+  atom *Root;
+  Graph *Leaflet;
+  KeySet *FragmentSet;
+  int ANOVAOrder;
+  int FragmentCounter;
+  int CurrentIndex;
+  double TEFactor;
+  int *ShortestPathList;
+  bool **UsedList;
+  bond **BondsPerSPList;
+  int *BondsPerSPCount;
 };
+/** Class containing the nodes of a graph.
+ */
+class Node
+{
+        int     id;                             //!< individual id of the node
+        char *Name;             //!< Name of the node for pretty printing
+};
+/********************************************** declarations *******************************/
+/** Class containing egdes in a Graph strructure.
+ */
+class Edge
+{
+        class Node *leftnode;           //!< pointer to first node
+        class Node *atomnode;           //!< pointer to second node
+};
+///** Graph class containing the graphs behind molecules.
+// */
+//class Graph
+//{
+//      NodeMap ListOfNodes;            //!< tree-list of all nodes in this graph
+//      EdgeMap ListOfEdges;            //!< tree-multi-list of all nodes, referenced to node id
+//};
+//
+///** Class describing subgraphs of the Class \a Graph.
+// * SubGraph has its own node and edge lists, however also a pointer to its father graph
+// * and hence access to its list as well.
+// */
+//class SubGraph : class Graph
+//{
+//      class Graph *FatherGraph;               //!< Graph whose subgraph we are
+//};
+//
+///** Class containing the nodes of a graph.
+// */
+//class Node
+//{
+//      int     id;                             //!< individual id of the node
+//      char *Name;             //!< Name of the node for pretty printing
+//};
+//
+///** Class containing egdes in a Graph strructure.
+// */
+//class Edge
+//{
+//      class Node *leftnode;           //!< pointer to first node
+//      class Node *atomnode;           //!< pointer to second node
+//};

src/helpers.cpp

-              r87e2e39
+              rc111db
 #include "helpers.hpp"
-#include "memoryusageobserver.hpp"
 /********************************************** helpful functions *********************************/
 …
 };
+/**
+ * Allocates a memory range using malloc().
+/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
+ * \param *symm 6-dim array of unique symmetric matrix components
+ * \return allocated NDIM*NDIM array with the symmetric matrix
+ */
+double * ReturnFullMatrixforSymmetric(double *symm)
+{
+  double *matrix = Malloc<double>(NDIM * NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
+  matrix[0] = symm[0];
+  matrix[1] = symm[1];
+  matrix[2] = symm[3];
+  matrix[3] = symm[1];
+  matrix[4] = symm[2];
+  matrix[5] = symm[4];
+  matrix[6] = symm[3];
+  matrix[7] = symm[4];
+  matrix[8] = symm[5];
+  return matrix;
+};
+/** Comparison function for GSL heapsort on distances in two molecules.
+ * \param *a
+ * \param *b
+ * \return <0, \a *a less than \a *b, ==0 if equal, >0 \a *a greater than \a *b
+ */
+int CompareDoubles (const void * a, const void * b)
+{
+  if (*(double *)a > *(double *)b)
+    return -1;
+  else if (*(double *)a < *(double *)b)
+    return 1;
+  else
+    return 0;
+};
+/** Allocates a memory range using malloc().
  * Prints the provided error message in case of a failure.
+ *
 …
  * Frees all memory registered by the memory observer and calls exit(225) afterwards.
  */
 void performCriticalExit() {
+static void performCriticalExit() {
   map<void*, size_t> pointers = MemoryUsageObserver::getInstance()->getPointersToAllocatedMemory();
   for (map<void*, size_t>::iterator runner = pointers.begin(); runner != pointers.end(); runner++) {

src/helpers.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #endif
-#include <iostream>
-#include <iomanip>
 #include <fstream>
-#include <sstream>
-#include <math.h>
-#include <string>
-#include "defs.hpp"
-#include "verbose.hpp"
 #include "memoryallocator.hpp"
 …
 char *FixedDigitNumber(const int FragmentNumber, const int digits);
 bool IsValidNumber( const char *string);
+int CompareDoubles (const void * a, const void * b);
+double * ReturnFullMatrixforSymmetric(double *cell_size);
 static void performCriticalExit();
 …
 };
-/******************************** Some templates for list management ***********************************/
-/** Adds linking of an item to a list.
- * \param *walker
- * \return true - adding succeeded, false - error in list
- */
-template <typename X> void link(X *walker, X *end)
+{
-  X *vorher = end->previous;
-  if (vorher != NULL)
-    vorher->next = walker;
-  end->previous = walker;
-  walker->previous = vorher;
-  walker->next = end;
-};
-/** Removes linking of an item in a list.
- * \param *walker
- * \return true - removing succeeded, false - given item not found in list
- */
-template <typename X> void unlink(X *walker)
+{
-  if (walker->next != NULL)
-    walker->next->previous = walker->previous;
-  if (walker->previous != NULL)
-    walker->previous->next = walker->next;
-};
-/** Adds new item before an item \a *end in a list.
- * \param *pointer   item to be added
- * \param *end  end of list
- * \return true - addition succeeded, false - unable to add item to list
- */
-template <typename X>  bool add(X *pointer, X *end)
+{
-  if (end != NULL) {
-    link(pointer, end);
-  } else {
-    pointer->previous = NULL;
-    pointer->next = NULL;
+  }
-  return true;
-};
-/** Finds item in list
- * \param *suche  search criteria
- * \param *start  begin of list
- * \param *end  end of list
- * \return X - if found, NULL - if not found
- */
-template <typename X, typename Y> X * find(Y *suche, X *start, X *end)
+{
-  X *walker = start;
-  while (walker->next != end) { // go through list
-    walker = walker->next; // step onward beforehand
-    if (*walker->sort == *suche) return (walker);
+  }
-  return NULL;
-};
-/** Removes an item from the list without check.
- * \param *walker item to be removed
- * \return true - removing succeeded, false - given item not found in list
- */
-template <typename X> void removewithoutcheck(X *walker)
+{
-  if (walker != NULL) {
-    unlink(walker);
-    delete(walker);
-    walker = NULL;
+  }
-};
-/** Removes an item from the list, checks if exists.
- * Checks beforehand if atom is really within molecule list.
- * \param *pointer   item to be removed
- * \param *start  begin of list
- * \param *end  end of list
- * \return true - removing succeeded, false - given item not found in list
- */
-template <typename X> bool remove(X *pointer, X *start, X *end)
+{
-  X *walker = find (pointer->sort, start, end);
-/*  while (walker->next != pointer) { // search through list
-    walker = walker->next;
-    if (walker == end) return false;  // item not found in list
-  }*/
-  // atom found, now unlink
-  if (walker != NULL)
-    removewithoutcheck(walker);
-  else
-    return false;
-  return true;
-};
-/** Cleans the whole list.
- * \param *start begin of list
- * \param *end end of list
- * \return true - list was cleaned successfully, false - error in list structure
- */
-template <typename X> bool cleanup(X *start, X *end)
+{
-  X *pointer = start->next;
-  X *walker;
-  while (pointer != end) { // go through list
-    walker = pointer; // mark current
-    pointer = pointer->next; // step onward beforehand
-    // remove walker
-    unlink(walker);
-    delete(walker);
-    walker = NULL;
+  }
-  return true;
-};
-/** Returns the first marker in a chain list.
- * \param *me one arbitrary item in chain list
- * \return poiner to first marker
- */
-template <typename X> X *GetFirst(X *me)
+{
-  X *Binder = me;
-  while(Binder->previous != NULL)
-    Binder = Binder->previous;
-  return Binder;
-};
-/** Returns the last marker in a chain list.
- * \param *me one arbitrary item in chain list
- * \return poiner to last marker
- */
-template <typename X> X *GetLast(X *me)
+{
-  X *Binder = me;
-  while(Binder->next != NULL)
-    Binder = Binder->next;
-  return Binder;
-};
 /** Frees a two-dimensional array.
  * \param *ptr pointer to array
 …
 };
+template <typename T> void Increment(T *value, T inc)
+{
+  *value += inc;
+};
+template <typename T> void AbsoluteValue(T *value, T abs)
+{
+  *value = abs;
+};

src/leastsquaremin.cpp

-              r87e2e39
+              rc111db
  */
+#include <iostream>
 #include "leastsquaremin.hpp"
+#include "vector.hpp"
 /** Determines sum of squared distances of \a X to all \a **vectors.

src/leastsquaremin.hpp

-              r87e2e39
+              rc111db
 #define LEASTSQUAREMIN_HPP_
+#include "vector.hpp"
+using namespace std;
+/*********************************************** includes ***********************************/
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <gsl/gsl_vector.h>
+/****************************************** forward declarations *****************************/
 class element;
 class molecule;
+class Vector;
+/********************************************** declarations *******************************/
 /** Parameter structure for least square minimsation.

src/linkedcell.cpp

-              r87e2e39
+              rc111db
+#include "atom.hpp"
+#include "helpers.hpp"
 #include "linkedcell.hpp"
 #include "molecules.hpp"
+#include "molecule.hpp"
 #include "tesselation.hpp"
+#include "vector.hpp"
 // ========================================================= class LinkedCell ===========================================
 …
     //cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
     set->GoToNext();
+  }
+  cout << "done."  << endl;
+  cout << Verbose(1) << "End of LinkedCell" << endl;
+};
+/** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
+ * \param *set LCNodeSet class with all LCNode's
+ * \param RADIUS edge length of cells
+ */
+LinkedCell::LinkedCell(LinkedNodes *set, double radius)
+{
+  class TesselPoint *Walker = NULL;
+  RADIUS = radius;
+  LC = NULL;
+  for(int i=0;i<NDIM;i++)
+    N[i] = 0;
+  index = -1;
+  max.Zero();
+  min.Zero();
+  cout << Verbose(1) << "Begin of LinkedCell" << endl;
+  if (set->empty()) {
+    cerr << "ERROR: set contains no linked cell nodes!" << endl;
+    return;
+  }
+  // 1. find max and min per axis of atoms
+  LinkedNodes::iterator Runner = set->begin();
+  for (int i=0;i<NDIM;i++) {
+    max.x[i] = (*Runner)->node->x[i];
+    min.x[i] = (*Runner)->node->x[i];
+  }
+  for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
+    Walker = *Runner;
+    for (int i=0;i<NDIM;i++) {
+      if (max.x[i] < Walker->node->x[i])
+        max.x[i] = Walker->node->x[i];
+      if (min.x[i] > Walker->node->x[i])
+        min.x[i] = Walker->node->x[i];
+    }
+  }
+  cout << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
+  // 2. find then number of cells per axis
+  for (int i=0;i<NDIM;i++) {
+    N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
+  }
+  cout << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl;
+  // 3. allocate the lists
+  cout << Verbose(2) << "Allocating cells ... ";
+  if (LC != NULL) {
+    cout << Verbose(1) << "ERROR: Linked Cell list is already allocated, I do nothing." << endl;
+    return;
+  }
+  LC = new LinkedNodes[N[0]*N[1]*N[2]];
+  for (index=0;index<N[0]*N[1]*N[2];index++) {
+    LC [index].clear();
+  }
+  cout << "done."  << endl;
+  // 4. put each atom into its respective cell
+  cout << Verbose(2) << "Filling cells ... ";
+  for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
+    Walker = *Runner;
+    for (int i=0;i<NDIM;i++) {
+      n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
+    }
+    index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
+    LC[index].push_back(Walker);
+    //cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
+  }
   cout << "done."  << endl;

src/linkedcell.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #include "defs.hpp"
-#include "helpers.hpp"
 #include "vector.hpp"
+/****************************************** forward declarations *****************************/
+class PointCloud;
 class TesselPoint;
+class PointCloud;
+/********************************************** definitions *********************************/
 #define LinkedNodes list<TesselPoint *>
+/********************************************** declarations *******************************/
 /** Linked Cell class for containing Vectors in real space efficiently.
 …
     LinkedCell();
     LinkedCell(PointCloud *set, double RADIUS);
+    LinkedCell(LinkedNodes *set, double radius);
     ~LinkedCell();
     LinkedNodes* GetCurrentCell();

src/memoryallocator.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 …
 #include "memoryusageobserver.hpp"
 /******************* wrappers for memory allocation functions ***********************/
+/********************************************** declarations *******************************/
+/**
+ * Allocates a memory range using malloc().
+/** Allocates a memory range using malloc().
  * Prints the provided error message in case of a failure.
+ *
 …
 template <> char* Malloc<char>(size_t size, const char* output);
+/**
+ * Allocates a memory range using calloc().
+/** Allocates a memory range using calloc().
  * Prints the provided error message in case of a failure.
+ *
 …
 };
+/**
+ * Reallocates a memory range using realloc(). If the provided pointer to the old
+/** Reallocates a memory range using realloc(). If the provided pointer to the old
  * memory range is NULL, malloc() is called instead.
  * Prints the provided error message in case of a failure (of either malloc() or realloc()).
 …
 };
+/**
+ * Frees allocated memory range using free().
+/** Frees allocated memory range using free().
+ *
  * \param pointer to the allocated memory range to free; may be NULL, this function is a no-op then
+ * \param *msg optional error message
  */
 template <typename X> void Free(X** buffer, const char *msg = NULL)

src/memoryusageobserver.cpp

-              r87e2e39
+              rc111db
+ *
  * \param pointer to the allocated piece of memory
+ * \param *msg optional error message
  */
 void MemoryUsageObserver::removeMemory(void* pointer, const char *msg) {
 …
   if (current == memoryUsers.end()) {
     cout << "WARNING: There is non-tracked memory to be freed. Pointer "
       << pointer << " is not registered by MemoryUsageObserver";
+      << pointer << " is not registered by MemoryUsageObserver: ";
     if (msg != NULL)
       cout << ": " << msg;
     cout << "." << endl;
+      cout << *msg;
+    cout << endl;
     return;
+  }

src/memoryusageobserver.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 …
 #include <string>
 #include <typeinfo>
+/********************************************** declarations *******************************/
 class MemoryUsageObserver {

src/moleculelist.cpp

-              r87e2e39
+              rc111db
  */
+#include "atom.hpp"
+#include "bond.hpp"
+#include "boundary.hpp"
 #include "config.hpp"
+#include "molecules.hpp"
+#include "element.hpp"
+#include "helpers.hpp"
+#include "linkedcell.hpp"
+#include "molecule.hpp"
 #include "memoryallocator.hpp"
+#include "periodentafel.hpp"
 /*********************************** Functions for class MoleculeListClass *************************/
 …
 /** Embedding merge of a given set of molecules into one.
  * Embedding merge inserts one molecule into the other.
  * \param *mol destination molecule
  * \param *srcmol source molecule
  * \return true - merge successful, false - merge failed (probably due to non-existant indices
  * \TODO find embedding center
+ * \param *mol destination molecule (fixed one)
+ * \param *srcmol source molecule (variable one, where atoms are taken from)
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices)
+ * \TODO linked cell dimensions for boundary points has to be as big as inner diameter!
  */
 bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
+{
+  if (srcmol == NULL)
+  if ((srcmol == NULL) || (mol == NULL)) {
+    cout << Verbose(1) << "ERROR: Either fixed or variable molecule is given as NULL." << endl;
     return false;
+  // calculate center for merge
+  srcmol->Center.CopyVector(mol->FindEmbeddingHole((ofstream *)&cout, srcmol));
+  srcmol->Center.Zero();
+  // perform simple merge
+  SimpleMerge(mol, srcmol);
+  }
+  // calculate envelope for *mol
+  LinkedCell *LCList = new LinkedCell(mol, 8.);
+  FindNonConvexBorder((ofstream *)&cout, mol, LCList, 4., NULL);
+  if (mol->TesselStruct == NULL) {
+    cout << Verbose(1) << "ERROR: Could not tesselate the fixed molecule." << endl;
+    return false;
+  }
+  delete(LCList);
+  LCList = new LinkedCell(mol->TesselStruct, 8.);  // re-create with boundary points only!
+  // prepare index list for bonds
+  srcmol->CountAtoms((ofstream *)&cout);
+  atom ** CopyAtoms = new atom*[srcmol->AtomCount];
+  for(int i=0;i<srcmol->AtomCount;i++)
+    CopyAtoms[i] = NULL;
+  // for each of the source atoms check whether we are in- or outside and add copy atom
+  atom *Walker = srcmol->start;
+  int nr=0;
+  while (Walker->next != srcmol->end) {
+    Walker = Walker->next;
+    cout << Verbose(2) << "INFO: Current Walker is " << *Walker << "." << endl;
+    if (!mol->TesselStruct->IsInnerPoint((ofstream *)&cout, Walker->x, LCList)) {
+      CopyAtoms[Walker->nr] = new atom(Walker);
+      mol->AddAtom(CopyAtoms[Walker->nr]);
+      nr++;
+    } else {
+      // do nothing
+    }
+  }
+  cout << Verbose(1) << nr << " of " << srcmol->AtomCount << " atoms have been merged.";
+  // go through all bonds and add as well
+  bond *Binder = srcmol->first;
+  while(Binder->next != srcmol->last) {
+    Binder = Binder->next;
+    cout << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
+    mol->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
+  }
+  delete(LCList);
   return true;
 };

src/parser.cpp

-              r87e2e39
+              rc111db
  * Sums over "E"-terms to create the "F"-terms
  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
  * \return true if summing was successful
  */
 bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
+bool MatrixContainer::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
+{
   // go through each order
   for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
+  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
+    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
     // then go per order through each suborder and pick together all the terms that contain this fragment
     for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
       for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
         if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
           //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
+      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
+        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
+          //cout << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
           // if the fragment's indices are all in the current fragment
           for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
             int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
+          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
+            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
             //cout << "Current index is " << k << "/" << m << "." << endl;
             if (m != -1) { // if it's not an added hydrogen
               for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
                 //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
                 if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
+              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
+                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
+                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                   m = l;
                   break;
 …
+              }
               //cout << "Corresponding index in CurrentFragment is " << m << "." << endl;
               if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
                 cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
+                cerr << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current force index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
                 return false;
+              }
               if (Order == SubOrder) { // equal order is always copy from Energies
                 for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
+                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
               } else {
                 for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
                   Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                for(int l=ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l--;)
+                  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
+              }
+            }
             //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
              //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
+             //cout << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+          }
         } else {
           //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
+          //cout << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
+        }
+      }
+    }
    //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
+   //cout << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
+  }
 …
  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
  * \parsm sign +1 or -1
  * \return true if summing was successful
  */
 bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign)
+bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySets, int Order, double sign)
+{
   // sum energy
   if (CorrectionFragments == NULL)
     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
         for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
           Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
+    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
+        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
+          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k];
   else
     for(int i=KeySet.FragmentsPerOrder[Order];i--;)
       for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
         for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
           Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
+    for(int i=KeySets.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySets.OrderSet[Order][i] ];j--;)
+        for(int k=ColumnCounter[ KeySets.OrderSet[Order][i] ];k--;)
+          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySets.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySets.OrderSet[Order][i] ][j][k]);
   return true;
 };
 …
 /** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
  *  \param sign +1 or -1
  * \return true if summing was successful
  */
 bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
+bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
+{
   int FragmentNr;
   // sum forces
   for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
     FragmentNr = KeySet.OrderSet[Order][i];
+  for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
+    FragmentNr = KeySets.OrderSet[Order][i];
     for(int l=0;l<RowCounter[ FragmentNr ];l++) {
       int j = Indices[ FragmentNr ][l];
 …
 /** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
  *  \param sign +1 or -1
  * \return true if summing was successful
  */
 bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
+bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
+{
   int FragmentNr;
   // sum forces
   for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
     FragmentNr = KeySet.OrderSet[Order][i];
+  for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
+    FragmentNr = KeySets.OrderSet[Order][i];
     for(int l=0;l<RowCounter[ FragmentNr ];l++) {
       int j = Indices[ FragmentNr ][l];
 …
  * Sums over "E"-terms to create the "F"-terms
  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
+ * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
  * \return true if summing was successful
  */
 bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
+bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
+{
   // go through each order
   for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
+  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
+    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
     // then go per order through each suborder and pick together all the terms that contain this fragment
     for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
       for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
         if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
           //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
+      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
+        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
+          //cout << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
           // if the fragment's indices are all in the current fragment
           for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
             int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
+          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
+            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
             //cout << "Current row index is " << k << "/" << m << "." << endl;
             if (m != -1) { // if it's not an added hydrogen
               for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
                 //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
                 if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
+              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
+                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
+                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                   m = l;
                   break;
 …
+              }
               //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
               if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
                 cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
+                cerr << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
                 return false;
+              }
               for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) {
                 int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
+              for(int l=0;l<ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l++) {
+                int n = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][l];
                 //cout << "Current column index is " << l << "/" << n << "." << endl;
                 if (n != -1) { // if it's not an added hydrogen
                   for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
                     //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
                     if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
+                  for (int p=0;p<ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
+                    //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
+                    if (n == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p]) {
                       n = p;
                       break;
 …
+                  }
                   //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
                   if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
                     cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
+                  if (n > ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
+                    cerr << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl;
                     return false;
+                  }
                   if (Order == SubOrder) { // equal order is always copy from Energies
                     //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                    //cout << "Adding " << MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
                   } else {
                     //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                     Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
+                    //cout << "Subtracting " << Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
+                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
+                  }
+                }
+              }
+            }
             //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
              //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
+             //cout << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+          }
         } else {
           //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
+          //cout << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
+        }
+      }
+    }
    //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
+   //cout << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
+  }

src/parser.hpp

-              r87e2e39
+              rc111db
 /** \file parsing.hpp
+/** \file parser.hpp
+ *
  * Definitions of various class functions for the parsing of value files.
 …
 #ifndef PARSING_HPP_
 #define PARSING_HPP_
+#ifndef PARSER_HPP_
+#define PARSER_HPP_
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 …
 #endif
+// ======================================= DEFINES ==========================================
+/****************************************** forward declarations *****************************/
+class EnergyMatrix;
+class ForceMatrix;
+class HessianMatrix;
+class KeySetsContainer;
+class MatrixContainer;
+/********************************************** definitions *********************************/
 #define EnergySuffix ".energy.all"
 …
 bool TestParams(int argc, char **argv);
+// ======================================= CLASS KeySetsContainer =============================
+class KeySetsContainer {
+  public:
+    int **KeySets;
+    int *AtomCounter;
+    int FragmentCounter;
+    int Order;
+    int *FragmentsPerOrder;
+    int **OrderSet;
+  KeySetsContainer();
+  ~KeySetsContainer();
+  bool ParseKeySets(const char *name, const int *ACounter, const int FCounter);
+  bool ParseManyBodyTerms();
+  bool Contains(const int GreaterSet, const int SmallerSet);
+};
 // ======================================= CLASS MatrixContainer =============================
 …
   //bool ParseIndices();
   //bool SumSubValues();
   bool SumSubManyBodyTerms(class MatrixContainer &Matrix, class KeySetsContainer &KeySet, int Order);
+  bool SumSubManyBodyTerms(class MatrixContainer &Matrix, class KeySetsContainer &KeySets, int Order);
   bool WriteTotalFragments(const char *name, const char *prefix);
   bool WriteLastMatrix(const char *name, const char *prefix, const char *suffix);
 …
   public:
     bool ParseIndices();
     bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign);
+    bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySets, int Order, double sign);
     bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
 };
 …
   public:
     bool ParseIndices(const char *name);
     bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
+    bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign);
     bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
 };
 …
     //~HessianMatrix();
     bool ParseIndices(char *name);
     bool SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order);
     bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
+    bool SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order);
+    bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign);
     bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
   private:
 …
 };
-// ======================================= CLASS KeySetsContainer =============================
-class KeySetsContainer {
-  public:
-    int **KeySets;
-    int *AtomCounter;
-    int FragmentCounter;
-    int Order;
-    int *FragmentsPerOrder;
-    int **OrderSet;
-  KeySetsContainer();
-  ~KeySetsContainer();
-  bool ParseKeySets(const char *name, const int *ACounter, const int FCounter);
-  bool ParseManyBodyTerms();
-  bool Contains(const int GreaterSet, const int SmallerSet);
-};
 // ======================================= END =============================================
 #endif /*PARSING_HPP_*/
+#endif /*PARSER_HPP_*/

src/periodentafel.cpp

-              r87e2e39
+              rc111db
 #include <fstream>
+#include "element.hpp"
 #include "helpers.hpp"
+#include "lists.hpp"
 #include "periodentafel.hpp"
 #include "verbose.hpp"

src/periodentafel.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 …
 #include "defs.hpp"
-#include "element.hpp"
+// ====================================== class definitions =========================
+/****************************************** forward declarations *****************************/
+class element;
+/********************************************** declarations *******************************/

src/stackclass.hpp

-              r87e2e39
+              rc111db
 #define STACKCLASS_HPP_
+using namespace std;
+/*********************************************** includes ***********************************/
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 #include "verbose.hpp"
 #include "memoryallocator.hpp"
+/****************************************** forward declarations *****************************/
 template <typename T> class StackClass;

src/tesselation.cpp

-              r87e2e39
+              rc111db
  */
+#include <fstream>
+#include "linkedcell.hpp"
 #include "tesselation.hpp"
+#include "memoryallocator.hpp"
+#include "tesselationhelpers.hpp"
+#include "vector.hpp"
+#include "verbose.hpp"
+class molecule;
 // ======================================== Points on Boundary =================================
 …
 BoundaryPointSet::~BoundaryPointSet()
+{
   cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  //cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
   if (!lines.empty())
     cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
 …
 ostream & operator <<(ostream &ost, BoundaryPointSet &a)
+{
   ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.node->Name << " at " << *a.node->node << "]";
   return ost;
+}
 …
         for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
           if ((*Runner).second == this) {
             cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+            //cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
             endpoints[i]->lines.erase(Runner);
             break;
+          }
       } else { // there's just a single line left
+        if (endpoints[i]->lines.erase(Nr))
+          cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+        if (endpoints[i]->lines.erase(Nr)) {
+          //cout << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+        }
+      }
       if (endpoints[i]->lines.empty()) {
         cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+        //cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         if (endpoints[i] != NULL) {
           delete(endpoints[i]);
 …
 /** 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.
+ * We sum the two angles of each height vector with respect to the center of the baseline.
  * If greater/equal M_PI than we are convex.
  * \param *out output stream for debugging
 …
   // get the two triangles
   if (triangles.size() != 2) {
     *out << Verbose(1) << "ERROR: Baseline " << *this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
+    *out << Verbose(1) << "ERROR: Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;
     return true;
+  }
 …
     NormalCheck.Scale(sign);
     sign = -sign;
+    BaseLineNormal.SubtractVector(&runner->second->NormalVector);   // we subtract as BaseLineNormal has to point inward in direction of [pi,2pi]
+    if (runner->second->NormalVector.NormSquared() > MYEPSILON)
+      BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
+    else {
+      *out << Verbose(1) << "CRITICAL: Triangle " << *runner->second << " has zero normal vector!" << endl;
+      exit(255);
+    }
     node = runner->second->GetThirdEndpoint(this);
     if (node != NULL) {
 …
       i++;
     } else {
       //*out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
+      //*out << Verbose(2) << "ERROR: I cannot find third node in triangle, something's wrong." << endl;
       return true;
+    }
 …
   //*out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
   if (NormalCheck.NormSquared() < MYEPSILON) {
     *out << Verbose(2) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
+    *out << Verbose(3) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
     return true;
+  }
+  BaseLineNormal.Scale(-1.);
   double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
   if ((angle - M_PI) > -MYEPSILON) {
     *out << Verbose(2) << "ACCEPT: Angle is greater than pi: convex." << endl;
+    *out << Verbose(3) << "ACCEPT: Angle is greater than pi: convex." << endl;
     return true;
   } else {
     *out << Verbose(2) << "REJECT: Angle is less than pi: concave." << endl;
+    *out << Verbose(3) << "REJECT: Angle is less than pi: concave." << endl;
     return false;
+  }
 …
 ostream & operator <<(ostream &ost, BoundaryLineSet &a)
+{
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]";
   return ost;
 };
 …
   for (int i = 0; i < 3; i++) {
     if (lines[i] != NULL) {
+      if (lines[i]->triangles.erase(Nr))
+        cout << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      if (lines[i]->triangles.erase(Nr)) {
+        //cout << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;
+      }
       if (lines[i]->triangles.empty()) {
           cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+          //cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
           delete (lines[i]);
           lines[i] = NULL;
 …
+    }
+  }
   cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+  //cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
 };
 …
  * 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.
+ * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
+ * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
+ * the first two basepoints) or not.
  * \param *out output stream for debugging
  * \param *MolCenter offset vector of line
 …
     exit(255);
+  }
   CrossPoint.SubtractVector(endpoints[i%3]->node->node);
+  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
   // check whether intersection is inside or not by comparing length of intersection and length of cross point
   if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
     return true;
   } else { // outside!
 …
   for(int i=0;i<3;i++)
     if (point == endpoints[i])
+      return true;
+  return false;
+};
+/** Checks whether point is any of the three endpoints this triangle contains.
+ * \param *point TesselPoint to test
+ * \return true - point is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
+{
+  for(int i=0;i<3;i++)
+    if (point == endpoints[i]->node)
       return true;
   return false;
 …
 };
+/** Checks whether three given \a *Points coincide with triangle's endpoints.
+ * \param *Points[3] pointer to BoundaryPointSet
+ * \return true - is the very triangle, false - is not
+ */
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
+{
+  return (((endpoints[0] == T->endpoints[0])
+            || (endpoints[0] == T->endpoints[1])
+            || (endpoints[0] == T->endpoints[2])
+          ) && (
+            (endpoints[1] == T->endpoints[0])
+            || (endpoints[1] == T->endpoints[1])
+            || (endpoints[1] == T->endpoints[2])
+          ) && (
+            (endpoints[2] == T->endpoints[0])
+            || (endpoints[2] == T->endpoints[1])
+            || (endpoints[2] == T->endpoints[2])
+          ));
+};
 /** Returns the endpoint which is not contained in the given \a *line.
  * \param *line baseline defining two endpoints
 …
 ostream &operator <<(ostream &ost, BoundaryTriangleSet &a)
+{
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
       << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
+      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
   return ost;
 };
 …
 TesselPoint::~TesselPoint()
+{
-  Free(&Name);
 };
 …
 ostream & operator << (ostream &ost, const TesselPoint &a)
+{
   ost << "[" << (a.Name) << "|" << &a << "]";
+  ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
   return ost;
 };
 …
   TrianglesOnBoundaryCount = 0;
   InternalPointer = PointsOnBoundary.begin();
+  LastTriangle = NULL;
+  TriangleFilesWritten = 0;
+}
+;
 …
       cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+  }
+  cout << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
+}
+;
 …
   Vector *Center = cloud->GetCenter(out);
   list<BoundaryTriangleSet*> *triangles = NULL;
+  bool AddFlag = false;
+  LinkedCell *BoundaryPoints = NULL;
   *out << Verbose(1) << "Begin of InsertStraddlingPoints" << endl;
   cloud->GoToFirst();
+  BoundaryPoints = new LinkedCell(this, 5.);
   while (!cloud->IsEnd()) {  // we only have to go once through all points, as boundary can become only bigger
+    LinkedCell BoundaryPoints(this, 5.);
+    if (AddFlag) {
+      delete(BoundaryPoints);
+      BoundaryPoints = new LinkedCell(this, 5.);
+      AddFlag = false;
+    }
     Walker = cloud->GetPoint();
     *out << Verbose(2) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
+    triangles = FindClosestTrianglesToPoint(out, Walker->node, &BoundaryPoints);
+    if (triangles == NULL) {
+      *out << Verbose(1) << "No triangles found, probably a tesselation point itself." << endl;
+    triangles = FindClosestTrianglesToPoint(out, Walker->node, BoundaryPoints);
+    BTS = triangles->front();
+    if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
+      *out << Verbose(2) << "No triangles found, probably a tesselation point itself." << endl;
       cloud->GoToNext();
       continue;
     } else {
-      BTS = triangles->front();
+    }
     *out << Verbose(2) << "Closest triangle is " << *BTS << "." << endl;
 …
         // add Walker to boundary points
         *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
+        AddFlag = true;
         if (AddBoundaryPoint(Walker,0))
           NewPoint = BPS[0];
 …
   } else {
     delete TPS[n];
     cout << Verbose(3) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
+    cout << Verbose(4) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
     TPS[n] = (InsertUnique.first)->second;
+  }
 …
   if (a->lines.find(b->node->nr) != a->lines.end()) {
     LineMap::iterator FindLine;
+    LineMap::iterator FindLine = a->lines.find(b->node->nr);
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
+    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+    cout << Verbose(5) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
+    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->triangles.size() < 2) {
         insertNewLine = false;
         cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
+        cout << Verbose(4) << "Using existing line " << *FindLine->second << endl;
         BPS[0] = FindLine->second->endpoints[0];
 …
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
+{
   cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
+  cout << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl;
   BPS[0] = a;
   BPS[1] = b;
 …
 };
 /** 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.
+/** Function adds triangle to global list.
+ * Furthermore, the triangle receives the next free id and id counter \a TrianglesOnBoundaryCount is increased.
  */
 void Tesselation::AddTesselationTriangle()
 …
   TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
   TrianglesOnBoundaryCount++;
+  // set as last new triangle
+  LastTriangle = BTS;
+  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
+};
+/** Function adds triangle to global list.
+ * Furthermore, the triangle number is set to \a nr.
+ * \param nr triangle number
+ */
+void Tesselation::AddTesselationTriangle(int nr)
+{
+  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
+  // add triangle to global map
+  TrianglesOnBoundary.insert(TrianglePair(nr, BTS));
+  // set as last new triangle
+  LastTriangle = BTS;
   // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
 …
           cout << Verbose(5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;
           RemoveTesselationLine(triangle->lines[i]);
+          triangle->lines[i] = NULL;
+      } else
+        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle." << endl;
+      } else {
+        cout << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle: ";
+        for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
+          cout << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t";
+        cout << endl;
+//        for (int j=0;j<2;j++) {
+//          cout << Verbose(5) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
+//          for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
+//            cout << "[" << *(LineRunner->second) << "] \t";
+//          cout << endl;
+//        }
+      }
+      triangle->lines[i] = NULL;  // free'd or not: disconnect
     } else
       cerr << "ERROR: This line " << i << " has already been free'd." << endl;
 …
   if (line == NULL)
     return;
   // get other endpoint number of finding copies of same line
+  // get other endpoint number for finding copies of same line
   if (line->endpoints[1] != NULL)
     Numbers[0] = line->endpoints[1]->Nr;
 …
         cout << Verbose(5) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;
         RemoveTesselationPoint(line->endpoints[i]);
+        line->endpoints[i] = NULL;
+      } else
+        cout << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to." << endl;
+      } else {
+        cout << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to: ";
+        for(LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
+          cout << "[" << *(LineRunner->second) << "] \t";
+        cout << endl;
+      }
+      line->endpoints[i] = NULL;  // free'd or not: disconnect
     } else
       cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
 …
+          }
           // 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;
+        }
+      }
 …
   *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
   return adjacentTriangleCount;
+};
+/** Checks whether the triangle consisting of the three points is already present.
+ * Searches for the points in Tesselation::PointsOnBoundary and checks their
+ * lines. If any of the three edges already has two triangles attached, false is
+ * returned.
+ * \param *out output stream for debugging
+ * \param *Candidates endpoints of the triangle candidate
+ * \return NULL - none found or pointer to triangle
+ */
+class BoundaryTriangleSet * Tesselation::GetPresentTriangle(ofstream *out, TesselPoint *Candidates[3])
+{
+  class BoundaryTriangleSet *triangle = NULL;
+  class BoundaryPointSet *Points[3];
+  // builds a triangle point set (Points) of the end points
+  for (int i = 0; i < 3; i++) {
+    PointMap::iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      Points[i] = FindPoint->second;
+    } else {
+      Points[i] = NULL;
+    }
+  }
+  // checks lines between the points in the Points for their adjacent triangles
+  for (int i = 0; i < 3; i++) {
+    if (Points[i] != NULL) {
+      for (int j = i; j < 3; j++) {
+        if (Points[j] != NULL) {
+          LineMap::iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
+          for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
+            TriangleMap *triangles = &FindLine->second->triangles;
+            for (TriangleMap::iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
+              if (FindTriangle->second->IsPresentTupel(Points)) {
+                if ((triangle == NULL) || (triangle->Nr > FindTriangle->second->Nr))
+                  triangle = FindTriangle->second;
+              }
+            }
+          }
+          // Only one of the triangle lines must be considered for the triangle count.
+          //*out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          //return adjacentTriangleCount;
+        }
+      }
+    }
+  }
+  return triangle;
 };
 …
   CandidateList *OptCandidates = new CandidateList();
   for (int k=0;k<NDIM;k++) {
+    Oben.Zero();
     Oben.x[k] = 1.;
     FirstPoint = MaxPoint[k];
 …
     ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
     FindSecondPointForTesselation(FirstPoint, NULL, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(FirstPoint, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     SecondPoint = OptCandidate;
     if (SecondPoint == NULL)  // have we found a second point?
       continue;
-    else
-      cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
     helper.CopyVector(FirstPoint->node);
 …
     // adding point 1 and point 2 and add the line between them
+    cout << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;
     AddTesselationPoint(FirstPoint, 0);
+    cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
     AddTesselationPoint(SecondPoint, 1);
     AddTesselationLine(TPS[0], TPS[1], 0);
 …
  * @param *LC LinkedCell structure with neighbouring points
  */
 bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC)
+bool Tesselation::FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC)
+{
   cout << Verbose(0) << "Begin of FindNextSuitableTriangle\n";
 …
     CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
     cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
     // construct old center
 …
         result = false;
+      }
     } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
+    } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
         AddTesselationPoint((*it)->point, 0);
         AddTesselationPoint(BaseRay->endpoints[0]->node, 1);
 …
     // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
     BaseRay = BLS[0];
+    if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) {
+      *out << Verbose(1) << "CRITICAL: Triangle " << *BTS << " has zero normal vector!" << endl;
+      exit(255);
+    }
+  }
 …
  * \param *out output stream for debugging
  * \param *Base line to be flipped
  * \return NULL - concave, otherwise endpoint that makes it concave
+ * \return NULL - convex, otherwise endpoint that makes it concave
  */
 class BoundaryPointSet *Tesselation::IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base)
 …
  * \param *out output stream for debugging
  * \param *Base line to be flipped
  * \return true - line was changed, false - same line as before
  */
 bool Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
+ * \return volume change due to flipping (0 - then no flipped occured)
+ */
+double Tesselation::PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base)
+{
   class BoundaryLineSet *OtherBase;
   Vector *ClosestPoint[2];
+  double volume;
   int m=0;
 …
   Distance.CopyVector(ClosestPoint[1]);
   Distance.SubtractVector(ClosestPoint[0]);
+  // calculate volume
+  volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->node, OtherBase->endpoints[0]->node->node, OtherBase->endpoints[1]->node->node, Base->endpoints[0]->node->node);
   // delete the temporary other base line and the closest points
 …
     if (Base->triangles.size() < 2) {
       *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
       return false;
+      return 0.;
+    }
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
 …
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
       *out << Verbose(2) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;
       FlipBaseline(out, Base);
       return true;
+      // calculate volume summand as a general tetraeder
+      return volume;
     } else {  // Base higher than OtherBase -> do nothing
       *out << Verbose(2) << "REJECT: Base line is higher: Nothing to do." << endl;
+      return false;
+    }
+  }
+};
+/** Returns the closest point on \a *Base with respect to \a *OtherBase.
+ * \param *out output stream for debugging
+ * \param *Base reference line
+ * \param *OtherBase other base line
+ * \return Vector on reference line that has closest distance
+ */
+Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
+{
+  // construct the plane of the two baselines (i.e. take both their directional vectors)
+  Vector Normal;
+  Vector Baseline, OtherBaseline;
+  Baseline.CopyVector(Base->endpoints[1]->node->node);
+  Baseline.SubtractVector(Base->endpoints[0]->node->node);
+  OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
+  OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
+  Normal.CopyVector(&Baseline);
+  Normal.VectorProduct(&OtherBaseline);
+  Normal.Normalize();
+  *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
+  // project one offset point of OtherBase onto this plane (and add plane offset vector)
+  Vector NewOffset;
+  NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
+  NewOffset.SubtractVector(Base->endpoints[0]->node->node);
+  NewOffset.ProjectOntoPlane(&Normal);
+  NewOffset.AddVector(Base->endpoints[0]->node->node);
+  Vector NewDirection;
+  NewDirection.CopyVector(&NewOffset);
+  NewDirection.AddVector(&OtherBaseline);
+  // calculate the intersection between this projected baseline and Base
+  Vector *Intersection = new Vector;
+  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
+  Normal.CopyVector(Intersection);
+  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
+  return Intersection;
+      return 0.;
+    }
+  }
 };
 …
  * \param *out output stream for debugging
  * \param *Base line to be flipped
  * \return true - flipping successful, false - something went awry
  */
 bool Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
+ * \return pointer to allocated new baseline - flipping successful, NULL - something went awry
+ */
+class BoundaryLineSet * Tesselation::FlipBaseline(ofstream *out, class BoundaryLineSet *Base)
+{
   class BoundaryLineSet *OldLines[4], *NewLine;
 …
   if (Base->triangles.size() < 2) {
     *out << Verbose(2) << "ERROR: Less than two triangles are attached to this baseline!" << endl;
     return false;
+    return NULL;
+  }
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
 …
   if (i<4) {
     *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
     return false;
+    return NULL;
+  }
   for (int j=0;j<4;j++)
     if (OldLines[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
       return false;
+      return NULL;
+    }
   for (int j=0;j<2;j++)
     if (OldPoints[j] == NULL) {
       *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
       return false;
+      return NULL;
+    }
 …
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
     BTS->GetNormalVector(BaseLineNormal);
     TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[0], BTS));
+    AddTesselationTriangle(OldTriangleNrs[0]);
     *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
 …
     BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
     BTS->GetNormalVector(BaseLineNormal);
     TrianglesOnBoundary.insert(TrianglePair(OldTriangleNrs[1], BTS));
+    AddTesselationTriangle(OldTriangleNrs[1]);
     *out << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl;
   } else {
     *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
     return false;
+    return NULL;
+  }
   *out << Verbose(1) << "End of FlipBaseline" << endl;
   return true;
+  return NewLine;
 };
 …
 /** Finds the second point of starting triangle.
  * \param *a first node
- * \param *Candidate pointer to candidate node on return
  * \param Oben vector indicating the outside
  * \param OptCandidate reference to recommended candidate on return
 …
  * \param *LC LinkedCell structure with neighbouring points
  */
 void Tesselation::FindSecondPointForTesselation(TesselPoint* a, TesselPoint* Candidate, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
+void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, LinkedCell *LC)
+{
   cout << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl;
   Vector AngleCheck;
+  class TesselPoint* Candidate = NULL;
   double norm = -1., angle;
   LinkedNodes *List = NULL;
 …
   cout << Verbose(1) << "Begin of FindThirdPointForTesselation" << endl;
   //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+  cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
   // construct center of circle
 …
     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
+      //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
       // check SearchDirection
 …
   TesselPoint *trianglePoints[3];
   TesselPoint *SecondPoint = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
   if (LinesOnBoundary.empty()) {
 …
   // check whether closest point is "too close" :), then it's inside
   if (trianglePoints[0] == NULL) {
     *out << Verbose(1) << "Is the only point, no one else is closeby." << endl;
+    *out << Verbose(2) << "Is the only point, no one else is closeby." << endl;
     return NULL;
+  }
   if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
+    *out << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
+    return NULL;
+  }
+  list<TesselPoint*> *connectedPoints = GetCircleOfConnectedPoints(out, trianglePoints[0]);
+  list<TesselPoint*> *connectedClosestPoints = GetNeighboursOnCircleOfConnectedPoints(out, connectedPoints, trianglePoints[0], x);
+  delete(connectedPoints);
+  trianglePoints[1] = connectedClosestPoints->front();
+  trianglePoints[2] = connectedClosestPoints->back();
+  for (int i=0;i<3;i++) {
+    if (trianglePoints[i] == NULL) {
+      *out << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+    }
+    //*out << Verbose(1) << "List of possible points:" << endl;
+    //*out << Verbose(2) << *trianglePoints[i] << endl;
+  }
+  list<BoundaryTriangleSet*> *triangles = FindTriangles(trianglePoints);
+  delete(connectedClosestPoints);
+    *out << Verbose(3) << "Point is right on a tesselation point, no nearest triangle." << endl;
+    PointMap::iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
+    triangles = new list<BoundaryTriangleSet*>;
+    if (PointRunner != PointsOnBoundary.end()) {
+      for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
+        for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
+          triangles->push_back(TriangleRunner->second);
+      triangles->sort();
+      triangles->unique();
+    } else {
+      PointRunner = PointsOnBoundary.find(SecondPoint->nr);
+      trianglePoints[0] = SecondPoint;
+      if (PointRunner != PointsOnBoundary.end()) {
+        for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
+          for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
+            triangles->push_back(TriangleRunner->second);
+        triangles->sort();
+        triangles->unique();
+      } else {
+        *out << Verbose(1) << "ERROR: I cannot find a boundary point to the tessel point " << *trianglePoints[0] << "." << endl;
+        return NULL;
+      }
+    }
+  } else {
+    list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(out, trianglePoints[0], x);
+    trianglePoints[1] = connectedClosestPoints->front();
+    trianglePoints[2] = connectedClosestPoints->back();
+    for (int i=0;i<3;i++) {
+      if (trianglePoints[i] == NULL) {
+        *out << Verbose(1) << "ERROR: IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+      }
+      //*out << Verbose(2) << "List of triangle points:" << endl;
+      //*out << Verbose(3) << *trianglePoints[i] << endl;
+    }
+    triangles = FindTriangles(trianglePoints);
+    *out << Verbose(2) << "List of possible triangles:" << endl;
+    for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      *out << Verbose(3) << **Runner << endl;
+    delete(connectedClosestPoints);
+  }
   if (triangles->empty()) {
+    *out << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
+    *out << Verbose(0) << "ERROR: There is no nearest triangle. Please check the tesselation structure.";
+    delete(triangles);
     return NULL;
   } else
 …
   class BoundaryTriangleSet *result = NULL;
   list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(out, x, LC);
+  Vector Center;
   if (triangles == NULL)
     return NULL;
+  if (x->ScalarProduct(&triangles->front()->NormalVector) < 0)
+    result = triangles->back();
+  else
+  if (triangles->size() == 1) { // there is no degenerate case
     result = triangles->front();
+    *out << Verbose(2) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
+  } else {
+    result = triangles->front();
+    result->GetCenter(&Center);
+    Center.SubtractVector(x);
+    *out << Verbose(2) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
+    if (Center.ScalarProduct(&result->NormalVector) < 0) {
+      result = triangles->back();
+      *out << Verbose(2) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
+      if (Center.ScalarProduct(&result->NormalVector) < 0) {
+        *out << Verbose(1) << "ERROR: Front and back side yield NormalVector in wrong direction!" << endl;
+      }
+    }
+  }
   delete(triangles);
   return result;
 …
+{
   class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, &Point, LC);
+  if (result == NULL)
+  Vector Center;
+  if (result == NULL) {// is boundary point or only point in point cloud?
+    *out << Verbose(1) << Point << " is the only point in vicinity." << endl;
+    return false;
+  }
+  result->GetCenter(&Center);
+  *out << Verbose(3) << "INFO: Central point of the triangle is " << Center << "." << endl;
+  Center.SubtractVector(&Point);
+  *out << Verbose(3) << "INFO: Vector from center to point to test is " << Center << "." << endl;
+  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
+    *out << Verbose(1) << Point << " is an inner point." << endl;
     return true;
+  if (Point.ScalarProduct(&result->NormalVector) < 0)
+    return true;
+  else
+  } else {
+    *out << Verbose(1) << Point << " is NOT an inner point." << endl;
     return false;
+  }
+}
 …
 bool Tesselation::IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC)
+{
+  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(out, Point->node, LC);
+  if (result == NULL)
+    return true;
+  if (Point->node->ScalarProduct(&result->NormalVector) < 0)
+    return true;
+  else
+    return false;
+  return IsInnerPoint(out, *(Point->node), LC);
+}
 …
  * @param *Point of which get all connected points
+ *
  * @return list of the all points linked to the provided one
  */
 list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
   list<TesselPoint*> *connectedPoints = new list<TesselPoint*>;
+ * @return set of the all points linked to the provided one
+ */
+set<TesselPoint*> * Tesselation::GetAllConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
 …
     ReferencePoint = PointRunner->second;
   } else {
     *out << Verbose(2) << "getCircleOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    *out << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
     ReferencePoint = NULL;
+  }
   // little trick so that we look just through lines connect to the BoundaryPoint
+  // little trick so that we look just through lines connect to the BoundaryPoint
   // OR fall-back to look through all lines if there is no such BoundaryPoint
   LineMap *Lines = &LinesOnBoundary;
 …
   LineMap::iterator findLines = Lines->begin();
   while (findLines != Lines->end()) {
     takePoint = false;
     if (findLines->second->endpoints[0]->Nr == Point->nr) {
       takePoint = true;
       current = findLines->second->endpoints[1]->node;
     } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
       takePoint = true;
       current = findLines->second->endpoints[0]->node;
+    }
     if (takePoint) {
       *out << Verbose(3) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is taken into the circle." << endl;
       connectedPoints->push_back(current);
+    }
     findLines++;
+   takePoint = false;
+   if (findLines->second->endpoints[0]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[1]->node;
+   } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
+     takePoint = true;
+     current = findLines->second->endpoints[0]->node;
+   }
+   if (takePoint) {
+     *out << Verbose(5) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;
+     connectedPoints->insert(current);
+   }
+   findLines++;
+  }
 …
     return NULL;
+  }
   return connectedPoints;
+}
+/** Gets the two neighbouring points with respect to a reference line to the provided point.
+};
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
  * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
  * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
 …
+ *
  * @param *out output stream for debugging
- * @param *connectedPoints list of connected points to the central \a *Point
  * @param *Point of which get all connected points
+ * @param *Reference Vector to be checked whether it is an inner point
+ *
+ * @return list of the two points linked to the provided one and closest to the point to be checked,
+ */
+list<TesselPoint*> * Tesselation::GetNeighboursOnCircleOfConnectedPoints(ofstream *out, list<TesselPoint*> *connectedPoints, TesselPoint* Point, Vector* Reference)
+ * @param *Reference Reference vector for zero angle or NULL for no preference
+ * @return list of the all points linked to the provided one
+ */
+list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector *Reference)
+{
   map<double, TesselPoint*> anglesOfPoints;
+  map<double, TesselPoint*>::iterator runner;
+  ;
+  Vector center, PlaneNormal, OrthogonalVector, helper, AngleZero;
+  if (connectedPoints->size() == 0) { // if have not found any points
+    *out << Verbose(1) << "ERROR: We have not found any connected points to " << *Point<< "." << endl;
+    return NULL;
+  }
+  set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(out, Point);
+  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
   // calculate central point
   for (list<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
+  for (set<TesselPoint*>::iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //*out << "Summed vectors " << center << "; number of points " << connectedPoints.size()
 …
   // construct one orthogonal vector
+  AngleZero.CopyVector(Reference);
+  if (Reference != NULL)
+    AngleZero.CopyVector(Reference);
+  else
+    AngleZero.CopyVector((*connectedPoints->begin())->node);
   AngleZero.SubtractVector(Point->node);
   AngleZero.ProjectOntoPlane(&PlaneNormal);
 …
   // go through all connected points and calculate angle
   for (list<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
+  for (set<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
     helper.ProjectOntoPlane(&PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     *out << Verbose(2) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
+    *out << Verbose(3) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
     anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+  }
+  list<TesselPoint*> *result = new list<TesselPoint*>;
+  runner = anglesOfPoints.begin();
+  result->push_back(runner->second);
+  runner = anglesOfPoints.end();
+  runner--;
+  result->push_back(runner->second);
+  *out << Verbose(2) << "List of closest points has " << result->size() << " elements, which are "
+    << *(result->front()) << " and " << *(result->back()) << endl;
+  return result;
+  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
+    connectedCircle->push_back(AngleRunner->second);
+  }
+  delete(connectedPoints);
+  return connectedCircle;
+}
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we walk along a closed path.
+ *
+ * @param *out output stream for debugging
+ * @param *Point of which get all connected points
+ * @return list of the all points linked to the provided one
+ */
+list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  map<double, TesselPoint*> anglesOfPoints;
+  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+  class BoundaryPointSet *ReferencePoint = NULL;
+  class BoundaryPointSet *CurrentPoint = NULL;
+  class BoundaryTriangleSet *triangle = NULL;
+  class BoundaryLineSet *CurrentLine = NULL;
+  class BoundaryLineSet *StartLine = NULL;
+  // find the respective boundary point
+  PointMap::iterator PointRunner = PointsOnBoundary.find(Point->nr);
+  if (PointRunner != PointsOnBoundary.end()) {
+    ReferencePoint = PointRunner->second;
+  } else {
+    *out << Verbose(2) << "ERROR: GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;
+    return NULL;
+  }
+  map <class BoundaryLineSet *, bool> TouchedLine;
+  map <class BoundaryTriangleSet *, bool> TouchedTriangle;
+  map <class BoundaryLineSet *, bool>::iterator LineRunner;
+  map <class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
+  for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
+    TouchedLine.insert( pair <class BoundaryLineSet *, bool>(Runner->second, false) );
+    for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++)
+      TouchedTriangle.insert( pair <class BoundaryTriangleSet *, bool>(Sprinter->second, false) );
+  }
+  if (!ReferencePoint->lines.empty()) {
+    for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++) {
+      LineRunner = TouchedLine.find(runner->second);
+      if (LineRunner == TouchedLine.end()) {
+        *out << Verbose(2) << "ERROR: I could not find " << *runner->second << " in the touched list." << endl;
+      } else if (!LineRunner->second) {
+        LineRunner->second = true;
+        connectedPath = new list<TesselPoint*>;
+        triangle = NULL;
+        CurrentLine = runner->second;
+        StartLine = CurrentLine;
+        CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+        *out << Verbose(3)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;
+        do {
+          // push current one
+          *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
+          connectedPath->push_back(CurrentPoint->node);
+          // find next triangle
+          for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
+            *out << Verbose(3) << "INFO: Inspecting triangle " << *Runner->second << "." << endl;
+            if ((Runner->second != triangle)) { // look for first triangle not equal to old one
+              triangle = Runner->second;
+              TriangleRunner = TouchedTriangle.find(triangle);
+              if (TriangleRunner != TouchedTriangle.end()) {
+                if (!TriangleRunner->second) {
+                  TriangleRunner->second = true;
+                  *out << Verbose(3) << "INFO: Connecting triangle is " << *triangle << "." << endl;
+                  break;
+                } else {
+                  *out << Verbose(3) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl;
+                  triangle = NULL;
+                }
+              } else {
+                *out << Verbose(2) << "ERROR: I could not find " << *triangle << " in the touched list." << endl;
+                triangle = NULL;
+              }
+            }
+          }
+          if (triangle == NULL)
+            break;
+          // find next line
+          for (int i=0;i<3;i++) {
+            if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
+              CurrentLine = triangle->lines[i];
+              *out << Verbose(3) << "INFO: Connecting line is " << *CurrentLine << "." << endl;
+              break;
+            }
+          }
+          LineRunner = TouchedLine.find(CurrentLine);
+          if (LineRunner == TouchedLine.end())
+            *out << Verbose(2) << "ERROR: I could not find " << *CurrentLine << " in the touched list." << endl;
+          else
+            LineRunner->second = true;
+          // find next point
+          CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
+        } while (CurrentLine != StartLine);
+        // last point is missing, as it's on start line
+        *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
+        if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
+          connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+        ListOfPaths->push_back(connectedPath);
+      } else {
+        *out << Verbose(3) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl;
+      }
+    }
+  } else {
+    *out << Verbose(1) << "ERROR: There are no lines attached to " << *ReferencePoint << "." << endl;
+  }
+  return ListOfPaths;
+}
+/** Gets all closed paths on the circle of points connected to the provided point by triangulation lines, if this very point is removed.
+ * From GetPathsOfConnectedPoints() extracts all single loops of intracrossing paths in the list of closed paths.
+ * @param *out output stream for debugging
+ * @param *Point of which get all connected points
+ * @return list of the closed paths
+ */
+list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(ofstream *out, TesselPoint* Point)
+{
+  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(out, Point);
+  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
+  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPoint*> *newPath = NULL;
+  int count = 0;
+  list<TesselPoint*>::iterator CircleRunner;
+  list<TesselPoint*>::iterator CircleStart;
+  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+    connectedPath = *ListRunner;
+    *out << Verbose(2) << "INFO: Current path is " << connectedPath << "." << endl;
+    // go through list, look for reappearance of starting Point and count
+    CircleStart = connectedPath->begin();
+    // go through list, look for reappearance of starting Point and create list
+    list<TesselPoint*>::iterator Marker = CircleStart;
+    for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
+      if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
+        // we have a closed circle from Marker to new Marker
+        *out << Verbose(3) << count+1 << ". closed path consists of: ";
+        newPath = new list<TesselPoint*>;
+        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        for (; CircleSprinter != CircleRunner; CircleSprinter++) {
+          newPath->push_back(*CircleSprinter);
+          *out << (**CircleSprinter) << " <-> ";
+        }
+        *out << ".." << endl;
+        count++;
+        Marker = CircleRunner;
+        // add to list
+        ListofClosedPaths->push_back(newPath);
+      }
+    }
+  }
+  *out << Verbose(3) << "INFO: " << count << " closed additional path(s) have been created." << endl;
+  // delete list of paths
+  while (!ListofPaths->empty()) {
+    connectedPath = *(ListofPaths->begin());
+    ListofPaths->remove(connectedPath);
+    delete(connectedPath);
+  }
+  delete(ListofPaths);
+  // exit
+  return ListofClosedPaths;
+};
+/** Gets all belonging triangles for a given BoundaryPointSet.
+ * \param *out output stream for debugging
+ * \param *Point BoundaryPoint
+ * \return pointer to allocated list of triangles
+ */
+set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(ofstream *out, class BoundaryPointSet *Point)
+{
+  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+  if (Point == NULL) {
+    *out << Verbose(1) << "ERROR: Point given is NULL." << endl;
+  } else {
+    // go through its lines and insert all triangles
+    for (LineMap::iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
+      for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
+      connectedTriangles->insert(TriangleRunner->second);
+    }
+  }
+  return connectedTriangles;
+};
 /** Removes a boundary point from the envelope while keeping it closed.
+ * We create new triangles and remove the old ones connected to the point.
+ * We remove the old triangles connected to the point and re-create new triangles to close the surface following this ansatz:
+ *  -# a closed path(s) of boundary points surrounding the point to be removed is constructed
+ *  -# on each closed path, we pick three adjacent points, create a triangle with them and subtract the middle point from the path
+ *  -# we advance two points (i.e. the next triangle will start at the ending point of the last triangle) and continue as before
+ *  -# the surface is closed, when the path is empty
+ * Thereby, we (hopefully) make sure that the removed points remains beneath the surface (this is checked via IsInnerPoint eventually).
  * \param *out output stream for debugging
  * \param *point point to be removed
 …
   class BoundaryLineSet *line = NULL;
   class BoundaryTriangleSet *triangle = NULL;
   Vector OldPoint, TetraederVector[3];
+  Vector OldPoint, NormalVector;
   double volume = 0;
-  int *numbers = NULL;
   int count = 0;
-  int i;
   if (point == NULL) {
 …
     return 0.;
+  }
+  list<TesselPoint*> *CircleofPoints = GetCircleOfConnectedPoints(out, point->node);
+  // remove all triangles
+  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(out, point->node);
+  list<TesselPoint*> *connectedPath = NULL;
+  // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
+  numbers = new int[count];
+  class BoundaryTriangleSet **Candidates = new BoundaryTriangleSet*[count];
+  i=0;
+  for (LineMap::iterator LineRunner = point->lines.begin(); (point != NULL) && (LineRunner != point->lines.end()); LineRunner++) {
+  map<class BoundaryTriangleSet *, int> Candidates;
+  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
+      Candidates[i] = triangle;
+      numbers[i++] = triangle->Nr;
+    }
+  }
+  for (int j=0;j<i;j++) {
+    RemoveTesselationTriangle(Candidates[j]);
+  }
+  delete[](Candidates);
+  *out << Verbose(1) << i << " triangles were removed." << endl;
+  // re-create all triangles by going through connected points list
+  list<TesselPoint*>::iterator CircleRunner = CircleofPoints->begin();
+  list<TesselPoint*>::iterator OtherCircleRunner = CircleofPoints->begin();
+  class TesselPoint *CentralNode = *CircleRunner;
+  // advance two with CircleRunner and one with OtherCircleRunner
+  CircleRunner++;
+  CircleRunner++;
+  OtherCircleRunner++;
+  i=0;
+  cout << Verbose(2) << "INFO: CentralNode is " << *CentralNode << "." << endl;
+  for (; (OtherCircleRunner != CircleofPoints->end()) && (CircleRunner != CircleofPoints->end()); (CircleRunner++), (OtherCircleRunner++)) {
+    cout << Verbose(3) << "INFO: CircleRunner's node is " << **CircleRunner << "." << endl;
+    cout << Verbose(3) << "INFO: OtherCircleRunner's node is " << **OtherCircleRunner << "." << endl;
+    *out << Verbose(4) << "Adding new triangle points."<< endl;
+    AddTesselationPoint(CentralNode, 0);
+    AddTesselationPoint(*OtherCircleRunner, 1);
+    AddTesselationPoint(*CircleRunner, 2);
+    *out << Verbose(4) << "Adding new triangle lines."<< endl;
+    AddTesselationLine(TPS[0], TPS[1], 0);
+    AddTesselationLine(TPS[0], TPS[2], 1);
+    AddTesselationLine(TPS[1], TPS[2], 2);
+    BTS = new class BoundaryTriangleSet(BLS, numbers[i]);
+    TrianglesOnBoundary.insert(TrianglePair(numbers[i], BTS));
+    *out << Verbose(4) << "Created triangle " << *BTS << "." << endl;
+    // calculate volume summand as a general tetraeder
+    for (int j=0;j<3;j++) {
+      TetraederVector[j].CopyVector(TPS[j]->node->node);
+      TetraederVector[j].SubtractVector(&OldPoint);
+    }
+    OldPoint.CopyVector(&TetraederVector[0]);
+    OldPoint.VectorProduct(&TetraederVector[1]);
+    volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
+    // advance number
+    i++;
+    if (i >= count)
+      *out << Verbose(2) << "WARNING: Maximum of numbers reached!" << endl;
+  }
+  *out << Verbose(1) << i << " triangles were created." << endl;
+  delete[](numbers);
+  return volume;
+};
+/** 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->triangles.size() < 2) {
+            counter++;
+            break;  // increase counter only once per edge
+          }
+      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+    }
+  }
+  // remove all triangles
+  count=0;
+  NormalVector.Zero();
+  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    *out << Verbose(3) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
+    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->first);
+    count++;
+  }
+  *out << Verbose(1) << count << " triangles were removed." << endl;
+  list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
+  map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
+  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  double angle;
+  double smallestangle;
+  Vector Point, Reference, OrthogonalVector;
+  if (count > 2) {  // less than three triangles, then nothing will be created
+    class TesselPoint *TriangleCandidates[3];
+    count = 0;
+    for ( ; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) {  // go through all closed paths
+      if (ListAdvance != ListOfClosedPaths->end())
+        ListAdvance++;
+      connectedPath = *ListRunner;
+      // re-create all triangles by going through connected points list
+      list<class BoundaryLineSet *> NewLines;
+      for (;!connectedPath->empty();) {
+        // search middle node with widest angle to next neighbours
+        EndNode = connectedPath->end();
+        smallestangle = 0.;
+        for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
+          cout << Verbose(3) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+          // construct vectors to next and previous neighbour
+          StartNode = MiddleNode;
+          if (StartNode == connectedPath->begin())
+            StartNode = connectedPath->end();
+          StartNode--;
+          //cout << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
+          Point.CopyVector((*StartNode)->node);
+          Point.SubtractVector((*MiddleNode)->node);
+          StartNode = MiddleNode;
+          StartNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          //cout << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
+          Reference.CopyVector((*StartNode)->node);
+          Reference.SubtractVector((*MiddleNode)->node);
+          OrthogonalVector.CopyVector((*MiddleNode)->node);
+          OrthogonalVector.SubtractVector(&OldPoint);
+          OrthogonalVector.MakeNormalVector(&Reference);
+          angle = GetAngle(Point, Reference, OrthogonalVector);
+          //if (angle < M_PI)  // no wrong-sided triangles, please?
+            if(fabs(angle - M_PI) < fabs(smallestangle - M_PI)) {  // get straightest angle (i.e. construct those triangles with smallest area first)
+              smallestangle = angle;
+              EndNode = MiddleNode;
+            }
+        }
+      } else { // no line
+        cout << Verbose(1) << "The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
+        result = true;
+        MiddleNode = EndNode;
+        if (MiddleNode == connectedPath->end()) {
+          cout << Verbose(1) << "CRITICAL: Could not find a smallest angle!" << endl;
+          exit(255);
+        }
+        StartNode = MiddleNode;
+        if (StartNode == connectedPath->begin())
+          StartNode = connectedPath->end();
+        StartNode--;
+        EndNode++;
+        if (EndNode == connectedPath->end())
+          EndNode = connectedPath->begin();
+        cout << Verbose(4) << "INFO: StartNode is " << **StartNode << "." << endl;
+        cout << Verbose(4) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;
+        cout << Verbose(4) << "INFO: EndNode is " << **EndNode << "." << endl;
+        *out << Verbose(3) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl;
+        TriangleCandidates[0] = *StartNode;
+        TriangleCandidates[1] = *MiddleNode;
+        TriangleCandidates[2] = *EndNode;
+        triangle = GetPresentTriangle(out, TriangleCandidates);
+        if (triangle != NULL) {
+          cout << Verbose(1) << "WARNING: New triangle already present, skipping!" << endl;
+          StartNode++;
+          MiddleNode++;
+          EndNode++;
+          if (StartNode == connectedPath->end())
+            StartNode = connectedPath->begin();
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+          continue;
+        }
+        *out << Verbose(5) << "Adding new triangle points."<< endl;
+        AddTesselationPoint(*StartNode, 0);
+        AddTesselationPoint(*MiddleNode, 1);
+        AddTesselationPoint(*EndNode, 2);
+        *out << Verbose(5) << "Adding new triangle lines."<< endl;
+        AddTesselationLine(TPS[0], TPS[1], 0);
+        AddTesselationLine(TPS[0], TPS[2], 1);
+        NewLines.push_back(BLS[1]);
+        AddTesselationLine(TPS[1], TPS[2], 2);
+        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+        BTS->GetNormalVector(NormalVector);
+        AddTesselationTriangle();
+        // calculate volume summand as a general tetraeder
+        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->node, TPS[1]->node->node, TPS[2]->node->node, &OldPoint);
+        // advance number
+        count++;
+        // prepare nodes for next triangle
+        StartNode = EndNode;
+        cout << Verbose(4) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl;
+        connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
+        if (connectedPath->size() == 2) { // we are done
+          connectedPath->remove(*StartNode); // remove the start node
+          connectedPath->remove(*EndNode); // remove the end node
+          break;
+        } else if (connectedPath->size() < 2) { // something's gone wrong!
+          cout << Verbose(1) << "CRITICAL: There are only two endpoints left!" << endl;
+          exit(255);
+        } else {
+          MiddleNode = StartNode;
+          MiddleNode++;
+          if (MiddleNode == connectedPath->end())
+            MiddleNode = connectedPath->begin();
+          EndNode = MiddleNode;
+          EndNode++;
+          if (EndNode == connectedPath->end())
+            EndNode = connectedPath->begin();
+        }
+      }
+    }
+  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;
+};
+/** Sort function for the candidate list.
+ */
+bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
+{
+  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->node);
+  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  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->node);
+  helper.SubtractVector(candidate1->point->node);
+  OrthogonalVector.CopyVector(&helper);
+  helper.VectorProduct(&BaseLineVector);
+  OrthogonalVector.SubtractVector(&helper);
+  OrthogonalVector.Normalize();
+  // calculate both angles and correct with in-plane vector
+  helper.CopyVector(candidate1->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  double phi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+  helper.CopyVector(candidate2->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  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;
+};
+/**
+ * Finds the point which is second closest to the provided one.
+ *
+ * @param Point to which to find the second closest other point
+ * @param linked cell structure
+ *
+ * @return point which is second closest to the provided one
+ */
+TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
+{
+  LinkedNodes *List = NULL;
+  TesselPoint* closestPoint = NULL;
+  TesselPoint* secondClosestPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  LC->SetIndexToVector(Point); // 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 << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper. Norm();
+            if (currentNorm < distance) {
+              // remember second point
+              secondDistance = distance;
+              secondClosestPoint = closestPoint;
+              // mark down new closest point
+              distance = currentNorm;
+              closestPoint = (*Runner);
+              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
+      // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
+      if (NewLines.size() > 1) {
+        list<class BoundaryLineSet *>::iterator Candidate;
+        class BoundaryLineSet *OtherBase = NULL;
+        double tmp, maxgain;
+        do {
+          maxgain = 0;
+          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+            tmp = PickFarthestofTwoBaselines(out, *Runner);
+            if (maxgain < tmp) {
+              maxgain = tmp;
+              Candidate = Runner;
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+          if (maxgain != 0) {
+            volume += maxgain;
+            cout << Verbose(3) << "Flipping baseline with highest volume" << **Candidate << "." << endl;
+            OtherBase = FlipBaseline(out, *Candidate);
+            NewLines.erase(Candidate);
+            NewLines.push_back(OtherBase);
+          }
+        } while (maxgain != 0.);
+      }
+  return secondClosestPoint;
+};
+/**
+ * Finds the point which is closest to the provided one.
+ *
+ * @param Point to which to find the closest other point
+ * @param SecondPoint the second closest other point on return, NULL if none found
+ * @param linked cell structure
+ *
+ * @return point which is closest to the provided one, NULL if none found
+ */
+TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
+{
+  LinkedNodes *List = NULL;
+  TesselPoint* closestPoint = NULL;
+  SecondPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  LC->SetIndexToVector(Point); // 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 << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper. Norm();
+            if (currentNorm < distance) {
+              secondDistance = distance;
+              SecondPoint = closestPoint;
+              distance = currentNorm;
+              closestPoint = (*Runner);
+              cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
+            } else if (currentNorm < secondDistance) {
+              secondDistance = currentNorm;
+              SecondPoint = (*Runner);
+              cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  return closestPoint;
+};
+      ListOfClosedPaths->remove(connectedPath);
+      delete(connectedPath);
+    }
+    *out << Verbose(1) << count << " triangles were created." << endl;
+  } else {
+    while (!ListOfClosedPaths->empty()) {
+      ListRunner = ListOfClosedPaths->begin();
+      connectedPath = *ListRunner;
+      ListOfClosedPaths->remove(connectedPath);
+      delete(connectedPath);
+    }
+    *out << Verbose(1) << "No need to create any triangles." << endl;
+  }
+  delete(ListOfClosedPaths);
+  *out << Verbose(1) << "Removed volume is " << volume << "." << endl;
+  return volume;
+};
 /**
 …
               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])
+                  )
+                ) {
+                if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
                   result->push_back(FindTriangle->second);
+                }
 …
 /**
+ * Finds all degenerated lines within the tesselation structure.
+ *
+ * @return map of keys of degenerated line pairs, each line occurs twice
+ *         in the list, once as key and once as value
+ */
+map<int, int> * Tesselation::FindAllDegeneratedLines()
+{
+  map<int, class BoundaryLineSet *> AllLines;
+  map<int, int> * DegeneratedLines = new map<int, int>;
+  // sanity check
+  if (LinesOnBoundary.empty()) {
+    cout << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
+    return DegeneratedLines;
+  }
+  LineMap::iterator LineRunner1;
+  pair<LineMap::iterator, bool> tester;
+  for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
+    tester = AllLines.insert( pair<int,BoundaryLineSet *> (LineRunner1->second->endpoints[0]->Nr, LineRunner1->second) );
+    if ((!tester.second) && (tester.first->second->endpoints[1]->Nr == LineRunner1->second->endpoints[1]->Nr)) { // found degenerated line
+      DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, tester.first->second->Nr) );
+      DegeneratedLines->insert ( pair<int, int> (tester.first->second->Nr, LineRunner1->second->Nr) );
+    }
+  }
+  AllLines.clear();
+  cout << Verbose(1) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
+  map<int,int>::iterator it;
+  for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++)
+      cout << Verbose(2) << (*it).first << " => " << (*it).second << endl;
+  return DegeneratedLines;
+}
+/**
  * Finds all degenerated triangles within the tesselation structure.
+ *
 …
  *         in the list, once as key and once as value
  */
+map<int, int> Tesselation::FindAllDegeneratedTriangles()
+{
+  map<int, int> DegeneratedTriangles;
+  // sanity check
+  if (LinesOnBoundary.empty()) {
+    cout << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure.";
+    return DegeneratedTriangles;
+  }
+  LineMap::iterator LineRunner1, LineRunner2;
+  for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
+    for (LineRunner2 = LinesOnBoundary.begin(); LineRunner2 != LinesOnBoundary.end(); ++LineRunner2) {
+      if ((LineRunner1->second != LineRunner2->second)
+        && (LineRunner1->second->endpoints[0] == LineRunner2->second->endpoints[0])
+        && (LineRunner1->second->endpoints[1] == LineRunner2->second->endpoints[1])
+      ) {
+        TriangleMap::iterator TriangleRunner1 = LineRunner1->second->triangles.begin(),
+          TriangleRunner2 = LineRunner2->second->triangles.begin();
+        for (; TriangleRunner1 != LineRunner1->second->triangles.end(); ++TriangleRunner1) {
+          for (; TriangleRunner2 != LineRunner2->second->triangles.end(); ++TriangleRunner2) {
+            if ((TriangleRunner1->second != TriangleRunner2->second)
+              && (TriangleRunner1->second->endpoints[0] == TriangleRunner2->second->endpoints[0])
+              && (TriangleRunner1->second->endpoints[1] == TriangleRunner2->second->endpoints[1])
+              && (TriangleRunner1->second->endpoints[2] == TriangleRunner2->second->endpoints[2])
+            ) {
+              DegeneratedTriangles[TriangleRunner1->second->Nr] = TriangleRunner2->second->Nr;
+              DegeneratedTriangles[TriangleRunner2->second->Nr] = TriangleRunner1->second->Nr;
+            }
+          }
+map<int, int> * Tesselation::FindAllDegeneratedTriangles()
+{
+  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
+  map<int, int> * DegeneratedTriangles = new map<int, int>;
+  TriangleMap::iterator TriangleRunner1, TriangleRunner2;
+  LineMap::iterator Liner;
+  class BoundaryLineSet *line1 = NULL, *line2 = NULL;
+  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+    // run over both lines' triangles
+    Liner = LinesOnBoundary.find(LineRunner->first);
+    if (Liner != LinesOnBoundary.end())
+      line1 = Liner->second;
+    Liner = LinesOnBoundary.find(LineRunner->second);
+    if (Liner != LinesOnBoundary.end())
+      line2 = Liner->second;
+    for (TriangleRunner1 = line1->triangles.begin(); TriangleRunner1 != line1->triangles.end(); ++TriangleRunner1) {
+      for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) {
+        if ((TriangleRunner1->second != TriangleRunner2->second)
+          && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr) );
+          DegeneratedTriangles->insert( pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr) );
+        }
+      }
+    }
+  }
+  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles.size() << " triangles." << endl;
+  delete(DegeneratedLines);
+  cout << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
   map<int,int>::iterator it;
   for (it = DegeneratedTriangles.begin(); it != DegeneratedTriangles.end(); it++)
+  for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       cout << Verbose(2) << (*it).first << " => " << (*it).second << endl;
 …
 void Tesselation::RemoveDegeneratedTriangles()
+{
+  map<int, int> DegeneratedTriangles = FindAllDegeneratedTriangles();
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles.begin();
+    TriangleKeyRunner != DegeneratedTriangles.end(); ++TriangleKeyRunner
+  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  TriangleMap::iterator finder;
+  BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
+  int count  = 0;
+  cout << Verbose(1) << "Begin of RemoveDegeneratedTriangles" << endl;
+  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+    TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
+    BoundaryTriangleSet *triangle = TrianglesOnBoundary.find(TriangleKeyRunner->first)->second,
+      *partnerTriangle = TrianglesOnBoundary.find(TriangleKeyRunner->second)->second;
+    finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
+    if (finder != TrianglesOnBoundary.end())
+      triangle = finder->second;
+    else
+      break;
+    finder = TrianglesOnBoundary.find(TriangleKeyRunner->second);
+    if (finder != TrianglesOnBoundary.end())
+      partnerTriangle = finder->second;
+    else
+      break;
     bool trianglesShareLine = false;
 …
       && (triangle->endpoints[0]->LinesCount > 2)
     ) {
+      // check whether we have to fix lines
+      BoundaryTriangleSet *Othertriangle = NULL;
+      BoundaryTriangleSet *OtherpartnerTriangle = NULL;
+      TriangleMap::iterator TriangleRunner;
+      for (int i = 0; i < 3; ++i)
+        for (int j = 0; j < 3; ++j)
+          if (triangle->lines[i] != partnerTriangle->lines[j]) {
+            // get the other two triangles
+            for (TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); ++TriangleRunner)
+              if (TriangleRunner->second != triangle) {
+                Othertriangle = TriangleRunner->second;
+              }
+            for (TriangleRunner = partnerTriangle->lines[i]->triangles.begin(); TriangleRunner != partnerTriangle->lines[i]->triangles.end(); ++TriangleRunner)
+              if (TriangleRunner->second != partnerTriangle) {
+                OtherpartnerTriangle = TriangleRunner->second;
+              }
+            /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]
+            // the line of triangle receives the degenerated ones
+            triangle->lines[i]->triangles.erase(Othertriangle->Nr);
+            triangle->lines[i]->triangles.insert( TrianglePair( partnerTriangle->Nr, partnerTriangle) );
+            for (int k=0;k<3;k++)
+              if (triangle->lines[i] == Othertriangle->lines[k]) {
+                Othertriangle->lines[k] = partnerTriangle->lines[j];
+                break;
+              }
+            // the line of partnerTriangle receives the non-degenerated ones
+            partnerTriangle->lines[j]->triangles.erase( partnerTriangle->Nr);
+            partnerTriangle->lines[j]->triangles.insert( TrianglePair( Othertriangle->Nr, Othertriangle) );
+            partnerTriangle->lines[j] = triangle->lines[i];
+          }
+      // erase the pair
+      count += (int) DegeneratedTriangles->erase(triangle->Nr);
       cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;
       RemoveTesselationTriangle(triangle);
+      count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
       cout << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;
       RemoveTesselationTriangle(partnerTriangle);
-      DegeneratedTriangles.erase(DegeneratedTriangles.find(partnerTriangle->Nr));
     } else {
       cout << Verbose(1) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle
 …
+    }
+  }
+  delete(DegeneratedTriangles);
+  cout << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl;
+  cout << Verbose(1) << "End of RemoveDegeneratedTriangles" << endl;
+}
+/** Gets the angle between a point and a reference relative to the provided center.
+ * We have two shanks point and reference between which the angle is calculated
+ * and by scalar product with OrthogonalVector we decide the interval.
+ * @param point to calculate the angle for
+ * @param reference to which to calculate the angle
+ * @param OrthogonalVector points in direction of [pi,2pi] interval
+ *
+ * @return angle between point and reference
+ */
+double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
+{
+  if (reference.IsZero())
+    return M_PI;
+  // calculate both angles and correct with in-plane vector
+  if (point.IsZero())
+    return M_PI;
+  double phi = point.Angle(&reference);
+  if (OrthogonalVector.ScalarProduct(&point) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+  cout << Verbose(3) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
+  return phi;
+}
+/** Adds an outside Tesselpoint to the envelope via (two) degenerated triangles.
+ * We look for the closest point on the boundary, we look through its connected boundary lines and
+ * seek the one with the minimum angle between its center point and the new point and this base line.
+ * We open up the line by adding a degenerated triangle, whose other side closes the base line again.
+ * \param *out output stream for debugging
+ * \param *point point to add
+ * \param *LC Linked Cell structure to find nearest point
+ */
+void Tesselation::AddBoundaryPointByDegeneratedTriangle(ofstream *out, class TesselPoint *point, LinkedCell *LC)
+{
+  *out << Verbose(2) << "Begin of AddBoundaryPointByDegeneratedTriangle" << endl;
+  // find nearest boundary point
+  class TesselPoint *BackupPoint = NULL;
+  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
+  class BoundaryPointSet *NearestBoundaryPoint = NULL;
+  PointMap::iterator PointRunner;
+  if (NearestPoint == point)
+    NearestPoint = BackupPoint;
+  PointRunner = PointsOnBoundary.find(NearestPoint->nr);
+  if (PointRunner != PointsOnBoundary.end()) {
+    NearestBoundaryPoint = PointRunner->second;
+  } else {
+    *out << Verbose(1) << "ERROR: I cannot find the boundary point." << endl;
+    return;
+  }
+  *out << Verbose(2) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl;
+  // go through its lines and find the best one to split
+  Vector CenterToPoint;
+  Vector BaseLine;
+  double angle, BestAngle = 0.;
+  class BoundaryLineSet *BestLine = NULL;
+  for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
+    BaseLine.CopyVector(Runner->second->endpoints[0]->node->node);
+    BaseLine.SubtractVector(Runner->second->endpoints[1]->node->node);
+    CenterToPoint.CopyVector(Runner->second->endpoints[0]->node->node);
+    CenterToPoint.AddVector(Runner->second->endpoints[1]->node->node);
+    CenterToPoint.Scale(0.5);
+    CenterToPoint.SubtractVector(point->node);
+    angle = CenterToPoint.Angle(&BaseLine);
+    if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
+      BestAngle = angle;
+      BestLine = Runner->second;
+    }
+  }
+  // remove one triangle from the chosen line
+  class BoundaryTriangleSet *TempTriangle = (BestLine->triangles.begin())->second;
+  BestLine->triangles.erase(TempTriangle->Nr);
+  int nr = -1;
+  for (int i=0;i<3; i++) {
+    if (TempTriangle->lines[i] == BestLine) {
+      nr = i;
+      break;
+    }
+  }
+  // create new triangle to connect point (connects automatically with the missing spot of the chosen line)
+  *out << Verbose(5) << "Adding new triangle points."<< endl;
+  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
+  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
+  AddTesselationPoint(point, 2);
+  *out << Verbose(5) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  BTS->GetNormalVector(TempTriangle->NormalVector);
+  BTS->NormalVector.Scale(-1.);
+  *out << Verbose(3) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl;
+  AddTesselationTriangle();
+  // create other side of this triangle and close both new sides of the first created triangle
+  *out << Verbose(5) << "Adding new triangle points."<< endl;
+  AddTesselationPoint((BestLine->endpoints[0]->node), 0);
+  AddTesselationPoint((BestLine->endpoints[1]->node), 1);
+  AddTesselationPoint(point, 2);
+  *out << Verbose(5) << "Adding new triangle lines."<< endl;
+  AddTesselationLine(TPS[0], TPS[1], 0);
+  AddTesselationLine(TPS[0], TPS[2], 1);
+  AddTesselationLine(TPS[1], TPS[2], 2);
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  BTS->GetNormalVector(TempTriangle->NormalVector);
+  *out << Verbose(3) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;
+  AddTesselationTriangle();
+  // add removed triangle to the last open line of the second triangle
+  for (int i=0;i<3;i++) { // look for the same line as BestLine (only it's its degenerated companion)
+    if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
+      if (BestLine == BTS->lines[i]){
+        *out << Verbose(1) << "CRITICAL: BestLine is same as found line, something's wrong here!" << endl;
+        exit(255);
+      }
+      BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) );
+      TempTriangle->lines[nr] = BTS->lines[i];
+      break;
+    }
+  }
+  // exit
+  *out << Verbose(2) << "End of AddBoundaryPointByDegeneratedTriangle" << endl;
+};
+/** Writes the envelope to file.
+ * \param *out otuput stream for debugging
+ * \param *filename basename of output file
+ * \param *cloud PointCloud structure with all nodes
+ */
+void Tesselation::Output(ofstream *out, const char *filename, PointCloud *cloud)
+{
+  ofstream *tempstream = NULL;
+  string NameofTempFile;
+  char NumberName[255];
+  if (LastTriangle != NULL) {
+    sprintf(NumberName, "-%04d-%s_%s_%s", (int)TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->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);
+      *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      WriteTecplotFile(out, tempstream, this, cloud, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoRaster3DOutput) {
+      string NameofTempFile(filename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+      NameofTempFile.erase(npos, 1);
+      NameofTempFile.append(Raster3DSuffix);
+      *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      WriteRaster3dFile(out, tempstream, this, cloud);
+      IncludeSphereinRaster3D(out, tempstream, this, cloud);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+  }
+  if (DoTecplotOutput || DoRaster3DOutput)
+    TriangleFilesWritten++;
+};

src/tesselation.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 …
 #include <set>
-#include "helpers.hpp"
-#include "linkedcell.hpp"
-#include "tesselationhelpers.hpp"
 #include "vector.hpp"
+/****************************************** forward declarations *****************************/
 class BoundaryPointSet;
 class BoundaryLineSet;
 class BoundaryTriangleSet;
+class LinkedCell;
 class TesselPoint;
 class PointCloud;
 class Tesselation;
+/********************************************** definitions *********************************/
+#define DoTecplotOutput 1
+#define DoRaster3DOutput 1
+#define DoVRMLOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
 // ======================================================= some template functions =========================================
 …
 #define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
+/********************************************** declarations *******************************/
 template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2)
 …
     void GetNormalVector(Vector &NormalVector);
+    void GetCenter(Vector *center);
     bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
     bool ContainsBoundaryLine(class BoundaryLineSet *line);
     bool ContainsBoundaryPoint(class BoundaryPointSet *point);
+    bool ContainsBoundaryPoint(class TesselPoint *point);
     class BoundaryPointSet *GetThirdEndpoint(class BoundaryLineSet *line);
     bool IsPresentTupel(class BoundaryPointSet *Points[3]);
     void GetCenter(Vector *center);
+    bool IsPresentTupel(class BoundaryTriangleSet *T);
     class BoundaryPointSet *endpoints[3];
 …
     void AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
     void AddTesselationTriangle();
+    void AddTesselationTriangle(int nr);
     void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle);
     void RemoveTesselationLine(class BoundaryLineSet *line);
     void RemoveTesselationPoint(class BoundaryPointSet *point);
-    bool IsInside(Vector *pointer);
     class BoundaryPointSet *GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2);
     // concave envelope
     void FindStartingTriangle(ofstream *out, const double RADIUS, class LinkedCell *LC);
     void FindSecondPointForTesselation(class TesselPoint* a, class TesselPoint* Candidate, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, class LinkedCell *LC);
+    void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, class LinkedCell *LC);
     void FindThirdPointForTesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, class LinkedCell *LC);
     bool FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC);
+    bool FindNextSuitableTriangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, LinkedCell *LC);
     int CheckPresenceOfTriangle(ofstream *out, class TesselPoint *Candidates[3]);
+    class BoundaryTriangleSet * GetPresentTriangle(ofstream *out, TesselPoint *Candidates[3]);
     // convex envelope
 …
     bool InsertStraddlingPoints(ofstream *out, PointCloud *cloud, LinkedCell *LC);
     double RemovePointFromTesselatedSurface(ofstream *out, class BoundaryPointSet *point);
     bool FlipBaseline(ofstream *out, class BoundaryLineSet *Base);
     bool PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base);
+    class BoundaryLineSet * FlipBaseline(ofstream *out, class BoundaryLineSet *Base);
+    double PickFarthestofTwoBaselines(ofstream *out, class BoundaryLineSet *Base);
     class BoundaryPointSet *IsConvexRectangle(ofstream *out, class BoundaryLineSet *Base);
+    map<int, int> FindAllDegeneratedTriangles();
+    map<int, int> * FindAllDegeneratedTriangles();
+    map<int, int> * FindAllDegeneratedLines();
     void RemoveDegeneratedTriangles();
+    list<TesselPoint*> * GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point);
+    list<TesselPoint*> * GetNeighboursOnCircleOfConnectedPoints(ofstream *out, list<TesselPoint*> *connectedPoints, TesselPoint* Point, Vector* Reference);
+    void AddBoundaryPointByDegeneratedTriangle(ofstream *out, class TesselPoint *point, LinkedCell *LC);
+    set<TesselPoint*> * GetAllConnectedPoints(ofstream *out, TesselPoint* Point);
+    set<BoundaryTriangleSet*> *GetAllTriangles(ofstream *out, class BoundaryPointSet *Point);
+    list<list<TesselPoint*> *> * GetPathsOfConnectedPoints(ofstream *out, TesselPoint* Point);
+    list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(ofstream *out, TesselPoint* Point);
+    list<TesselPoint*> * GetCircleOfConnectedPoints(ofstream *out, TesselPoint* Point, Vector *Reference = NULL);
     list<BoundaryTriangleSet*> *FindTriangles(TesselPoint* Points[3]);
     list<BoundaryTriangleSet*> * FindClosestTrianglesToPoint(ofstream *out, Vector *x, LinkedCell* LC);
 …
     void PrintAllBoundaryTriangles(ofstream *out);
+    // store envelope in file
+    void Output(ofstream *out, const char *filename, PointCloud *cloud);
     PointMap PointsOnBoundary;
 …
     class BoundaryLineSet *BLS[3];
     class BoundaryTriangleSet *BTS;
+    class BoundaryTriangleSet *LastTriangle;
+    int TriangleFilesWritten;
   private:
 …
 };
-bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3]);
-bool SortCandidates(class CandidateForTesselation* candidate1, class CandidateForTesselation* candidate2);
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC);
-TesselPoint* FindSecondClosestPoint(const Vector*, LinkedCell*);
-double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector);
-Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase);
 #endif /* TESSELATION_HPP_ */

src/tesselationhelpers.cpp

-              r87e2e39
+              rc111db
  */
+#include <fstream>
+#include "linkedcell.hpp"
+#include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
+#include "vector.hpp"
+#include "verbose.hpp"
 double DetGet(gsl_matrix *A, int inPlace) {
 …
   return result;
+};
+/** Gets the angle between a point and a reference relative to the provided center.
+ * We have two shanks point and reference between which the angle is calculated
+ * and by scalar product with OrthogonalVector we decide the interval.
+ * @param point to calculate the angle for
+ * @param reference to which to calculate the angle
+ * @param OrthogonalVector points in direction of [pi,2pi] interval
+ *
+ * @return angle between point and reference
+ */
+double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
+{
+  if (reference.IsZero())
+    return M_PI;
+  // calculate both angles and correct with in-plane vector
+  if (point.IsZero())
+    return M_PI;
+  double phi = point.Angle(&reference);
+  if (OrthogonalVector.ScalarProduct(&point) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+  cout << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
+  return phi;
+}
+/** Calculates the volume of a general tetraeder.
+ * \param *a first vector
+ * \param *a first vector
+ * \param *a first vector
+ * \param *a first vector
+ * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot  (a-b)) \f$
+ */
+double CalculateVolumeofGeneralTetraeder(Vector *a, Vector *b, Vector *c, Vector *d)
+{
+  Vector Point, TetraederVector[3];
+  double volume;
+  TetraederVector[0].CopyVector(a);
+  TetraederVector[1].CopyVector(b);
+  TetraederVector[2].CopyVector(c);
+  for (int j=0;j<3;j++)
+    TetraederVector[j].SubtractVector(d);
+  Point.CopyVector(&TetraederVector[0]);
+  Point.VectorProduct(&TetraederVector[1]);
+  volume = 1./6. * fabs(Point.ScalarProduct(&TetraederVector[2]));
+  return volume;
+};
+/** 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->triangles.size() < 2) {
+            counter++;
+            break;  // increase counter only once per edge
+          }
+        }
+      } else { // no line
+        cout << Verbose(1) << "The line between " << *nodes[i] << " and " << *nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl;
+        result = true;
+      }
+    }
+  if ((!result) && (counter > 1)) {
+    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
+    result = true;
+  }
+  return result;
+};
+/** Sort function for the candidate list.
+ */
+bool SortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2)
+{
+  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->node);
+  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  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->node);
+  helper.SubtractVector(candidate1->point->node);
+  OrthogonalVector.CopyVector(&helper);
+  helper.VectorProduct(&BaseLineVector);
+  OrthogonalVector.SubtractVector(&helper);
+  OrthogonalVector.Normalize();
+  // calculate both angles and correct with in-plane vector
+  helper.CopyVector(candidate1->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  double phi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+  helper.CopyVector(candidate2->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  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;
+};
+/**
+ * Finds the point which is second closest to the provided one.
+ *
+ * @param Point to which to find the second closest other point
+ * @param linked cell structure
+ *
+ * @return point which is second closest to the provided one
+ */
+TesselPoint* FindSecondClosestPoint(const Vector* Point, LinkedCell* LC)
+{
+  LinkedNodes *List = NULL;
+  TesselPoint* closestPoint = NULL;
+  TesselPoint* secondClosestPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  LC->SetIndexToVector(Point); // 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 << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper. Norm();
+            if (currentNorm < distance) {
+              // remember second point
+              secondDistance = distance;
+              secondClosestPoint = closestPoint;
+              // mark down new closest point
+              distance = currentNorm;
+              closestPoint = (*Runner);
+              //cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *secondClosestPoint << "." << endl;
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  return secondClosestPoint;
+};
+/**
+ * Finds the point which is closest to the provided one.
+ *
+ * @param Point to which to find the closest other point
+ * @param SecondPoint the second closest other point on return, NULL if none found
+ * @param linked cell structure
+ *
+ * @return point which is closest to the provided one, NULL if none found
+ */
+TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC)
+{
+  LinkedNodes *List = NULL;
+  TesselPoint* closestPoint = NULL;
+  SecondPoint = NULL;
+  double distance = 1e16;
+  double secondDistance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  LC->SetIndexToVector(Point); // 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(3) << "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 << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(Point);
+            helper.SubtractVector((*Runner)->node);
+            double currentNorm = helper. Norm();
+            if (currentNorm < distance) {
+              secondDistance = distance;
+              SecondPoint = closestPoint;
+              distance = currentNorm;
+              closestPoint = (*Runner);
+              //cout << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
+            } else if (currentNorm < secondDistance) {
+              secondDistance = currentNorm;
+              SecondPoint = (*Runner);
+              //cout << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  return closestPoint;
+};
+/** Returns the closest point on \a *Base with respect to \a *OtherBase.
+ * \param *out output stream for debugging
+ * \param *Base reference line
+ * \param *OtherBase other base line
+ * \return Vector on reference line that has closest distance
+ */
+Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase)
+{
+  // construct the plane of the two baselines (i.e. take both their directional vectors)
+  Vector Normal;
+  Vector Baseline, OtherBaseline;
+  Baseline.CopyVector(Base->endpoints[1]->node->node);
+  Baseline.SubtractVector(Base->endpoints[0]->node->node);
+  OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
+  OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
+  Normal.CopyVector(&Baseline);
+  Normal.VectorProduct(&OtherBaseline);
+  Normal.Normalize();
+  *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
+  // project one offset point of OtherBase onto this plane (and add plane offset vector)
+  Vector NewOffset;
+  NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
+  NewOffset.SubtractVector(Base->endpoints[0]->node->node);
+  NewOffset.ProjectOntoPlane(&Normal);
+  NewOffset.AddVector(Base->endpoints[0]->node->node);
+  Vector NewDirection;
+  NewDirection.CopyVector(&NewOffset);
+  NewDirection.AddVector(&OtherBaseline);
+  // calculate the intersection between this projected baseline and Base
+  Vector *Intersection = new Vector;
+  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
+  Normal.CopyVector(Intersection);
+  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
+  return Intersection;
+};
+/** Creates the objects in a VRML file.
+ * \param *out output stream for debugging
+ * \param *vrmlfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, PointCloud *cloud)
+{
+  TesselPoint *Walker = NULL;
+  int i;
+  Vector *center = cloud->GetCenter(out);
+  if (vrmlfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *vrmlfile << "#VRML V2.0 utf8" << endl;
+    *vrmlfile << "#Created by molecuilder" << endl;
+    *vrmlfile << "#All atoms as spheres" << endl;
+    cloud->GoToFirst();
+    while (!cloud->IsEnd()) {
+      Walker = cloud->GetPoint();
+      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Walker->node->x[i]-center->x[i] << " ";
+      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+      cloud->GoToNext();
+    }
+    *vrmlfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *vrmlfile << "1" << endl << "  "; // 1 is triangle type
+      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->node->x[j]-center->x[j] << " ";
+        *vrmlfile << "\t";
+      }
+      *vrmlfile << "1. 0. 0." << endl;  // red as colour
+      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
+  }
+  delete(center);
+};
+/** Writes additionally the current sphere (i.e. the last triangle to file).
+ * \param *out output stream for debugging
+ * \param *rasterfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
+{
+  Vector helper;
+  // include the current position of the virtual sphere in the temporary raster3d file
+  Vector *center = cloud->GetCenter(out);
+  // make the circumsphere's center absolute again
+  helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
+  helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
+  helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
+  helper.Scale(1./3.);
+  helper.SubtractVector(center);
+  // and add to file plus translucency object
+  *rasterfile << "# current virtual sphere\n";
+  *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+  *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
+  *rasterfile << "9\n  terminating special property\n";
+  delete(center);
+};
+/** Creates the objects in a raster3d file (renderable with a header.r3d).
+ * \param *out output stream for debugging
+ * \param *rasterfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud)
+{
+  TesselPoint *Walker = NULL;
+  int i;
+  Vector *center = cloud->GetCenter(out);
+  if (rasterfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+    *rasterfile << "@header.r3d" << endl;
+    *rasterfile << "# All atoms as spheres" << endl;
+    cloud->GoToFirst();
+    while (!cloud->IsEnd()) {
+      Walker = cloud->GetPoint();
+      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Walker->node->x[i]-center->x[i] << " ";
+      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+      cloud->GoToNext();
+    }
+    *rasterfile << "# All tesselation triangles" << endl;
+    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
+      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->node->x[j]-center->x[j] << " ";
+        *rasterfile << "\t";
+      }
+      *rasterfile << "1. 0. 0." << endl;  // red as colour
+      //*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
+    }
+    *rasterfile << "9\n#  terminating special property\n";
+  } else {
+    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
+  IncludeSphereinRaster3D(out, rasterfile, Tess, cloud);
+  delete(center);
+};
+/** This function creates the tecplot file, displaying the tesselation of the hull.
+ * \param *out output stream for debugging
+ * \param *tecplot output stream for tecplot data
+ * \param N arbitrary number to differentiate various zones in the tecplot format
+ */
+void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, PointCloud *cloud, int N)
+{
+  if ((tecplot != NULL) && (TesselStruct != NULL)) {
+    // write header
+    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
+    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
+    *tecplot << "ZONE T=\"" << N << "-";
+    for (int i=0;i<3;i++)
+      *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
+    *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+    int i=0;
+    for (cloud->GoToFirst(); !cloud->IsEnd(); cloud->GoToNext(), i++);
+    int *LookupList = new int[i];
+    for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++)
+      LookupList[i] = -1;
+    // print atom coordinates
+    *out << Verbose(2) << "The following triangles were created:";
+    int Counter = 1;
+    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->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
+    }
+    *tecplot << endl;
+    // print connectivity
+    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
+      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
+      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+    }
+    delete[] (LookupList);
+    *out << endl;
+  }
+};
+/** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
+ * Sets BoundaryPointSet::value equal to the number of connected lines that are not convex.
+ * \param *out output stream for debugging
+ * \param *TesselStruct pointer to Tesselation structure
+ */
+void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct)
+{
+  class BoundaryPointSet *point = NULL;
+  class BoundaryLineSet *line = NULL;
+  //*out << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
+  // calculate remaining concavity
+  for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+    point = PointRunner->second;
+    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
+    point->value = 0;
+    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+      line = LineRunner->second;
+      //*out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+      if (!line->CheckConvexityCriterion(out))
+        point->value += 1;
+    }
+  }
+  //*out << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;
+};
+/** Checks whether each BoundaryLineSet in the Tesselation has two triangles.
+ * \param *out output stream for debugging
+ * \param *TesselStruct
+ * \return true - all have exactly two triangles, false - some not, list is printed to screen
+ */
+bool CheckListOfBaselines(ofstream *out, Tesselation *TesselStruct)
+{
+  LineMap::iterator testline;
+  bool result = false;
+  int counter = 0;
+  *out << Verbose(1) << "Check: List of Baselines with not two connected triangles:" << endl;
+  for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
+    if (testline->second->triangles.size() != 2) {
+      *out << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl;
+      counter++;
+    }
+  }
+  if (counter == 0) {
+    *out << Verbose(1) << "None." << endl;
+    result = true;
+  }
+  return result;
+}

src/tesselationhelpers.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+/*********************************************** includes ***********************************/
 // include config.h
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
-#define HULLEPSILON 1e-7
 #include <gsl/gsl_linalg.h>
 …
 #include <gsl/gsl_vector.h>
+#include "vector.hpp"
+#include <iostream>
+#include "defs.hpp"
+/****************************************** forward declarations *****************************/
+class BoundaryPointSet;
+class BoundaryLineSet;
+class BoundaryTriangleSet;
+class LinkedCell;
+class TesselPoint;
+class PointCloud;
+class Tesselation;
+class Vector;
+/********************************************** definitions *********************************/
+#define HULLEPSILON 1e-10
+/********************************************** declarations *******************************/
 double DetGet(gsl_matrix *A, int inPlace);
 …
 double MinIntersectDistance(const gsl_vector * x, void *params);
 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4);
+double CalculateVolumeofGeneralTetraeder(Vector *a, Vector *b, Vector *c, Vector *d);
+double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector);
+bool CheckLineCriteriaForDegeneratedTriangle(class BoundaryPointSet *nodes[3]);
+bool SortCandidates(class CandidateForTesselation* candidate1, class CandidateForTesselation* candidate2);
+TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, LinkedCell* LC);
+TesselPoint* FindSecondClosestPoint(const Vector*, LinkedCell*);
+Vector * GetClosestPointBetweenLine(ofstream *out, class BoundaryLineSet *Base, class BoundaryLineSet *OtherBase);
+void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, PointCloud *cloud, int N);
+void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud);
+void IncludeSphereinRaster3D(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, PointCloud *cloud);
+void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, PointCloud *cloud);
+void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct);
+bool CheckListOfBaselines(ofstream *out, Tesselation *TesselStruct);
 #endif /* TESSELATIONHELPERS_HPP_ */

src/unittests/Makefile.am

-              r87e2e39
+              rc111db
 INCLUDES = -I$(top_srcdir)/src
 noinst_PROGRAMS =  MemoryAllocatorUnitTest MemoryUsageObserverUnitTest VectorUnitTest
+noinst_PROGRAMS =  ActOnAllTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest VectorUnitTest
 TESTS = VectorUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest
+TESTS = ActOnAllTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest VectorUnitTest
 check_PROGRAMS = $(TESTS)
+ActOnAllTest_SOURCES = ActOnAllTest.hpp ActOnAllUnitTest.cpp ActOnAllUnitTest.hpp memoryallocator.hpp
+ActOnAllTest_CXXFLAGS = $(CPPUNIT_CFLAGS)
+ActOnAllTest_LDFLAGS = $(CPPUNIT_LIBS) -ldl
+ActOnAllTest_LDADD = ../libmolecuilder.a
 VectorUnitTest_SOURCES = defs.hpp helpers.hpp leastsquaremin.hpp memoryallocator.hpp memoryusageobserver.hpp vectorunittest.cpp vectorunittest.hpp vector.hpp verbose.hpp

src/unittests/tesselationunittest.cpp

Property mode changed from 100755 to 100644

-              r87e2e39
+              rc111db
   delete(TesselStruct);
   for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
     delete((*Runner)->Name);
+    delete[]((*Runner)->Name);
     delete((*Runner)->node);
     delete(*Runner);
 …
     for (set<BoundaryTriangleSet*>::iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); TriangleRunner++)
       CPPUNIT_ASSERT_EQUAL( true, (*TriangleRunner)->ContainsBoundaryPoint(Walker) );
+    delete(triangles);
+  }
+}

src/unittests/tesselationunittest.hpp

Property mode changed from 100755 to 100644

-              r87e2e39
+              rc111db
 #define TESSELATIONUNITTEST_HPP_
+/*********************************************** includes ***********************************/
 #include <cppunit/extensions/HelperMacros.h>
+#include "linkedcell.hpp"
 #include "tesselation.hpp"

src/unittests/vectorunittest.cpp

-              r87e2e39
+              rc111db
 #include <cppunit/ui/text/TestRunner.h>
+#include "defs.hpp"
+#include "vector.hpp"
 #include "vectorunittest.hpp"
-#include "vector.hpp"
-#include "defs.hpp"
 /********************************************** Test classes **************************************/

src/vector.cpp

-              r87e2e39
+              rc111db
   Direction.CopyVector(LineVector);
   Direction.SubtractVector(Origin);
+  Direction.Normalize();
   //*out << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
 …
   helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
+  if (factor < MYEPSILON) { // Origin is in-plane
+    //*out << Verbose(2) << "Origin of line is in-plane, simple." << endl;
+    CopyVector(Origin);
+    return true;
+  }
   //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
   Direction.Scale(factor);
 …
 };
+/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
+ * \param *M matrix of box
+ * \param *Minv inverse matrix
+ */
+void Vector::WrapPeriodically(const double *M, const double *Minv)
+{
+  MatrixMultiplication(Minv);
+  // truncate to [0,1] for each axis
+  for (int i=0;i<NDIM;i++) {
+    x[i] += 0.5;  // set to center of box
+    while (x[i] >= 1.)
+      x[i] -= 1.;
+    while (x[i] < 0.)
+      x[i] += 1.;
+  }
+  MatrixMultiplication(M);
+};
 /** Do a matrix multiplication.
  * \param *matrix NDIM_NDIM array
  */
 void Vector::MatrixMultiplication(double *M)
+void Vector::MatrixMultiplication(const double *M)
+{
   Vector C;
 …
  * \param *matrix NDIM_NDIM array
  */
 void Vector::InverseMatrixMultiplication(double *A)
+void Vector::InverseMatrixMultiplication(const double *A)
+{
   Vector C;

src/vector.hpp

-              r87e2e39
+              rc111db
 using namespace std;
+#include "helpers.hpp"
+/*********************************************** includes ***********************************/
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_multimin.h>
+class Vector;
+#include "defs.hpp"
+/********************************************** declarations *******************************/
 /** Single vector.
 …
   void Scale(double *factor);
   void Scale(double factor);
   void MatrixMultiplication(double *M);
+  void MatrixMultiplication(const double *M);
   double * InverseMatrix(double *A);
   void InverseMatrixMultiplication(double *M);
+  void InverseMatrixMultiplication(const double *M);
   void KeepPeriodic(ofstream *out, double *matrix);
   void LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors);
 …
   bool Output(ofstream *out) const;
   bool IsInParallelepiped(Vector offset, double *parallelepiped);
+  void WrapPeriodically(const double *M, const double *Minv);
 };

src/verbose.cpp

r87e2e39	rc111db
1		#include "molecules.hpp"
	1	using namespace std;
	2
	3	#include "verbose.hpp"
2	4
3	5	/** Prints the tabs according to verbosity stored in the temporary constructed class.

src/verbose.hpp

-              r87e2e39
+              rc111db
 #ifndef VERBOSE_HPP_
 #define VERBOSE_HPP_
+using namespace std;
+/*********************************************** includes ***********************************/
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <iostream>
 /************************************* Class Verbose & Binary *******************************/

tests/Makefile.am

r87e2e39	rc111db
1	1	AUTOM4TE = autom4te
2		EXTRA_DIST = testsuite.at $(TESTSUITE) atlocal.in
	2	EXTRA_DIST = testsuite.at $(TESTSUITE) atlocal.in regression
3	3	TESTSUITE = $(srcdir)/testsuite
4	4

tests/testsuite.at

-              r87e2e39
+              rc111db
 # Process with autom4te to create an -*- Autotest -*- test suite.
+#
+# see regression/... subdirs wherein for each of the cases and each check (enumerated) is a pre and post dir.
+# In pre initial files are placed, in post results can be found to be checked by diff in this testsuite.
 AT_INIT([Molecular Builder])
 …
 AT_BANNER([MoleCuilder - standard options])
 AT_SETUP([Standard Options])
+AT_KEYWORDS([options])
 AT_CHECK([pwd],[ignore],[ignore])
 AT_CHECK([../../molecuilder -v], 0, [stdout], [ignore])
 …
 AT_CLEANUP
 AT_BANNER([MoleCuilder - molecular config creation from xyz file and atom adding])
+AT_SETUP([Simple configuration])
+# 1. create a fake element database with the only element we need
+AT_DATA([elements.db],[#        Covalent        radius  of      each    element in      Angstroem       from    CSD     (binding        is:     [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t]   with    t       =       0.4A
+#Element        Name    Symbol  Period  Group   Block   Atomic  Number  AtomicWeight    Covalent        Radius  vdW     Radius
+Hydrogen        H       1       1       s       1       1.008   0.23    1.09
+])
+# 2. create some simplest molecular geometry
+AT_KEYWORDS([Atom handling])
+# 1. create some simplest molecular geometry
+AT_SETUP([Simple configuration - xyz file generation])
 AT_DATA([test.xyz], [[1
  # test configuration, created by molecuilder test suite
 H       10.     10.     10.
 ]])
+# 3. make sure config is empty and not remnant from last test with broken dirs
+AT_DATA([test.conf], [])
+AT_CHECK([../../molecuilder test.conf -e ./ -p test.xyz], 0, [ignore], [ignore])
+AT_CHECK([fgrep "Ion_Type1_1" test.conf], 0, [Ion_Type1_1       10.000000000    10.000000000    10.000000000    0 # Number in molecule 0
+AT_CHECK([file=test.xyz; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/1/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP
+# 2. parsing an xyz
+AT_SETUP([Simple configuration - parsing xyz file])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/2/pre/test.xyz .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -p test.xyz], 0, [ignore], [ignore])
+AT_CHECK([fgrep "Ion_Type1_1" test.conf], 0, [Ion_Type1_1       10.000000000    10.000000000    10.000000000    0 # molecule nr 0
 ], [ignore])
+AT_CHECK([cp test.conf main_pcp_linux], 0, [ignore], [ignore])
+AT_DATA([input], [aa 10. 10. 10. 1
+s
+q
+])
+AT_CHECK([../../molecuilder -e ./ <input], 0, [ignore], [ignore])
+AT_CHECK([diff main_pcp_linux test.conf], 0, [ignore], [ignore])
+# 4. test some more configuration
+AT_CHECK([../../molecuilder test.conf -e ./ -t -s -b -F -E -c -b -a -U -T -u], 0, [ignore], [stderr])
+AT_CHECK([fgrep -c "Not enough or invalid" stderr], 0, [10
+AT_CHECK([file=test.conf; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/2/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf.in; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/2/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf.xyz; diff -I '.*Created by molecuilder.*' $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/2/post/$file], 0, [ignore], [ignore])
+# 3. add atom
+AT_CLEANUP
+AT_SETUP([Simple configuration - adding atom])
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -a 1 10. 10. 10.], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/3/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf.in; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/3/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf.xyz; diff -I '.*Created by molecuilder.*' $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/3/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP
+# 4. change the element
+AT_SETUP([Simple configuration - Changing element])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/4/pre/test.conf test.conf], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -E 0 6], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/4/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP
+# 5. remove atom
+AT_SETUP([Simple configuration - Atom removal])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/5/pre/test.conf .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -r 0], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/5/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf.in; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/5/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=test.conf.xyz; diff -I '.*Created by molecuilder.*' $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/5/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP
+# 6. test some more configuration that all desire parameters and count how many complain
+AT_SETUP([Simple configuration - invalid commands on empty configs])
+AT_CHECK([../../molecuilder empty.conf -e ${abs_top_srcdir}/src/ -t -s -b -E -c -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "Not enough or invalid" stderr], 0, [1
+], [ignore])
+AT_CLEANUP
+# 7. test some more configuration that all need parameters and count how many complain
+AT_SETUP([Simple configuration - invalid commands on present configs])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Simple_configuration/7/pre/test.conf .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -t -s -b -E -c -b -a -U -T -u], 255, [ignore], [stderr])
+AT_CHECK([fgrep -c "Not enough or invalid" stderr], 0, [9
+], [ignore])
+AT_CLEANUP
+AT_BANNER([MoleCuilder - Graph routines test])
+AT_KEYWORDS([graph])
+AT_SETUP([Graph - DFS analysis])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Graph/1/pre/test.conf .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -D 2.], 0, [stdout], [stderr])
+AT_CHECK([fgrep -c "No rings were detected in the molecular structure." stdout], 0, [1
 ], [ignore])
 AT_CLEANUP
 AT_BANNER([MoleCuilder - Fragmentation and Re-fragmentation test])
+AT_SETUP([Fragmentation])
+# 1. create a fake element database with the only two elements we need
+AT_DATA([elements.db],[#  Covalent  radius  of  each  element in  Angstroem from  CSD (binding  is: [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with  t = 0.4A
+#Element  Name  Symbol  Period  Group Block Atomic  Number  AtomicWeight  Covalent  Radius  vdW Radius
+Hydrogen  H 1 1 s 1 1.008 0.23  1.09
+Carbon  C       2       14      p       6       12.011  0.68    1.70
+])
+# 2. create molecular geometry
+AT_DATA([test.xyz], [[11
+ # test configuration, created by molecuilder test suite
+C       9.782085945     3.275186040     3.535886037
+C       8.532785963     4.158586027     3.535886037
+C       7.283585982     3.275186040     3.535886037
+H       9.782085945     2.645886050     2.645886050
+H       9.782085945     2.645886050     4.425886024
+H       10.672039608    3.904536878     3.535886037
+H       8.532785963     4.787886018     2.645886050
+H       8.532785963     4.787886018     4.425886024
+H       6.393632318     3.904536877     3.535886037
+H       7.283585982     2.645886050     2.645886050
+H       7.283585982     2.645886050     4.425886024
+]])
+# 3. make sure config is empty and not remnant from last test with broken dirs
+AT_DATA([test.conf], [])
+# 4. create the config and check it
+AT_CHECK([../../molecuilder test.conf -e ./ -p test.xyz], 0, [ignore], [ignore])
+AT_CHECK([fgrep "Ion_Type1_4" test.conf], 0, [Ion_Type1_4       8.532785963     4.787886018     2.645886050     0 # Number in molecule 6
+AT_KEYWORDS([fragmentation])
+# 1. check config
+AT_SETUP([Fragmentation - Checking present config])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Fragmentation/1/pre/test.conf .], 0)
+AT_CHECK([fgrep "Ion_Type1_4" test.conf], 0, [Ion_Type1_4       8.532785963     4.787886018     2.645886050     0 # molecule nr 6
 ], [ignore])
 AT_CHECK([fgrep "Ion_Type2_4" test.conf], 1, [ignore], [ignore])
+# 5. fragment the molecule and check the number of configs
+AT_CHECK([../../molecuilder test.conf -e ./ -f 1.55 2], 0, [ignore], [ignore], [mkdir std; mv BondFragment*.conf* std/])
+AT_CLEANUP
+# 2. fragment the molecule and check the number of configs
+AT_SETUP([Fragmentation - Fragmentation])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Fragmentation/2/pre/test.conf .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -f 1.55 2], 0, [ignore], [ignore], [mkdir std; mv BondFragment*.conf* std/])
 AT_CHECK([mkdir std; mv BondFragment*.conf* std/], 0)
 AT_CHECK([ls -l std/BondFragment*.conf | wc -l], 0, [5
 ], [ignore])
+# 3a. check whether parsing of BondFragment files and re-rwriting config files is working (exit code is 2 as we don't need to continue wrt to ...OrderAtSite)
+AT_CHECK([../../molecuilder test.conf -e ./ -f 1.55 2], 2, [ignore], [ignore], [mkdir new; mv BondFragment*.conf* new/])
+AT_CLEANUP
+# 3. check whether parsing of BondFragment files and re-rwriting config files is working (exit code is 2 as we don't need to continue wrt to ...OrderAtSite)
+AT_SETUP([Fragmentation - Fragmentation is at MaxOrder])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Fragmentation/3/pre/* .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -f 1.55 2], 2, [ignore], [ignore], [mkdir new; mv BondFragment*.conf* new/])
 AT_CHECK([mkdir new; mv BondFragment*.conf* new/], 0)
-# 6. compare both dirs by diff'ing
-AT_CHECK([diff -I '.*Created by molecuilder.*' std/ new/], 0, [], [])
 AT_CLEANUP
+AT_BANNER([MoleCuilder - Tesselation test])
+AT_KEYWORDS([Tesselation])
+# 1. Non convex tesselation
+AT_SETUP([Tesselation - Non-Convex Envelope])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/1/pre/* .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -N 4. NonConvexEnvelope], 0, [stdout], [stderr])
+AT_CHECK([file=NonConvexEnvelope.dat; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/1/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=NonConvexEnvelope.r3d; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/1/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP
+# 2. convex tesselation (where the non-convex is already convex)
+AT_SETUP([Tesselation - Convex Envelope])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/2/pre/* .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -o ConvexEnvelope NonConvexEnvelope], 0, [stdout], [stderr])
+AT_CHECK([file=ConvexEnvelope.dat; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/2/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=ConvexEnvelope.r3d; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/2/post/$file], 0, [ignore], [ignore])
+AT_CHECK([fgrep "RESULT: The summed volume is 16.401577 angstrom^3" stdout], 0, [ignore], [ignore])
+AT_CHECK([diff ConvexEnvelope.dat NonConvexEnvelope.dat], 0, [ignore], [ignore])
+AT_CLEANUP
+# 3. Big Non convex tesselation
+AT_SETUP([Tesselation - Big non-Convex Envelope])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/3/pre/* .], 0)
+AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -N 4. NonConvexEnvelope], 0, [stdout], [stderr])
+AT_CHECK([file=NonConvexEnvelope.dat; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/3/post/$file], 0, [ignore], [ignore])
+AT_CHECK([file=NonConvexEnvelope.r3d; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/3/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP
+# 4. Big convex tesselation - is not working yet
+#AT_SETUP([Tesselation - big convex Envelope])
+#AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/4/pre/* .], 0)
+#AT_CHECK([../../molecuilder test.conf -e ${abs_top_srcdir}/src/ -o ConvexEnvelope NonConvexEnvelope], 0, [stdout], [stderr])
+#AT_CHECK([file=ConvexEnvelope.dat; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/4/post/$file], 0, [ignore], [ignore])
+#AT_CHECK([file=ConvexEnvelope.r3d; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/regression/Tesselation/4/post/$file], 0, [ignore], [ignore])
+#AT_CHECK([fgrep "RESULT: The summed volume is 16.401577 angstrom^3" stdout], 0, [ignore], [ignore])
+#AT_CLEANUP

Context Navigation

Changes in / [87e2e39:c111db]

Legend:

Download in other formats: