Changes in src/boundary.cpp [b998c3:e359a8]

File:

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

src/boundary.cpp

@@ -10 +10 @@
 #include "element.hpp"
 #include "helpers.hpp"
-#include "info.hpp"
 #include "linkedcell.hpp"
 #include "log.hpp"
@@ -34 +33 @@
 double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
 {
-        Info FunctionInfo(__func__);
   // get points on boundary of NULL was given as parameter
   bool BoundaryFreeFlag = false;
@@ -55 +53 @@
   } else {
     BoundaryPoints = BoundaryPtr;
-    Log() << Verbose(0) << "Using given boundary points set." << endl;
+    Log() << Verbose(1) << "Using given boundary points set." << endl;
   }
   // determine biggest "diameter" of cluster for each axis
@@ -69 +67 @@
           //Log() << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
           for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-              //Log() << Verbose(1) << "Current runner is " << *(runner->second.second) << "." << endl;
+              //Log() << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
               // seek for the neighbours pair where the Othercomponent sign flips
               Neighbour = runner;
@@ -84 +82 @@
                   DistanceVector.CopyVector(&runner->second.second->x);
                   DistanceVector.SubtractVector(&Neighbour->second.second->x);
-                  //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
+                  //Log() << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
                   OldComponent) - DistanceVector.x[Othercomponent] / fabs(
@@ -93 +91 @@
                     OtherNeighbour = BoundaryPoints[axis].end();
                   OtherNeighbour--;
-                  //Log() << Verbose(1) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
+                  //Log() << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
                   // now we have found the pair: Neighbour and OtherNeighbour
                   OtherVector.CopyVector(&runner->second.second->x);
                   OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
-                  //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
-                  //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
+                  //Log() << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
+                  //Log() << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
                   // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
                   w1 = fabs(OtherVector.x[Othercomponent]);
@@ -105 +103 @@
                       * OtherVector.x[component]) / (w1 + w2));
                   // mark if it has greater diameter
-                  //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
+                  //Log() << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
                   GreatestDiameter[component] = (GreatestDiameter[component]
                       > tmp) ? GreatestDiameter[component] : tmp;
                 } //else
-              //Log() << Verbose(1) << "Saw no sign flip, probably top or bottom node." << endl;
+              //Log() << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
             }
         }
@@ -137 +135 @@
 Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
 {
-        Info FunctionInfo(__func__);
   atom *Walker = NULL;
   PointMap PointsOnBoundary;
@@ -152 +149 @@
   double angle = 0.;
 
+  Log() << Verbose(1) << "Finding all boundary points." << endl;
   // 3a. Go through every axis
   for (int axis = 0; axis < NDIM; axis++) {
@@ -178 +176 @@
         angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
 
-      //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+      //Log() << Verbose(2) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
       if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
         angle = 2. * M_PI - angle;
       }
-      Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
+      Log() << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
       BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
       if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
@@ -212 +210 @@
     // printing all inserted for debugging
     //    {
-    //      Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+    //      Log() << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     //      int i=0;
     //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
     //        if (runner != BoundaryPoints[axis].begin())
-    //          Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
+    //          Log() << Verbose(2) << ", " << i << ": " << *runner->second.second;
     //        else
-    //          Log() << Verbose(0) << i << ": " << *runner->second.second;
+    //          Log() << Verbose(2) << i << ": " << *runner->second.second;
     //        i++;
     //      }
-    //      Log() << Verbose(0) << endl;
+    //      Log() << Verbose(2) << endl;
     //    }
     // 3c. throw out points whose distance is less than the mean of left and right neighbours
@@ -251 +249 @@
           SideA.SubtractVector(MolCenter);
           SideA.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(1) << "SideA: " << SideA << endl;
+          //          Log() << Verbose(0) << "SideA: " << SideA << endl;
 
           SideB.CopyVector(&right->second.second->x);
           SideB.SubtractVector(MolCenter);
           SideB.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(1) << "SideB: " << SideB << endl;
+          //          Log() << Verbose(0) << "SideB: " << SideB << endl;
 
           SideC.CopyVector(&left->second.second->x);
           SideC.SubtractVector(&right->second.second->x);
           SideC.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(1) << "SideC: " << SideC << endl;
+          //          Log() << Verbose(0) << "SideC: " << SideC << endl;
 
           SideH.CopyVector(&runner->second.second->x);
           SideH.SubtractVector(MolCenter);
           SideH.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(1) << "SideH: " << SideH << endl;
+          //          Log() << Verbose(0) << "SideH: " << SideH << endl;
 
           // calculate each length
@@ -279 +277 @@
           const double delta = SideC.Angle(&SideH);
           const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
-          //Log() << Verbose(1) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+          //Log() << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
           Log() << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
           if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
@@ -305 +303 @@
 void FindConvexBorder(const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
 {
-        Info FunctionInfo(__func__);
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = NULL;
+
+  Log() << Verbose(1) << "Begin of FindConvexBorder" << endl;
 
   if (TesselStruct != NULL) // free if allocated
@@ -318 +317 @@
       BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
   } else {
-      Log() << Verbose(0) << "Using given boundary points set." << endl;
+      Log() << Verbose(1) << "Using given boundary points set." << endl;
   }
 
@@ -324 +323 @@
   for (int axis=0; axis < NDIM; axis++)
     {
-      Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      Log() << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
       int i=0;
       for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
         if (runner != BoundaryPoints[axis].begin())
-          Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
+          Log() << Verbose(2) << ", " << i << ": " << *runner->second.second;
         else
-          Log() << Verbose(0) << i << ": " << *runner->second.second;
+          Log() << Verbose(2) << i << ": " << *runner->second.second;
         i++;
       }
-      Log() << Verbose(0) << endl;
+      Log() << Verbose(2) << endl;
     }
 
@@ -342 +341 @@
           eLog() << Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl;
 
-  Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
+  Log() << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
   // now we have the whole set of edge points in the BoundaryList
 
@@ -348 +347 @@
   //  Log() << Verbose(1) << "Listing PointsOnBoundary:";
   //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
-  //    Log() << Verbose(0) << " " << *runner->second;
+  //    Log() << Verbose(1) << " " << *runner->second;
   //  }
-  //  Log() << Verbose(0) << endl;
+  //  Log() << Verbose(1) << endl;
 
   // 3a. guess starting triangle
@@ -360 +359 @@
   // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
   if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
-    eLog() << Verbose(1) << "Insertion of straddling points failed!" << endl;
-
-  Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
+    Log() << Verbose(1) << "Insertion of straddling points failed!" << endl;
+
+  Log() << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
 
   // 4. Store triangles in tecplot file
@@ -412 +411 @@
 //    Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
 
-  Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
+  Log() << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
 
   // 4. Store triangles in tecplot file
@@ -438 +437 @@
   if (BoundaryFreeFlag)
     delete[] (BoundaryPoints);
+
+  Log() << Verbose(1) << "End of FindConvexBorder" << endl;
 };
 
@@ -449 +450 @@
 bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
 {
-        Info FunctionInfo(__func__);
   int i=0;
   char number[MAXSTRINGSIZE];
@@ -460 +460 @@
   PointMap::iterator PointRunner;
   while (!TesselStruct->PointsOnBoundary.empty()) {
-    Log() << Verbose(1) << "Remaining points are: ";
+    Log() << Verbose(2) << "Remaining points are: ";
     for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
-      Log() << Verbose(0) << *(PointSprinter->second) << "\t";
-    Log() << Verbose(0) << endl;
+      Log() << Verbose(2) << *(PointSprinter->second) << "\t";
+    Log() << Verbose(2) << endl;
 
     PointRunner = TesselStruct->PointsOnBoundary.begin();
@@ -503 +503 @@
 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
 {
-        Info FunctionInfo(__func__);
   double volume = 0;
   class BoundaryPointSet *point = NULL;
@@ -517 +516 @@
   int run = 0;
 
+  Log() << Verbose(0) << "Begin of ConvexizeNonconvexEnvelope" << endl;
+
   // check whether there is something to work on
   if (TesselStruct == NULL) {
@@ -538 +539 @@
       for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
         line = LineRunner->second;
-        Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+        Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
         if (!line->CheckConvexityCriterion()) {
           // remove the point if needed
@@ -603 +604 @@
 
   // end
-  Log() << Verbose(0) << "Volume is " << volume << "." << endl;
+  Log() << Verbose(1) << "Volume is " << volume << "." << endl;
+  Log() << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
   return volume;
 };
@@ -617 +619 @@
 double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration)
 {
-        Info FunctionInfo(__func__);
   bool IsAngstroem = configuration->GetIsAngstroem();
   double volume = 0.;
@@ -624 +625 @@
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+  Log() << Verbose(1)
+      << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
+      << endl;
   for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
     { // go through every triangle, calculate volume of its pyramid with CoG as peak
@@ -638 +642 @@
       const double h = x.Norm(); // distance of CoG to triangle
       const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
-      Log() << Verbose(1) << "Area of triangle is " << setprecision(10) << G << " "
+      Log() << Verbose(2) << "Area of triangle is " << setprecision(10) << G << " "
           << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
           << h << " and the volume is " << PyramidVolume << " "
@@ -660 +664 @@
 void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
 {
-        Info FunctionInfo(__func__);
   // 4. Store triangles in tecplot file
   if (filename != NULL) {
@@ -668 +671 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(tecplot, TesselStruct, mol, -1);
+      WriteTecplotFile(tecplot, TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -695 +698 @@
 void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
 {
-        Info FunctionInfo(__func__);
   bool IsAngstroem = true;
   double *GreatestDiameter = NULL;
@@ -796 +798 @@
 molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
 {
-        Info FunctionInfo(__func__);
   molecule *Filling = new molecule(filler->elemente);
   Vector CurrentPosition;
@@ -813 +814 @@
   double phi[NDIM];
   class Tesselation *TesselStruct[List->ListOfMolecules.size()];
+
+  Log() << Verbose(0) << "Begin of FillBoxWithMolecule" << endl;
 
   i=0;
@@ -874 +877 @@
           for (int i=0;i<NDIM;i++)
             FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
-          Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
+          Log() << Verbose(3) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
 
           // go through all atoms
@@ -935 +938 @@
         delete(TesselStruct[i]);
   }
+  Log() << Verbose(0) << "End of FillBoxWithMolecule" << endl;
+
   return Filling;
 };
@@ -946 +951 @@
  * \param RADIUS radius of the virtual sphere
  * \param *filename filename prefix for output of vertex data
- * \return true - tesselation successful, false - tesselation failed
- */
-bool FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
+ */
+void FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
 {
-        Info FunctionInfo(__func__);
   bool freeLC = false;
-  bool status = false;
-  CandidateForTesselation *baseline;
+  LineMap::iterator baseline;
   LineMap::iterator testline;
-  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
+  bool OneLoopWithoutSuccessFlag = false;  // marks whether we went once through all baselines without finding any without two triangles
   bool TesselationFailFlag = false;
-  BoundaryTriangleSet *T = NULL;
-
+
+  Log() << Verbose(1) << "Entering search for non convex hull. " << endl;
   if (TesselStruct == NULL) {
     Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
@@ -968 +970 @@
   }
 
+  Log() << Verbose(0) << "Begin of FindNonConvexBorder\n";
+
   // initialise Linked Cell
   if (LCList == NULL) {
@@ -978 +982 @@
 
   // 2. expand from there
-  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
-    // 2a. fill all new OpenLines
-    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
-    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
-      Log() << Verbose(2) << *(Runner->second) << endl;
-
-    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
-      baseline = Runner->second;
-      if (baseline->pointlist.empty()) {
-        T = (((baseline->BaseLine->triangles.begin()))->second);
-        Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
-        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
+  baseline = TesselStruct->LinesOnBoundary.begin();
+  baseline++; // skip first line
+  while ((baseline != TesselStruct->LinesOnBoundary.end()) || (OneLoopWithoutSuccessFlag)) {
+    if (baseline->second->triangles.size() == 1) {
+      CheckListOfBaselines(TesselStruct);
+      // 3. find next triangle
+      TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*(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)
+        eLog() << Verbose(2) << "FindNextSuitableTriangle failed." << endl;
+
+      // write temporary envelope
+      if (filename != NULL) {
+        if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
+          TesselStruct->Output(filename, mol);
+        }
       }
-    }
-
-    // 2b. search for smallest ShortestAngle among all candidates
-    double ShortestAngle = 4.*M_PI;
-    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
-    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
-      Log() << Verbose(2) << *(Runner->second) << endl;
-
-    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
-      if (Runner->second->ShortestAngle < ShortestAngle) {
-        baseline = Runner->second;
-        ShortestAngle = baseline->ShortestAngle;
-        //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
+      if (TesselationFailFlag) {
+        baseline = TesselStruct->LinesOnBoundary.begin();
+        OneLoopWithoutSuccessFlag = false;
+        Log() << Verbose(2) << "Baseline set to begin." << endl;
       }
-    }
-    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+    } else {
+      //Log() << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->triangles.size() << " triangles adjacent" << endl;
+      if (baseline->second->triangles.size() != 2) {
+        eLog() << Verbose(0) << "TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
+        performCriticalExit();
+      }
+    }
+
+    if ((baseline == TesselStruct->LinesOnBoundary.end()) && (OneLoopWithoutSuccessFlag)) {
+      baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
       OneLoopWithoutSuccessFlag = false;
-    else {
-      TesselStruct->AddCandidateTriangle(*baseline);
-    }
-
-    // write temporary envelope
-    if (filename != NULL) {
-      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
-        TesselStruct->Output(filename, mol);
-      }
-    }
-  }
-//  // check envelope for consistency
-//  status = CheckListOfBaselines(TesselStruct);
-//
-//  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
-//  //->InsertStraddlingPoints(mol, LCList);
+    }
+    baseline++;
+  }
+  // check envelope for consistency
+  CheckListOfBaselines(TesselStruct);
+
+  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
+  //->InsertStraddlingPoints(mol, LCList);
 //  mol->GoToFirst();
 //  class TesselPoint *Runner = NULL;
@@ -1038 +1037 @@
 //  }
 
-//  // Purges surplus triangles.
-//  TesselStruct->RemoveDegeneratedTriangles();
+  // Purges surplus triangles.
+  TesselStruct->RemoveDegeneratedTriangles();
 
   // check envelope for consistency
-  status = CheckListOfBaselines(TesselStruct);
-
-  // store before correction
-  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
-
-  // correct degenerated polygons
-  TesselStruct->CorrectAllDegeneratedPolygons();
-
-  // check envelope for consistency
-  status = CheckListOfBaselines(TesselStruct);
+  CheckListOfBaselines(TesselStruct);
 
   // write final envelope
@@ -1059 +1049 @@
   if (freeLC)
     delete(LCList);
-
-  return status;
+  Log() << Verbose(0) << "End of FindNonConvexBorder\n";
 };
 
@@ -1072 +1061 @@
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
 {
-        Info FunctionInfo(__func__);
   Vector *Center = new Vector;
   Center->Zero();