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

File:

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

src/boundary.cpp

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