Changes in / [c0917c:a33931]

Location:

src

Files:

• Makefile.am (modified) (2 diffs)
• atom.cpp (modified) (1 diff)
• boundary.cpp (modified) (14 diffs)
• boundary.hpp (modified) (2 diffs)
• builder.cpp (modified) (8 diffs)
• config.cpp (modified) (1 diff)
• linkedcell.cpp (modified) (1 diff)
• linkedcell.hpp (modified) (1 diff)
• moleculelist.cpp (modified) (2 diffs)
• molecules.cpp (modified) (5 diffs)
• molecules.hpp (modified) (1 diff)
• tesselation.cpp (modified) (74 diffs)
• tesselation.hpp (modified) (7 diffs)
• tesselationhelpers.cpp (modified) (1 diff)
• tesselationhelpers.hpp (modified) (3 diffs)
• tesselationunittest.cpp (modified) (1 prop)
• tesselationunittest.hpp (modified) (1 prop)
• vector.cpp (modified) (2 diffs)
• vectorunittest.cpp (modified) (1 diff)

src/Makefile.am

@@ -2 +2 @@
 HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp memoryallocator.hpp memoryusageobserver.cpp memoryusageobserver.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
 
-noinst_PROGRAMS =  VectorUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest
+noinst_PROGRAMS =  VectorUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest TesselationUnitTest
 
 bin_PROGRAMS = molecuilder joiner analyzer
@@ -11 +11 @@
 analyzer_SOURCES = analyzer.cpp datacreator.cpp element.cpp helpers.cpp memoryallocator.hpp memoryusageobserver.cpp memoryusageobserver.hpp periodentafel.cpp parser.cpp verbose.cpp helpers.hpp periodentafel.hpp parser.hpp datacreator.hpp 
 
-TESTS = VectorUnitTest MemoryAllocatorUnitTest
+TESTS = VectorUnitTest MemoryAllocatorUnitTest MemoryUsageObserverUnitTest TesselationUnitTest
 check_PROGRAMS = $(TESTS)
 VectorUnitTest_SOURCES = defs.hpp helpers.cpp helpers.hpp leastsquaremin.cpp leastsquaremin.hpp memoryallocator.hpp memoryusageobserver.cpp memoryusageobserver.hpp vectorunittest.cpp vectorunittest.hpp vector.cpp vector.hpp verbose.cpp verbose.hpp 
 VectorUnitTest_CXXFLAGS = $(CPPUNIT_CFLAGS)
 VectorUnitTest_LDFLAGS = $(CPPUNIT_LIBS) -ldl
+TesselationUnitTest_SOURCES = defs.hpp helpers.cpp helpers.hpp leastsquaremin.cpp leastsquaremin.hpp linkedcell.cpp linkedcell.hpp tesselation.cpp tesselation.hpp tesselationhelpers.cpp tesselationhelpers.hpp tesselationunittest.cpp tesselationunittest.hpp vector.cpp vector.hpp verbose.cpp verbose.hpp
+TesselationUnitTest_CXXFLAGS = $(CPPUNIT_CFLAGS)
+TesselationUnitTest_LDFLAGS = $(CPPUNIT_LIBS) -ldl
 
 MemoryAllocatorUnitTest_SOURCES = defs.hpp helpers.cpp helpers.hpp memoryallocatorunittest.cpp memoryallocatorunittest.hpp memoryallocator.hpp memoryusageobserver.cpp memoryusageobserver.hpp verbose.cpp verbose.hpp 

src/atom.cpp

@@ -62 +62 @@
 atom::~atom()
 {
-  Free(&ComponentNr);
+  Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
+  Free<char *>(&Name, "atom::~atom: *Name");
 };
 

src/boundary.cpp

@@ -113 +113 @@
 ;
 
-/** 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.
@@ -537 +391 @@
 
         // 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;
+        }
       }
     }
@@ -577 +433 @@
 
   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;
 };
 
@@ -609 +501 @@
   class BoundaryLineSet *line = NULL;
   bool Concavity;
+  char dummy[MAXSTRINGSIZE];
   PointMap::iterator PointRunner, PointAdvance;
   LineMap::iterator LineRunner, LineAdvance;
@@ -621 +514 @@
   }
 
-  //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
@@ -636 +532 @@
         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
@@ -658 +558 @@
       // 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.
@@ -794 +647 @@
 
   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.
@@ -958 +841 @@
     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++;
@@ -1080 +963 @@
  * \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;
@@ -1103 +982 @@
   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);
@@ -1111 +991 @@
   }
 
+  // 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;
@@ -1184 +1019 @@
     }
 
-    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)
@@ -1265 +1058 @@
   *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

@@ -16 +16 @@
 
 #define DEBUG 1
-#define DoSingleStepOutput 0
+#define DoSingleStepOutput 1
 #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* >
@@ -38 +32 @@
 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

@@ -1018 +1018 @@
       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;
@@ -1075 +1076 @@
 
     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;
 
@@ -1626 +1642 @@
                 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;
               }
@@ -1654 +1670 @@
             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 ((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':
@@ -1710 +1733 @@
                 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]);
@@ -1716 +1739 @@
                 argptr+=3;
               }
+              break;
             case 'T':
               ExitFlag = 1;
@@ -1724 +1748 @@
                 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]);
@@ -1865 +1889 @@
             case 'o':
               ExitFlag = 1;
-              if ((argptr >= argc) || (argv[argptr][0] == '-')){
+              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;

src/config.cpp

@@ -1011 +1011 @@
       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/linkedcell.cpp

@@ -96 +96 @@
     //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

@@ -44 +44 @@
     LinkedCell();
     LinkedCell(PointCloud *set, double RADIUS);
+    LinkedCell(LinkedNodes *set, double radius);
     ~LinkedCell();
     LinkedNodes* GetCurrentCell();

src/moleculelist.cpp

@@ -5 +5 @@
  */
 
+#include "boundary.hpp"
 #include "config.hpp"
 #include "molecules.hpp"
@@ -292 +293 @@
 /** 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/molecules.cpp

@@ -687 +687 @@
 {
   int length = 0;
-  char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
+  const char *molname = strrchr(filename, '/');
+  if (molname != NULL)
+    molname += sizeof(char);  // search for filename without dirs
+  else
+    molname = filename; // contains no slashes
   char *endname = strchr(molname, '.');
   if ((endname == NULL) || (endname < molname))
@@ -1193 +1197 @@
 /** Evaluates the potential energy used for constrained molecular dynamics.
  * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \frac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$
- *     where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}$ is minimum distance between 
+ *     where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}\f$ is minimum distance between
  *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
  * Note that for the second term we have to solve the following linear system:
@@ -1587 +1591 @@
  * \param endstep stating final configuration in molecule::Trajectories
  * \param &config configuration structure
+ * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential()
  * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
  */
-bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration) 
+bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity)
 {
   molecule *mol = NULL;
@@ -1598 +1603 @@
   atom **PermutationMap = NULL;
   atom *Walker = NULL, *Sprinter = NULL;
-  MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
+  if (!MapByIdentity)
+    MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
+  else {
+    PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
+    Walker = start;
+    while (Walker->next != end) {
+      Walker = Walker->next;
+      PermutationMap[Walker->nr] = Walker;   // always pick target with the smallest distance
+    }
+  }
 
   // check whether we have sufficient space in Trajectories for each atom
@@ -3747 +3761 @@
  * \param *&SortIndex Mapping array of size molecule::AtomCount
  * \return true - success, false - failure of SortIndex alloc
+ * \todo do we really need this still as the IonType may appear in any order due to recent changes
  */
 bool molecule::CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex)

src/molecules.hpp

@@ -182 +182 @@
   double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
   void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
-  bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration); 
+  bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity);
         
   bool CheckBounds(const Vector *x) const;

src/tesselation.cpp

@@ -6 +6 @@
  */
 
+#include "helpers.hpp"
+#include "linkedcell.hpp"
+#include "memoryallocator.hpp"
 #include "tesselation.hpp"
-#include "memoryallocator.hpp"
+#include "tesselationhelpers.hpp"
+#include "vector.hpp"
+
+class molecule;
 
 // ======================================== Points on Boundary =================================
@@ -37 +43 @@
 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;
@@ -67 +73 @@
 ostream & operator <<(ostream &ost, BoundaryPointSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.node->Name << " at " << *a.node->node << "]";
   return ost;
 }
@@ -125 +131 @@
         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]);
@@ -168 +175 @@
 
 /** 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
@@ -178 +185 @@
   // 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;
   }
@@ -201 +208 @@
     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) {
@@ -211 +223 @@
       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;
     }
@@ -217 +229 @@
   //*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;
   }
@@ -262 +275 @@
 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;
 };
@@ -332 +345 @@
   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;
@@ -341 +355 @@
     }
   }
-  cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+  //cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
 };
 
@@ -361 +375 @@
  * 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
@@ -395 +409 @@
     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!
@@ -426 +440 @@
   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;
@@ -451 +477 @@
 };
 
+/** 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
@@ -485 +532 @@
 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;
 };
@@ -505 +552 @@
 TesselPoint::~TesselPoint()
 {
-  Free(&Name);
+  Free<char *>(&Name, "~TesselPoint - *Name");
 };
 
@@ -512 +559 @@
 ostream & operator << (ostream &ost, const TesselPoint &a)
 {
-  ost << "[" << (a.Name) << "|" << &a << "]";
+  ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
   return ost;
 };
@@ -569 +616 @@
   TrianglesOnBoundaryCount = 0;
   InternalPointer = PointsOnBoundary.begin();
+  LastTriangle = NULL;
+  TriangleFilesWritten = 0;
 }
 ;
@@ -585 +634 @@
       cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
   }
+  cout << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
 }
 ;
@@ -1052 +1102 @@
   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;
@@ -1090 +1147 @@
         // add Walker to boundary points
         *out << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl;
+        AddFlag = true;
         if (AddBoundaryPoint(Walker,0))
           NewPoint = BPS[0];
@@ -1180 +1238 @@
   } 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;
   }
@@ -1197 +1255 @@
 
   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];
@@ -1232 +1291 @@
 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;
@@ -1242 +1301 @@
 };
 
-/** 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()
@@ -1252 +1311 @@
   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
@@ -1272 +1351 @@
           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;
@@ -1294 +1383 @@
   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;
@@ -1321 +1410 @@
         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;
@@ -1390 +1483 @@
           }
           // 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;
         }
       }
@@ -1400 +1493 @@
   *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;
 };
 
@@ -1461 +1603 @@
   CandidateList *OptCandidates = new CandidateList();
   for (int k=0;k<NDIM;k++) {
+    Oben.Zero();
     Oben.x[k] = 1.;
     FirstPoint = MaxPoint[k];
@@ -1469 +1612 @@
     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);
@@ -1496 +1637 @@
 
     // 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);
@@ -1569 +1712 @@
  * @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";
@@ -1604 +1747 @@
     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
@@ -1691 +1834 @@
         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);
@@ -1732 +1875 @@
     // 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);
+    }
+
   }
 
@@ -1750 +1898 @@
  * \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)
@@ -1827 +1975 @@
  * \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;
@@ -1852 +2001 @@
   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
@@ -1866 +2018 @@
     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++) {
@@ -1876 +2028 @@
     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.;
+    }
+  }
 };
 
@@ -1930 +2042 @@
  * \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;
@@ -1946 +2058 @@
   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++) {
@@ -1982 +2094 @@
   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;
     }
 
@@ -2024 +2136 @@
     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;
 
@@ -2032 +2144 @@
     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;
 };
 
@@ -2046 +2158 @@
 /** 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
@@ -2053 +2164 @@
  * \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;
@@ -2191 +2303 @@
   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
@@ -2216 +2328 @@
     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
+      //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 
       // check SearchDirection
@@ -2392 +2505 @@
   TesselPoint *trianglePoints[3];
   TesselPoint *SecondPoint = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
 
   if (LinesOnBoundary.empty()) {
@@ -2402 +2516 @@
   // 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
@@ -2443 +2581 @@
   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;
@@ -2466 +2616 @@
 {
   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;
+  }
 }
 
@@ -2483 +2645 @@
 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);
 }
 
@@ -2496 +2652 @@
  * @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;
@@ -2510 +2666 @@
     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;
@@ -2521 +2677 @@
   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++;
   }
 
@@ -2543 +2699 @@
     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
@@ -2553 +2711 @@
  *
  * @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()
@@ -2586 +2741 @@
 
   // 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);
@@ -2594 +2752 @@
 
   // 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
@@ -2625 +2983 @@
   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) {
@@ -2645 +3001 @@
     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;
+};
+
+
 
 /**
@@ -2939 +3230 @@
               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);
                 }
@@ -2970 +3248 @@
 
 /**
+ * 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.
  *
@@ -2975 +3290 @@
  *         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;
 
@@ -3026 +3333 @@
 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;
@@ -3044 +3364 @@
       && (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
@@ -3055 +3408 @@
     }
   }
+  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

@@ -34 +34 @@
 class PointCloud;
 class Tesselation;
+
+#define DoTecplotOutput 1
+#define DoRaster3DOutput 1
+#define DoVRMLOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
 
 // ======================================================= some template functions =========================================
@@ -115 +122 @@
 
     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];
@@ -196 +205 @@
     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
@@ -215 +225 @@
     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);
@@ -235 +250 @@
     void PrintAllBoundaryTriangles(ofstream *out);
 
+    // store envelope in file
+    void Output(ofstream *out, const char *filename, PointCloud *cloud);
 
     PointMap PointsOnBoundary;
@@ -257 +274 @@
     class BoundaryLineSet *BLS[3];
     class BoundaryTriangleSet *BTS;
+    class BoundaryTriangleSet *LastTriangle;
+    int TriangleFilesWritten;
 
   private:
@@ -264 +283 @@
 };
 
-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

@@ -386 +386 @@
 
   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

@@ -18 +18 @@
 #endif
 
-#define HULLEPSILON 1e-7
-
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
@@ -26 +24 @@
 #include <gsl/gsl_vector.h>
 
+#include "defs.hpp"
+#include "tesselation.hpp"
 #include "vector.hpp"
+
+#define HULLEPSILON 1e-10
 
 double DetGet(gsl_matrix *A, int inPlace);
@@ -35 +37 @@
 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/vector.cpp

@@ -225 +225 @@
   Direction.CopyVector(LineVector);
   Direction.SubtractVector(Origin);
+  Direction.Normalize();
   //*out << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
@@ -234 +235 @@
   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);

src/vectorunittest.cpp

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