Changes in src/boundary.cpp [44fd95:196a5a]

File:

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

src/boundary.cpp

@@ -2 +2 @@
 #include "boundary.hpp"
 
-#define DEBUG 1
+#define DEBUG 0
 
 // ======================================== Points on Boundary =================================
@@ -1422 +1422 @@
  * @param c vector third point of triangle
  * @param Direction vector indicates up/down
+ * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
  * @param Halfplaneindicator double indicates whether Direction is up or down
+ * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
  * @param alpha double angle at a
  * @param beta double, angle at b
@@ -1430 +1432 @@
  */
 
-  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector* Direction, double HalfplaneIndicator,
-      double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector* Direction, Vector* AlternativeDirection,
+      double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
   {
     Vector TempNormal, helper;
     double Restradius;
-    cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
+
     *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
-    Center->Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+    Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
     // Here we calculated center of circumscribing circle, using barycentric coordinates
-    cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
 
     TempNormal.CopyVector(&a);
@@ -1446 +1447 @@
     helper.SubtractVector(&c);
     TempNormal.VectorProduct(&helper);
-    if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+    if (fabs(HalfplaneIndicator) < MYEPSILON)
       {
-        TempNormal.Scale(-1);
+        if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+          {
+            TempNormal.Scale(-1);
+          }
+      }
+    else
+      {
+        if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+          {
+            TempNormal.Scale(-1);
+          }
       }
 
     TempNormal.Normalize();
     Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-    cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
     TempNormal.Scale(Restradius);
-    cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
 
     Center->AddVector(&TempNormal);
-    cout << Verbose(4) << "Center of sphere of circumference is " << *Center << ".\n";
-    cout << Verbose(3) << "End of Get_center_of_sphere.\n";
   }
   ;
@@ -1474 +1481 @@
  * @param a first point
  * @param b second point
+ * *param c atom old third point of old triangle
  * @param Candidate base point along whose bonds to start looking from
  * @param Parent point to avoid during search as its wrong direction
@@ -1487 +1495 @@
  */
 
-void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* Candidate, atom* Parent,
+void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
     atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
 {
-  cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
-  cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
-  cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
-  cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
-  cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
-  cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
+  //cout << "ReferencePoint is " << ReferencePoint.x[0] << " "<< ReferencePoint.x[1] << " "<< ReferencePoint.x[2] << " "<< endl;
   /* OldNormal is normal vector on the old triangle
    * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
@@ -1509 +1512 @@
   double CurrentEpsilon = 0.1;
   double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-  double BallAngle;
+  double BallAngle, AlternativeSign;
   atom *Walker; // variable atom point
 
 
-  if (a != Candidate and b != Candidate)
-    {
-      cout << Verbose(3) << "We have a unique candidate!" << endl;
-      dif_a.CopyVector(&(a->x));
-      dif_a.SubtractVector(&(Candidate->x));
-      dif_b.CopyVector(&(b->x));
-      dif_b.SubtractVector(&(Candidate->x));
-      AngleCheck.CopyVector(&(Candidate->x));
-      AngleCheck.SubtractVector(&(a->x));
-      AngleCheck.ProjectOntoPlane(Chord);
-
-      SideA = dif_b.Norm();
-      SideB = dif_a.Norm();
-      SideC = Chord->Norm();
-      //Chord->Scale(-1);
-
-      alpha = Chord->Angle(&dif_a);
-      beta = M_PI - Chord->Angle(&dif_b);
-      gamma = dif_a.Angle(&dif_b);
-
-      cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-
-      if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-        cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+  dif_a.CopyVector(&(a->x));
+  dif_a.SubtractVector(&(Candidate->x));
+  dif_b.CopyVector(&(b->x));
+  dif_b.SubtractVector(&(Candidate->x));
+  AngleCheck.CopyVector(&(Candidate->x));
+  AngleCheck.SubtractVector(&(a->x));
+  AngleCheck.ProjectOntoPlane(Chord);
+
+  SideA = dif_b.Norm();
+  SideB = dif_a.Norm();
+  SideC = Chord->Norm();
+  //Chord->Scale(-1);
+
+  alpha = Chord->Angle(&dif_a);
+  beta = M_PI - Chord->Angle(&dif_b);
+  gamma = dif_a.Angle(&dif_b);
+
+
+  if (a != Candidate and b != Candidate and c != Candidate)
+    {
 
       Umkreisradius = SideA / 2.0 / sin(alpha);
@@ -1545 +1544 @@
       if (Umkreisradius < RADIUS) //Checking whether ball will at least rest on points.
         {
-          cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
+          cout << Verbose(1) << "Candidate is "<< *Candidate << endl;
           sign = AngleCheck.ScalarProduct(direction1);
-          if (fabs(sign) < MYEPSILON)
-            sign = 0;
-          else
-            sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-          if (sign >= 0)
+          if (fabs(sign)<MYEPSILON)
             {
-            cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
-
-            Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, sign, alpha, beta, gamma, RADIUS, Umkreisradius);
-
-            AngleCheck.CopyVector(&ReferencePoint);
-            AngleCheck.Scale(-1);
-            //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-            AngleCheck.AddVector(&Mittelpunkt);
-            //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-            cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
-
-            BallAngle = AngleCheck.Angle(OldNormal);
-            cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
-  //          cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
-  //          cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-            sign = AngleCheck.ScalarProduct(direction1);
-            if (fabs(sign) < MYEPSILON)
-              sign = 0;
-            else
-              sign /= fabs(sign);
-
-
-              if (Storage[0]< -1.5) // first Candidate at all
+              if (AngleCheck.ScalarProduct(OldNormal)<0)
                 {
-
-                  cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
-                  Opt_Candidate = Candidate;
-                  Storage[0] = sign;
-                  Storage[1] = BallAngle;
-                  cout << " angle " << Storage[1] << " and Up/Down "
-                    << Storage[0] << endl;
+                  sign =0;
+                  AlternativeSign=1;
                 }
               else
                 {
-                  if ( Storage[1] > BallAngle)
+                  sign =0;
+                  AlternativeSign=-1;
+                }
+            }
+          else
+            {
+              sign /= fabs(sign);
+            }
+
+
+
+          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+
+          AngleCheck.CopyVector(&ReferencePoint);
+          AngleCheck.Scale(-1);
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          AngleCheck.AddVector(&Mittelpunkt);
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+
+          BallAngle = AngleCheck.Angle(OldNormal);
+
+          //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+
+          cout << Verbose(1) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+
+          if (AngleCheck.ScalarProduct(direction1) >=0)
+            {
+              if (Storage[0]< -1.5) // first Candidate at all
+                {
+
+                  cout << "Next better candidate is " << *Candidate << " with ";
+                  Opt_Candidate = Candidate;
+                  Storage[0] = sign;
+                  Storage[1] = AlternativeSign;
+                  Storage[2] = BallAngle;
+                  cout << "Angle is " << Storage[2] << ", Halbraum ist "
+                    << Storage[0] << endl;
+
+
+                }
+              else
+                {
+                  if ( Storage[2] > BallAngle)
                     {
-                      cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
+                      cout << "Next better candidate is " << *Candidate << " with ";
                       Opt_Candidate = Candidate;
                       Storage[0] = sign;
-                      Storage[1] = BallAngle;
-                      cout << " angle " << Storage[1] << " and Up/Down "
+                      Storage[1] = AlternativeSign;
+                      Storage[2] = BallAngle;
+                      cout << "Angle is " << Storage[2] << ", Halbraum ist "
                         << Storage[0] << endl;
                     }
                   else
                     {
-                      if (DEBUG)
-                        {
-                          cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
-                        }
+                      //if (DEBUG)
+                        cout << "Looses to better candidate" << endl;
+
                     }
                 }
@@ -1607 +1618 @@
             else
               {
-              if (DEBUG)
-                {
-                  cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
-                }
+                //if (DEBUG)
+                  cout << "Refused due to bad direction of ball centre." << endl;
               }
           }
         else
           {
-            if (DEBUG)
-              {
-                cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
-              }
+            //if (DEBUG)
+              cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
           }
       }
     else
       {
-        if (DEBUG)
-          {
-            cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
-          }
+        //if (DEBUG)
+          cout << "identisch mit Ursprungslinie" << endl;
+
       }
 
 
 
-  if (RecursionLevel < 7) // Five is the recursion level threshold.
+  if (RecursionLevel < 9) // Seven is the recursion level threshold.
     {
       for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
@@ -1643 +1649 @@
           else
             {
-              Find_next_suitable_point_via_Angle_of_Sphere(a, b, Walker, Candidate, RecursionLevel
+              Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel
                   + 1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS,
                   mol); //call function again
@@ -1649 +1655 @@
         }
     }
-  cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
 }
 ;
@@ -1771 +1776 @@
           d1->ProjectOntoPlane(&AngleCheckReference);
           sign = AngleCheck.ScalarProduct(d1);
-          if (fabs(sign) < MYEPSILON)
-            sign = 0;
-          else
-            sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
+          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
 
 
@@ -1932 +1934 @@
     const double& RADIUS, int N)
 {
-  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
+  cout << Verbose(1) << "Looking for next suitable triangle \n";
   Vector direction1;
   Vector helper;
@@ -1938 +1940 @@
   ofstream *tempstream = NULL;
   char filename[255];
-  double tmp;
-  //atom* Walker;
-
-  double Storage[2];
+  atom* OldThirdPoint;
+
+  double Storage[3];
   Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[1] = 9999999.; // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[2] = 9999999.;
   atom* Opt_Candidate = NULL;
   Vector Opt_Mittelpunkt;
 
-  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-
+  cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
   helper.CopyVector(&(Line.endpoints[0]->node->x));
   for (int i = 0; i < 3; i++)
@@ -1955 +1956 @@
           && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr)
         {
+          OldThirdPoint = T.endpoints[i]->node;
           helper.SubtractVector(&T.endpoints[i]->node->x);
           break;
@@ -1968 +1970 @@
       direction1.Scale(-1);
     }
-  cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
 
   Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
-  cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
-
-
-  Vector Umkreismittelpunkt, a, b, c;
-  double alpha, beta, gamma;
-  a.CopyVector(&(T.endpoints[0]->node->x));
-  b.CopyVector(&(T.endpoints[1]->node->x));
-  c.CopyVector(&(T.endpoints[2]->node->x));
-  a.SubtractVector(&(T.endpoints[1]->node->x));
-  b.SubtractVector(&(T.endpoints[2]->node->x));
-  c.SubtractVector(&(T.endpoints[0]->node->x));
-
-  alpha = a.Angle(&c) - M_PI/2.;
-  beta = b.Angle(&a);
-  gamma = c.Angle(&b) - M_PI/2.;
-  if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-    cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-
-  cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-  Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
-  Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-  cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
-  if (DEBUG)
-    cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
-  tmp = 0;
-  for (int i=0;i<NDIM;i++) {
-    helper.CopyVector(&T.endpoints[i]->node->x);
-    helper.SubtractVector(&Umkreismittelpunkt);
-    if (DEBUG)
-      cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
-    if (tmp == 0) // set on first time for comparison against next ones
-      tmp = helper.Norm();
-    if (fabs(helper.Norm() - tmp) > MYEPSILON)
-      cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
-  }
+
+
+    Vector Umkreismittelpunkt, a, b, c;
+    double alpha, beta, gamma;
+    a.CopyVector(&(T.endpoints[0]->node->x));
+    b.CopyVector(&(T.endpoints[1]->node->x));
+    c.CopyVector(&(T.endpoints[2]->node->x));
+    a.SubtractVector(&(T.endpoints[1]->node->x));
+    b.SubtractVector(&(T.endpoints[2]->node->x));
+    c.SubtractVector(&(T.endpoints[0]->node->x));
+
+    alpha = M_PI - a.Angle(&c);
+    beta = M_PI - b.Angle(&a);
+    gamma = M_PI - c.Angle(&b);
+
+    Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
+    //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+    Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+    cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+    cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
+    cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
+
 
   cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
 
-  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node,
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
       Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
       &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node,
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
       Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
       &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
 
 
-  if ((TrianglesOnBoundaryCount % 10) == 0)
-    {
-    cout << Verbose(1) << "Writing temporary non convex hull to file testEnvelope-" << TriangleFilesWritten << "d.dat.\n";
+      cout << "Letzter Winkel bei " << TrianglesOnBoundaryCount << " Winkel ist " << Storage[2] << endl;
+
+
+  if ((TrianglesOnBoundaryCount % 1) == 0)
+    {
     sprintf(filename, "testEnvelope-%d.dat", TriangleFilesWritten);
     tempstream = new ofstream(filename, ios::trunc);
@@ -2029 +2019 @@
   }
 
-  if (TrianglesOnBoundaryCount >1000 )
-    {
-      cerr << Verbose(0)
-          << "WARNING: Already more than " << TrianglesOnBoundaryCount-1 << "triangles, construction probably has crashed." << endl;
+  if (TrianglesOnBoundaryCount >189 and TrianglesOnBoundaryCount < 200 )
+    {
+      cout << Verbose(1)
+          << "No new Atom found, triangle construction will crash" << endl;
       write_tecplot_file(out, tecplot, this, mol, TriangleFilesWritten);
-      cout << Verbose(2) << "This is currently added candidate" << Opt_Candidate << endl;
+      cout << "This is currently added candidate" << Opt_Candidate << endl;
     }
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
 
-  cout << Verbose(2) << " Optimal candidate is " << *Opt_Candidate << endl;
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
 
   AddTrianglePoint(Opt_Candidate, 0);
@@ -2050 +2040 @@
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   AddTriangleToLines();
-
-  BTS->GetNormalVector(T.NormalVector);
-
-//  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
-//      (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0))
-//    {
-//      BTS->NormalVector.Scale(-1);
-//    };
-  cout << Verbose(2) << "New triangle with " << *BTS << "and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-  cout << Verbose(2) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
-  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
+  cout << "New triangle with " << *BTS << endl;
+  cout << "We have "<< TrianglesOnBoundaryCount << endl;
+  cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
+
+  BTS->GetNormalVector(BTS->NormalVector);
+
+  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
+      (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) ||
+      (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) )
+
+    {
+      BTS->NormalVector.Scale(-1);
+    };
+
 }
 ;
 
 void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
-    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[2],
+    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[3],
     molecule* mol, double RADIUS)
 {
-  cout << Verbose(2)
-      << "Begin of Find_second_point_for_Tesselation, recursive level "
-      << RecursionLevel << endl;
+  cout << Verbose(1)
+      << "Looking for second point of starting triangle, recursive level "
+      << RecursionLevel << endl;;
   int i;
   Vector AngleCheck;
   atom* Walker;
-  double norm = -1., angle;
+  double norm = -1.;
 
   // check if we only have one unique point yet ...
   if (a != Candidate)
     {
-    cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
       AngleCheck.CopyVector(&(Candidate->x));
       AngleCheck.SubtractVector(&(a->x));
@@ -2086 +2078 @@
       if (norm < RADIUS)
         {
-        angle = AngleCheck.Angle(&Oben);
-          if (angle < Storage[0])
+          if (AngleCheck.Angle(&Oben) < Storage[0])
             {
-              //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-              cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
+              //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[2]);
+              cout << "Next better candidate is " << *Candidate
+                  << " with distance " << norm << ".\n";
               Opt_Candidate = Candidate;
               Storage[0] = AngleCheck.Angle(&Oben);
-              //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+              //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
             }
           else
             {
-              cout << "Looses with angle " << angle << " to a better candidate "
+              cout << Verbose(1) << "Supposedly looses to a better candidate "
                   << *Opt_Candidate << endl;
             }
@@ -2103 +2095 @@
       else
         {
-          cout << "Refused due to Radius " << norm
+          cout << Verbose(1) << *Candidate << " refused due to Radius " << norm
               << endl;
         }
@@ -2122 +2114 @@
         };
     };
-  cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level "
-      << RecursionLevel << endl;
 }
 ;
@@ -2129 +2119 @@
 void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
-  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+  cout << Verbose(1) << "Looking for starting triangle \n";
   int i = 0;
   atom* Walker;
   atom* FirstPoint;
   atom* SecondPoint;
-  atom* max_index[NDIM];
-  double max_coordinate[NDIM];
+  atom* max_index[3];
+  double max_coordinate[3];
   Vector Oben;
   Vector helper;
@@ -2148 +2138 @@
       max_coordinate[i] = -1;
     }
-  cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount
+  cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount
       << " Atoms \n";
 
@@ -2166 +2156 @@
     }
 
-  cout << Verbose(2) << "Found maximum coordinates: ";
-  for (int i=0;i<NDIM;i++)
-    cout << i << ": " << *max_index[i] << "\t";
-  cout << endl;
+  cout << Verbose(1) << "Found maximum coordinates. " << endl;
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
   const int k = 1;
@@ -2175 +2162 @@
   FirstPoint = max_index[k];
 
-  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
-  double Storage[2];
+  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << ": "
+      << FirstPoint->x.x[0] << endl;
+  double Storage[3];
   atom* Opt_Candidate = NULL;
   Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
   Storage[1] = 999999.; // This will be an angle looking for the third point.
+  Storage[2] = 999999.;
+  cout << Verbose(1) << "Number of Bonds: "
+      << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
 
   Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0,
       Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
   SecondPoint = Opt_Candidate;
-  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+  cout << Verbose(1) << "Found second point is " << *SecondPoint << ".\n";
 
   helper.CopyVector(&(FirstPoint->x));
@@ -2194 +2185 @@
   Storage[0] = -2.; // This will indicate the quadrant.
   Storage[1] = 9999999.; // This will be an angle looking for the third point.
+  Storage[2] = 9999999.;
 
   Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
@@ -2199 +2191 @@
   // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
 
-  cout << Verbose(2) << "Looking for third point candidates \n";
+  cout << Verbose(1) << "Looking for third point candidates \n";
+  cout << Verbose(1) << " In direction " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << " " << endl;
   // look in one direction of baseline for initial candidate
   Opt_Candidate = NULL;
   CenterOfFirstLine.CopyVector(&Chord);
-  CenterOfFirstLine.Scale(-0.5);
+  CenterOfFirstLine.Scale(0.5);
   CenterOfFirstLine.AddVector(&(SecondPoint->x));
 
-  cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
-  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, 0,
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0,
       &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
   // look in other direction of baseline for possible better candidate
-  cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
-  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, FirstPoint, SecondPoint, 0,
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0,
       &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
   cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
 
   // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+
+  cout << Verbose(1) << "The found starting triangle consists of "
+      << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
+      << "." << endl;
 
   // Finally, we only have to add the found points
@@ -2231 +2226 @@
   Oben.Scale(-1.);
   BTS->GetNormalVector(Oben);
-  cout << Verbose(0) << "==> The found starting triangle consists of "
-      << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
-      << " with normal vector " << BTS->NormalVector << ".\n";
-  cout << Verbose(1) << "End of Find_starting_triangle\n";
 }
 ;
@@ -2246 +2237 @@
   const double RADIUS = 6.;
   LineMap::iterator baseline;
-  cout << Verbose(0) << "Begin of Find_non_convex_border\n";
-  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
-
   Tess->Find_starting_triangle(mol, RADIUS);
 
   baseline = Tess->LinesOnBoundary.begin();
-  while ((baseline != Tess->LinesOnBoundary.end()) || (!flag))
+  while (baseline != Tess->LinesOnBoundary.end())
     {
       if (baseline->second->TrianglesCount == 1)
         {
+          cout << Verbose(1) << "Begin of Tesselation ... " << endl;
           Tess->Find_next_suitable_triangle(out, tecplot, mol,
               *(baseline->second),
               *(((baseline->second->triangles.begin()))->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
-          flag = true;
+          cout << Verbose(1) << "End of Tesselation ... " << endl;
         }
       else
         {
-          cout << Verbose(1) << "Line " << &baseline << " has "
+          cout << Verbose(1) << "There is a line with "
               << baseline->second->TrianglesCount << " triangles adjacent"
               << endl;
@@ -2269 +2258 @@
       N++;
       baseline++;
-      if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
-        baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-        flag = false;
-      }
-    }
-  cout << Verbose(1) << "Writing final tecplot file\n";
+    }
   write_tecplot_file(out, tecplot, Tess, mol, -1);
 
-  cout << Verbose(0) << "End of Find_non_convex_border\n";
-}
-;
+}
+;