Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 91e7e4a510c48a116682e9998422b928cf2f1117)
+++ src/boundary.cpp	(revision 57066a6de6f2d289266e1e6e8466859414ebf216)
@@ -112,150 +112,4 @@
 ;
 
-/** 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.
@@ -536,8 +390,10 @@
 
         // 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;
+        }
       }
     }
@@ -644,4 +500,5 @@
   class BoundaryLineSet *line = NULL;
   bool Concavity;
+  char dummy[MAXSTRINGSIZE];
   PointMap::iterator PointRunner, PointAdvance;
   LineMap::iterator LineRunner, LineAdvance;
@@ -656,10 +513,13 @@
   }
 
-  //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
@@ -671,16 +531,20 @@
         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
@@ -693,12 +557,17 @@
       // 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");
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  //StoreTrianglesinFile(out, mol, filename, "-third");
 
   // third step: IsConvexRectangle
@@ -717,5 +586,5 @@
       point = TesselStruct->IsConvexRectangle(out, line);
       if (point != NULL)
-        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+        volume += TesselStruct->RemovePointFromTesselatedSurface(out, point);
     }
     LineRunner = LineAdvance;
@@ -723,7 +592,8 @@
 
   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;
@@ -993,5 +863,5 @@
     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++;
@@ -1115,14 +985,11 @@
  * \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;
@@ -1141,4 +1008,5 @@
   bool failflag = false;
 
+  // initialise Linked Cell
   if (LCList == NULL) {
     LCList = new LinkedCell(mol, 2.*RADIUS);
@@ -1146,71 +1014,25 @@
   }
 
+  // 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++;
+  baseline++; // skip first line
   while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (flag)) {
     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.
+      // 3. find next triangle
+      failflag = 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.
       flag = flag || failflag;
       if (!failflag)
         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) || (mol->TesselStruct->TrianglesOnBoundary.size() > 8520 && mol->TesselStruct->TrianglesOnBoundary.size() <= 8530)))) { // 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();
     } else {
       //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->triangles.size() << " triangles adjacent" << endl;
@@ -1220,36 +1042,28 @@
 
     N++;
-    baseline++;
     if ((baseline == mol->TesselStruct->LinesOnBoundary.end()) && (flag)) {
       baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+      baseline++;
       flag = false;
     }
   }
+  // 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();
-
-  // 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;
@@ -1264,35 +1078,7 @@
     *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;
+  // write final envelope
+  CalculateConcavityPerBoundaryPoint(out, mol->TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "");
 
   if (freeLC)
@@ -1300,4 +1086,5 @@
   *out << Verbose(0) << "End of FindNonConvexBorder\n";
 };
+
 
 /** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 91e7e4a510c48a116682e9998422b928cf2f1117)
+++ src/boundary.hpp	(revision 57066a6de6f2d289266e1e6e8466859414ebf216)
@@ -16,12 +16,6 @@
 
 #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,5 +32,5 @@
 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);
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision 91e7e4a510c48a116682e9998422b928cf2f1117)
+++ src/tesselation.cpp	(revision 57066a6de6f2d289266e1e6e8466859414ebf216)
@@ -6,5 +6,11 @@
  */
 
+#include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "tesselation.hpp"
+#include "tesselationhelpers.hpp"
+#include "vector.hpp"
+
+class molecule;
 
 // ======================================== Points on Boundary =================================
@@ -36,5 +42,5 @@
 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;
@@ -66,5 +72,5 @@
 ostream & operator <<(ostream &ost, BoundaryPointSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  ost << "[" << a.Nr << "|" << a.node->Name << " at " << *a.node->node << "]";
   return ost;
 }
@@ -124,14 +130,15 @@
         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]);
@@ -167,5 +174,5 @@
 
 /** 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
@@ -177,5 +184,5 @@
   // 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;
   }
@@ -200,29 +207,35 @@
     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) {
-      //*out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
+      *out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
       helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
-      //*out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
+      *out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
     } else {
-      //*out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
+      //*out << Verbose(2) << "ERROR: I cannot find third node in triangle, something's wrong." << endl;
       return true;
     }
   }
-  //*out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
+  *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;
   }
@@ -261,5 +274,5 @@
 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;
 };
@@ -331,8 +344,9 @@
   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;
@@ -340,5 +354,5 @@
     }
   }
-  cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
+  //cout << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl;
 };
 
@@ -450,4 +464,25 @@
 };
 
+/** 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
@@ -484,6 +519,6 @@
 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;
 };
@@ -504,5 +539,4 @@
 TesselPoint::~TesselPoint()
 {
-  Free((void **)&Name, "TesselPoint::~TesselPoint: *Name");
 };
 
@@ -511,5 +545,5 @@
 ostream & operator << (ostream &ost, const TesselPoint &a)
 {
-  ost << "[" << (a.Name) << "|" << &a << "]";
+  ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]";
   return ost;
 };
@@ -568,4 +602,6 @@
   TrianglesOnBoundaryCount = 0;
   InternalPointer = PointsOnBoundary.begin();
+  LastTriangle = NULL;
+  TriangleFilesWritten = 0;
 }
 ;
@@ -584,4 +620,5 @@
       cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
   }
+  cout << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;
 }
 ;
@@ -1232,5 +1269,5 @@
 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
 {
-  cout << Verbose(4) << "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;
@@ -1252,4 +1289,7 @@
   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
@@ -1524,4 +1564,5 @@
   CandidateList *OptCandidates = new CandidateList();
   for (int k=0;k<NDIM;k++) {
+    Oben.Zero();
     Oben.x[k] = 1.;
     FirstPoint = MaxPoint[k];
@@ -1532,10 +1573,8 @@
     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);
@@ -1559,5 +1598,7 @@
 
     // 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);
@@ -1632,5 +1673,5 @@
  * @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";
@@ -1667,5 +1708,5 @@
     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
@@ -1754,5 +1795,5 @@
         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);
@@ -1795,4 +1836,9 @@
     // 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);
+    }
+
   }
 
@@ -1813,5 +1859,5 @@
  * \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)
@@ -1890,10 +1936,11 @@
  * \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;
@@ -1915,4 +1962,7 @@
   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
@@ -1929,5 +1979,5 @@
     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++) {
@@ -1939,51 +1989,11 @@
     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.;
+    }
+  }
 };
 
@@ -1993,7 +2003,7 @@
  * \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;
@@ -2009,5 +2019,5 @@
   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++) {
@@ -2045,15 +2055,15 @@
   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;
     }
 
@@ -2099,9 +2109,9 @@
   } 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;
 };
 
@@ -2109,5 +2119,4 @@
 /** 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
@@ -2116,8 +2125,9 @@
  * \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;
@@ -2254,5 +2264,5 @@
   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
@@ -2279,4 +2289,5 @@
     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
+      //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 
       // check SearchDirection
@@ -2455,4 +2466,5 @@
   TesselPoint *trianglePoints[3];
   TesselPoint *SecondPoint = NULL;
+  list<BoundaryTriangleSet*> *triangles = NULL;
 
   if (LinesOnBoundary.empty()) {
@@ -2465,28 +2477,54 @@
   // 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*> *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) << "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(2) << "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(1) << "List of triangle points:" << endl;
+      *out << Verbose(2) << *trianglePoints[i] << endl;
+    }
+
+    triangles = FindTriangles(trianglePoints);
+    *out << Verbose(1) << "List of possible triangles:" << endl;
+    for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      *out << Verbose(2) << **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
@@ -2504,13 +2542,25 @@
   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;
@@ -2527,10 +2577,22 @@
 {
   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;
+  }
 }
 
@@ -2544,11 +2606,5 @@
 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);
 }
 
@@ -2705,15 +2761,20 @@
   }
 
-  map <class BoundaryLineSet *, bool> Touched;
-  map <class BoundaryLineSet *, bool>::iterator line;
-  for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++)
-    Touched.insert( pair <class BoundaryLineSet *, bool>(runner->second, false) );
+  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++) {
-      line = Touched.find(runner->second);
-      if (line == Touched.end()) {
+      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 (!line->second) {
-        line->second = true;
+      } else if (!LineRunner->second) {
+        LineRunner->second = true;
         connectedPath = new list<TesselPoint*>;
         triangle = NULL;
@@ -2728,23 +2789,39 @@
 
           // find next triangle
-          for (TriangleMap::iterator TriangleRunner = CurrentLine->triangles.begin(); TriangleRunner != CurrentLine->triangles.end(); TriangleRunner++) {
-            if (TriangleRunner->second != triangle) { // look for first triangle not equal to old one
-              triangle = TriangleRunner->second;
-              break;
+          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;
             }
           }
-          line = Touched.find(CurrentLine);
-          if (line == Touched.end())
+          LineRunner = TouchedLine.find(CurrentLine);
+          if (LineRunner == TouchedLine.end())
             *out << Verbose(2) << "ERROR: I could not find " << *CurrentLine << " in the touched list." << endl;
           else
-            line->second = true;
+            LineRunner->second = true;
           // find next point
           CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
@@ -2753,5 +2830,6 @@
         // last point is missing, as it's on start line
         *out << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;
-        connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
+        if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
+          connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
 
         ListOfPaths->push_back(connectedPath);
@@ -2853,5 +2931,10 @@
 
 /** 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
@@ -2861,5 +2944,5 @@
   class BoundaryLineSet *line = NULL;
   class BoundaryTriangleSet *triangle = NULL;
-  Vector OldPoint, TetraederVector[3];
+  Vector OldPoint, NormalVector;
   double volume = 0;
   int count = 0;
@@ -2887,5 +2970,5 @@
     count+=LineRunner->second->triangles.size();
   map<class BoundaryTriangleSet *, int> Candidates;
-  for (LineMap::iterator LineRunner = point->lines.begin(); (point != NULL) && (LineRunner != point->lines.end()); LineRunner++) {
+  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++) {
@@ -2897,6 +2980,8 @@
   // 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++;
@@ -2907,4 +2992,8 @@
   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];
@@ -2915,47 +3004,134 @@
 
       connectedPath = *ListRunner;
-      // initialize the path to start
-      list<TesselPoint*>::iterator CircleRunner = connectedPath->begin();
-      list<TesselPoint*>::iterator OtherCircleRunner = connectedPath->begin();
-      class TesselPoint *CentralNode = *CircleRunner;
 
       // re-create all triangles by going through connected points list
-      // advance two with CircleRunner and one with OtherCircleRunner
-      CircleRunner++;
-      CircleRunner++;
-      OtherCircleRunner++;
-      cout << Verbose(3) << "INFO: CentralNode is " << *CentralNode << "." << endl;
-      for (; (OtherCircleRunner != connectedPath->end()) && (CircleRunner != connectedPath->end()); (CircleRunner++), (OtherCircleRunner++)) {
-        cout << Verbose(4) << "INFO: CircleRunner's node is " << **CircleRunner << "." << endl;
-        cout << Verbose(4) << "INFO: OtherCircleRunner's node is " << **OtherCircleRunner << "." << endl;
-        *out << Verbose(3) << "INFO: Creating triangle " << CentralNode->Name << ", " << (*OtherCircleRunner)->Name << " and " << (*CircleRunner)->Name << "." << endl;
+      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;
+            }
+        }
+        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;
-        TriangleCandidates[0] = CentralNode;
-        TriangleCandidates[1] = *OtherCircleRunner;
-        TriangleCandidates[2] = *CircleRunner;
-        triangle = GetPresentTriangle(out, TriangleCandidates);
-        AddTesselationPoint(CentralNode, 0);
-        AddTesselationPoint(*OtherCircleRunner, 1);
-        AddTesselationPoint(*CircleRunner, 2);
+        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
-        for (int j=0;j<3;j++) {
-          TetraederVector[j].CopyVector(TPS[j]->node->node);
-          TetraederVector[j].SubtractVector(&OldPoint);
+        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();
         }
-        OldPoint.CopyVector(&TetraederVector[0]);
-        OldPoint.VectorProduct(&TetraederVector[1]);
-        volume += 1./6. * fabs(OldPoint.ScalarProduct(&TetraederVector[2]));
-        // advance number
-        NumberRunner++;
-        count++;
-        if (NumberRunner == Candidates.end())
-          *out << Verbose(2) << "WARN: Maximum of numbers reached!" << 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;
+            }
+          }
+          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.);
+      }
+
       ListOfClosedPaths->remove(connectedPath);
       delete(connectedPath);
@@ -2971,204 +3147,12 @@
     *out << Verbose(1) << "No need to create any triangles." << endl;
   }
-
   delete(ListOfClosedPaths);
 
+  *out << Verbose(1) << "Removed volume is " << volume << "." << endl;
+
   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 (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;
-            }
-          }
-        } 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(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;
-};
 
 /**
@@ -3207,18 +3191,5 @@
               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);
                 }
@@ -3238,4 +3209,41 @@
 
 /**
+ * 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.
  *
@@ -3243,44 +3251,36 @@
  *         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;
 
@@ -3294,11 +3294,22 @@
 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;
+
+  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;
@@ -3312,9 +3323,42 @@
       && (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
@@ -3323,31 +3367,157 @@
     }
   }
+  delete(DegeneratedTriangles);
+
+  cout << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << 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(4) << "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++;
+};
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision 91e7e4a510c48a116682e9998422b928cf2f1117)
+++ src/tesselation.hpp	(revision 57066a6de6f2d289266e1e6e8466859414ebf216)
@@ -34,4 +34,11 @@
 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,4 +122,5 @@
 
     void GetNormalVector(Vector &NormalVector);
+    void GetCenter(Vector *center);
     bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
     bool ContainsBoundaryLine(class BoundaryLineSet *line);
@@ -120,5 +128,5 @@
     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];
@@ -200,12 +208,11 @@
     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]);
@@ -216,9 +223,11 @@
     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();
+    void AddBoundaryPointByDegeneratedTriangle(ofstream *out, class TesselPoint *point, LinkedCell *LC);
 
     set<TesselPoint*> * GetAllConnectedPoints(ofstream *out, TesselPoint* Point);
@@ -239,4 +248,6 @@
     void PrintAllBoundaryTriangles(ofstream *out);
 
+    // store envelope in file
+    void Output(ofstream *out, const char *filename, PointCloud *cloud);
 
     PointMap PointsOnBoundary;
@@ -261,4 +272,6 @@
     class BoundaryLineSet *BLS[3];
     class BoundaryTriangleSet *BTS;
+    class BoundaryTriangleSet *LastTriangle;
+    int TriangleFilesWritten;
 
   private:
@@ -268,10 +281,4 @@
 };
 
-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_ */
Index: src/tesselationhelpers.cpp
===================================================================
--- src/tesselationhelpers.cpp	(revision 91e7e4a510c48a116682e9998422b928cf2f1117)
+++ src/tesselationhelpers.cpp	(revision 57066a6de6f2d289266e1e6e8466859414ebf216)
@@ -388,4 +388,31 @@
 };
 
+/** 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.
@@ -412,2 +439,391 @@
 };
 
+
+/** 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(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;
+};
+
+/** 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;
+  }
+};
Index: src/tesselationhelpers.hpp
===================================================================
--- src/tesselationhelpers.hpp	(revision 91e7e4a510c48a116682e9998422b928cf2f1117)
+++ src/tesselationhelpers.hpp	(revision 57066a6de6f2d289266e1e6e8466859414ebf216)
@@ -25,4 +25,5 @@
 
 #include "defs.hpp"
+#include "tesselation.hpp"
 #include "vector.hpp"
 
@@ -37,4 +38,16 @@
 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);
 
 
