Changes in src/boundary.cpp [f66195:29812d]

File:

• src/boundary.cpp (modified) (17 diffs)

src/boundary.cpp

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