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

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"
@@ -33 +34 @@
 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;
@@ -53 +55 @@
   } else {
     BoundaryPoints = BoundaryPtr;
-    Log() << Verbose(1) << "Using given boundary points set." << endl;
+    Log() << Verbose(0) << "Using given boundary points set." << endl;
   }
   // determine biggest "diameter" of cluster for each axis
@@ -67 +69 @@
           //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(2) << "Current runner is " << *(runner->second.second) << "." << endl;
+              //Log() << Verbose(1) << "Current runner is " << *(runner->second.second) << "." << endl;
               // seek for the neighbours pair where the Othercomponent sign flips
               Neighbour = runner;
@@ -82 +84 @@
                   DistanceVector.CopyVector(&runner->second.second->x);
                   DistanceVector.SubtractVector(&Neighbour->second.second->x);
-                  //Log() << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
+                  //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
                   OldComponent) - DistanceVector.x[Othercomponent] / fabs(
@@ -91 +93 @@
                     OtherNeighbour = BoundaryPoints[axis].end();
                   OtherNeighbour--;
-                  //Log() << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
+                  //Log() << Verbose(1) << "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(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;
+                  //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;
                   // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
                   w1 = fabs(OtherVector.x[Othercomponent]);
@@ -103 +105 @@
                       * OtherVector.x[component]) / (w1 + w2));
                   // mark if it has greater diameter
-                  //Log() << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
+                  //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
                   GreatestDiameter[component] = (GreatestDiameter[component]
                       > tmp) ? GreatestDiameter[component] : tmp;
                 } //else
-              //Log() << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
+              //Log() << Verbose(1) << "Saw no sign flip, probably top or bottom node." << endl;
             }
         }
@@ -135 +137 @@
 Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
 {
+        Info FunctionInfo(__func__);
   atom *Walker = NULL;
   PointMap PointsOnBoundary;
@@ -149 +152 @@
   double angle = 0.;
 
-  Log() << Verbose(1) << "Finding all boundary points." << endl;
   // 3a. Go through every axis
   for (int axis = 0; axis < NDIM; axis++) {
@@ -176 +178 @@
         angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
 
-      //Log() << Verbose(2) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+      //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
       if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
         angle = 2. * M_PI - angle;
       }
-      Log() << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
+      Log() << Verbose(1) << "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
@@ -210 +212 @@
     // printing all inserted for debugging
     //    {
-    //      Log() << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+    //      Log() << Verbose(1) << "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(2) << ", " << i << ": " << *runner->second.second;
+    //          Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
     //        else
-    //          Log() << Verbose(2) << i << ": " << *runner->second.second;
+    //          Log() << Verbose(0) << i << ": " << *runner->second.second;
     //        i++;
     //      }
-    //      Log() << Verbose(2) << endl;
+    //      Log() << Verbose(0) << endl;
     //    }
     // 3c. throw out points whose distance is less than the mean of left and right neighbours
@@ -249 +251 @@
           SideA.SubtractVector(MolCenter);
           SideA.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(0) << "SideA: " << SideA << endl;
+          //          Log() << Verbose(1) << "SideA: " << SideA << endl;
 
           SideB.CopyVector(&right->second.second->x);
           SideB.SubtractVector(MolCenter);
           SideB.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(0) << "SideB: " << SideB << endl;
+          //          Log() << Verbose(1) << "SideB: " << SideB << endl;
 
           SideC.CopyVector(&left->second.second->x);
           SideC.SubtractVector(&right->second.second->x);
           SideC.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(0) << "SideC: " << SideC << endl;
+          //          Log() << Verbose(1) << "SideC: " << SideC << endl;
 
           SideH.CopyVector(&runner->second.second->x);
           SideH.SubtractVector(MolCenter);
           SideH.ProjectOntoPlane(&AxisVector);
-          //          Log() << Verbose(0) << "SideH: " << SideH << endl;
+          //          Log() << Verbose(1) << "SideH: " << SideH << endl;
 
           // calculate each length
@@ -277 +279 @@
           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(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) << " 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) {
@@ -303 +305 @@
 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
@@ -317 +318 @@
       BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
   } else {
-      Log() << Verbose(1) << "Using given boundary points set." << endl;
+      Log() << Verbose(0) << "Using given boundary points set." << endl;
   }
 
@@ -323 +324 @@
   for (int axis=0; axis < NDIM; axis++)
     {
-      Log() << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      Log() << Verbose(1) << "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(2) << ", " << i << ": " << *runner->second.second;
+          Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
         else
-          Log() << Verbose(2) << i << ": " << *runner->second.second;
+          Log() << Verbose(0) << i << ": " << *runner->second.second;
         i++;
       }
-      Log() << Verbose(2) << endl;
+      Log() << Verbose(0) << endl;
     }
 
@@ -341 +342 @@
           eLog() << Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl;
 
-  Log() << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
+  Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
   // now we have the whole set of edge points in the BoundaryList
 
@@ -347 +348 @@
   //  Log() << Verbose(1) << "Listing PointsOnBoundary:";
   //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
-  //    Log() << Verbose(1) << " " << *runner->second;
+  //    Log() << Verbose(0) << " " << *runner->second;
   //  }
-  //  Log() << Verbose(1) << endl;
+  //  Log() << Verbose(0) << endl;
 
   // 3a. guess starting triangle
@@ -359 +360 @@
   // 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))
-    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;
+    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;
 
   // 4. Store triangles in tecplot file
@@ -411 +412 @@
 //    Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
 
-  Log() << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
+  Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
 
   // 4. Store triangles in tecplot file
@@ -437 +438 @@
   if (BoundaryFreeFlag)
     delete[] (BoundaryPoints);
-
-  Log() << Verbose(1) << "End of FindConvexBorder" << endl;
 };
 
@@ -450 +449 @@
 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(2) << "Remaining points are: ";
+    Log() << Verbose(1) << "Remaining points are: ";
     for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
-      Log() << Verbose(2) << *(PointSprinter->second) << "\t";
-    Log() << Verbose(2) << endl;
+      Log() << Verbose(0) << *(PointSprinter->second) << "\t";
+    Log() << Verbose(0) << 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;
@@ -516 +517 @@
   int run = 0;
 
-  Log() << Verbose(0) << "Begin of ConvexizeNonconvexEnvelope" << endl;
-
   // check whether there is something to work on
   if (TesselStruct == NULL) {
@@ -539 +538 @@
       for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
         line = LineRunner->second;
-        Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+        Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
         if (!line->CheckConvexityCriterion()) {
           // remove the point if needed
@@ -604 +603 @@
 
   // end
-  Log() << Verbose(1) << "Volume is " << volume << "." << endl;
-  Log() << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
+  Log() << Verbose(0) << "Volume is " << volume << "." << endl;
   return volume;
 };
@@ -619 +617 @@
 double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration)
 {
+        Info FunctionInfo(__func__);
   bool IsAngstroem = configuration->GetIsAngstroem();
   double volume = 0.;
@@ -625 +624 @@
 
   // 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
@@ -642 +638 @@
       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(2) << "Area of triangle is " << setprecision(10) << G << " "
+      Log() << Verbose(1) << "Area of triangle is " << setprecision(10) << G << " "
           << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
           << h << " and the volume is " << PyramidVolume << " "
@@ -664 +660 @@
 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) {
@@ -671 +668 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(tecplot, TesselStruct, mol, 0);
+      WriteTecplotFile(tecplot, TesselStruct, mol, -1);
       tecplot->close();
       delete(tecplot);
@@ -698 +695 @@
 void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
 {
+        Info FunctionInfo(__func__);
   bool IsAngstroem = true;
   double *GreatestDiameter = NULL;
@@ -798 +796 @@
 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;
@@ -814 +813 @@
   double phi[NDIM];
   class Tesselation *TesselStruct[List->ListOfMolecules.size()];
-
-  Log() << Verbose(0) << "Begin of FillBoxWithMolecule" << endl;
 
   i=0;
@@ -877 +874 @@
           for (int i=0;i<NDIM;i++)
             FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
-          Log() << Verbose(3) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
+          Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
 
           // go through all atoms
@@ -938 +935 @@
         delete(TesselStruct[i]);
   }
-  Log() << Verbose(0) << "End of FillBoxWithMolecule" << endl;
-
   return Filling;
 };
@@ -951 +946 @@
  * \param RADIUS radius of the virtual sphere
  * \param *filename filename prefix for output of vertex data
- */
-void FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
+ * \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)
 {
+        Info FunctionInfo(__func__);
   bool freeLC = false;
-  LineMap::iterator baseline;
+  bool status = false;
+  CandidateForTesselation *baseline;
   LineMap::iterator testline;
-  bool OneLoopWithoutSuccessFlag = false;  // marks whether we went once through all baselines without finding any without two triangles
+  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
   bool TesselationFailFlag = false;
-
-  Log() << Verbose(1) << "Entering search for non convex hull. " << endl;
+  BoundaryTriangleSet *T = NULL;
+
   if (TesselStruct == NULL) {
     Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
@@ -970 +968 @@
   }
 
-  Log() << Verbose(0) << "Begin of FindNonConvexBorder\n";
-
   // initialise Linked Cell
   if (LCList == NULL) {
@@ -982 +978 @@
 
   // 2. expand from there
-  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);
-        }
+  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.
       }
-      if (TesselationFailFlag) {
-        baseline = TesselStruct->LinesOnBoundary.begin();
-        OneLoopWithoutSuccessFlag = false;
-        Log() << Verbose(2) << "Baseline set to begin." << endl;
+    }
+
+    // 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;
       }
-    } 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 ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+      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);
       }
     }
-
-    if ((baseline == TesselStruct->LinesOnBoundary.end()) && (OneLoopWithoutSuccessFlag)) {
-      baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-      OneLoopWithoutSuccessFlag = false;
-    }
-    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);
+  }
+//  // 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);
 //  mol->GoToFirst();
 //  class TesselPoint *Runner = NULL;
@@ -1037 +1038 @@
 //  }
 
-  // Purges surplus triangles.
-  TesselStruct->RemoveDegeneratedTriangles();
+//  // Purges surplus triangles.
+//  TesselStruct->RemoveDegeneratedTriangles();
 
   // check envelope for consistency
-  CheckListOfBaselines(TesselStruct);
+  status = CheckListOfBaselines(TesselStruct);
+
+  // store before correction
+  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
+
+  // correct degenerated polygons
+  TesselStruct->CorrectAllDegeneratedPolygons();
+
+  // check envelope for consistency
+  status = CheckListOfBaselines(TesselStruct);
 
   // write final envelope
@@ -1049 +1059 @@
   if (freeLC)
     delete(LCList);
-  Log() << Verbose(0) << "End of FindNonConvexBorder\n";
+
+  return status;
 };
 
@@ -1061 +1072 @@
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
 {
+        Info FunctionInfo(__func__);
   Vector *Center = new Vector;
   Center->Zero();