Changeset a8bcea6 for src/boundary.cpp

Timestamp:

Dec 4, 2008, 3:15:00 PM (17 years ago)

Author:

Christian Neuen <neuen@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

f714979

Parents:

f683fe

Message:

several changes, now output is created, quality unknown

File:

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

src/boundary.cpp

@@ -456 +456 @@
 
 
+/*
+ * 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
+ */
+void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol)
+{
+  // 8. Store triangles in tecplot file
+  if (tecplot != NULL) {
+    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
+    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
+    *tecplot << "ZONE T=\"TRIANGLES\", N=" <<  TesselStruct->PointsOnBoundaryCount << ", E=" <<  TesselStruct->TrianglesOnBoundaryCount << ", 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;
+    atom *Walker = NULL;
+    for (PointMap::iterator target =  TesselStruct->PointsOnBoundary.begin(); target !=  TesselStruct->PointsOnBoundary.end(); target++) {
+      Walker = target->second->node;
+      LookupList[Walker->nr] = Counter++;
+      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << 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 volume of a cluster.
  * Determines first the convex envelope, then tesselates it and calculates its volume.
@@ -478 +514 @@
   double a,b,c;
 
-  Find_non_convex_border(*TesselStruct, mol);
+  Find_non_convex_border(out, tecplot, *TesselStruct, mol);
   /*
   // 1. calculate center of gravity
@@ -560 +596 @@
 
 
-
-  // 8. Store triangles in tecplot file
-  if (tecplot != NULL) {
-    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES\", N=" <<  TesselStruct->PointsOnBoundaryCount << ", E=" <<  TesselStruct->TrianglesOnBoundaryCount << ", 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;
-    atom *Walker = NULL;
-    for (PointMap::iterator target =  TesselStruct->PointsOnBoundary.begin(); target !=  TesselStruct->PointsOnBoundary.end(); target++) {
-      Walker = target->second->node;
-      LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << 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;
-  }
+  write_tecplot_file(out, tecplot, TesselStruct, mol);
+
 
   // free reference lists
@@ -699 +709 @@
 Tesselation::~Tesselation()
 {
+        cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
   for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
     delete(runner->second);
@@ -1028 +1039 @@
  *  with all neighbours of the candidate.
  */
-void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule* mol)
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
 {
   /* d2 is normal vector on the triangle
@@ -1045 +1056 @@
   AngleCheck.ProjectOntoPlane(Chord);
 
-  if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+  if (problem)
   {
-
-    if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
+          cout << "Atom number" << Candidate->nr << endl;
+          Candidate->x.Output((ofstream *)&cout);
+          cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
+  }
+
+  if (a != Candidate and b != Candidate)
+  {
+
+          if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+          {
+                  if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[0]) //This will give absolute preference to those in "right-hand" quadrants
+                  {
+                          Opt_Candidate = Candidate;
+                          Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+                          Storage[1]=AngleCheck.Angle(d2);
+                  }
+                  else
+                  {
+                          if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]< AngleCheck.Angle(d2)) or \
+                                          (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]> AngleCheck.Angle(d2)))
+                                  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+                          {
+                                  Opt_Candidate = Candidate;
+                                  Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+                                  Storage[1]=AngleCheck.Angle(d2);
+                          }
+                          else
+                          {
+                                  if (problem)
+                                          cout << "Loses to better candidate" << endl;
+                          }
+                  }
+          }
+          else
+          {
+                  if (problem)
+                          cout << "erfuellt Dreiecksbedingung fuer sehne nicht" <<endl;
+          }
+  }
+  else
+  {
+          if (problem)
+                          cout << "identisch mit Ursprungslinie" << endl;
+  }
+
+  if (n<5) // Five is the recursion level threshold.
+  {
+          for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
       {
-        Storage[0]=(double)Candidate->nr;
-        Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-        Storage[2]=AngleCheck.Angle(d2);
+                  Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+
+                  Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
       }
-    else
-      {
-        if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]>0 &&  Storage[2]< AngleCheck.Angle(d2)) or \
-            (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]<0 &&  Storage[2]> AngleCheck.Angle(d2)))
-          //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-          {
-            Storage[0]=(double)Candidate->nr;
-            Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-            Storage[2]=AngleCheck.Angle(d2);
-          }
-      }
-  }
-
-  if (n<5) // Five is the recursion level threshold.
-    {
-      for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
-        {
-                Walker= mol->start; // go through all atoms
-
-          while (Walker->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr))
-            { // until atom found which belongs to bond
-              Walker = Walker->next;
-            }
-
-
-          Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Storage, RADIUS, mol);  //call function again
-        }
-    }
+  }
 };
 
@@ -1087 +1117 @@
  * this function fins a triangle to a line, adjacent to an existing one.
  */
-void Tesselation::Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
-{
-        printf("Looking for next suitable triangle \n");
+void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS)
+{
+        cout << Verbose(1) << "Looking for next suitable triangle \n";
   Vector direction1;
   Vector helper;
@@ -1095 +1125 @@
   atom* Walker;
 
-  double *Storage;
-  Storage = new double[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
-
-
+  double Storage[2];
+  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
+  atom* Opt_Candidate = NULL;
+
+
+  cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
   helper.CopyVector(&(Line.endpoints[0]->node->x));
   for (int i =0; i<3; i++)
@@ -1115 +1145 @@
   direction1.CopyVector(&Line.endpoints[0]->node->x);
   direction1.SubtractVector(&Line.endpoints[1]->node->x);
-  direction1.VectorProduct(T.NormalVector);
+  direction1.VectorProduct(&(T.NormalVector));
 
   if (direction1.ScalarProduct(&helper)>0)
@@ -1125 +1155 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
-printf("Looking for third point candidates for triangle \n");
-  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, T.NormalVector, Storage, RADIUS, mol);
-
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol,0);
+
+  if (Opt_Candidate==NULL)
+  {
+          cout << Verbose(1) << "No new Atom found, triangle construction will crash" << endl;
+          write_tecplot_file(out, tecplot, this, mol);
+          tecplot->flush();
+          Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol, 1);
+
+  }
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
-  // Next Triangle is Line, atom with number in Storage[0]
-
-  Walker= mol->start;
-  while (Walker->nr != (int)Storage[0])
-    {
-      Walker = Walker->next;
-    }
-
-  AddPoint(Walker);
-
-  BPS[0] = new class BoundaryPointSet(Walker);
+
+  cout << Verbose(1) << "Adding exactly one Walker for reasons completely unknown to me ... " << endl;
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+  AddPoint(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+
+  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(Walker);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
   BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
@@ -1151 +1189 @@
   TrianglesOnBoundaryCount++;
 
+  cout << Verbose(1) << "Checking for already present lines ... " << endl;
   for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
     {
@@ -1160 +1199 @@
         }
     }
-  BTS->GetNormalVector(*BTS->NormalVector);
-
-  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
-      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
+  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)) )
     {
-      BTS->NormalVector->Scale(-1);
-    }
-
-};
-
-
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, double Storage[3], molecule* mol)
-{
-        printf("Looking for second point of starting triangle \n");
-  int i;
-  Vector *AngleCheck;
-  atom* Walker;
-
-  AngleCheck->CopyVector(&(Candidate->x));
-  AngleCheck->SubtractVector(&(a->x));
-  if (AngleCheck->Angle(&Oben) < Storage[1])
+      BTS->NormalVector.Scale(-1);
+    };
+
+};
+
+
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, atom*& Opt_Candidate, double Storage[2], molecule* mol, double RADIUS)
+{
+        cout << Verbose(1) << "Looking for second point of starting triangle, recursive level "<< n <<endl;;
+        int i;
+        Vector AngleCheck;
+        atom* Walker;
+
+        AngleCheck.CopyVector(&(Candidate->x));
+        AngleCheck.SubtractVector(&(a->x));
+        if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
     {
-          //printf("Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-      Storage[0]=(double)(Candidate->nr);
-      Storage[1]=AngleCheck->Angle(&Oben);
-      //printf("Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+          //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+      Opt_Candidate=Candidate;
+      Storage[0]=AngleCheck.Angle(&Oben);
+      //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
     };
-  printf("%d \n", n);
-  if (n<5)
+        if (n<5)
     {
       for (i = 0; i< mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
-        {
-          Walker = mol->start;
-          while (Candidate->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr))
-            {
-              Walker = Walker->next;
-            };
-
-          Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Storage, mol);
-            };
+      {
+          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+          Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Opt_Candidate, Storage, mol, RADIUS);
+      };
     };
 
@@ -1208 +1243 @@
 void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
-        printf("Looking for starting triangle \n");
+        cout << Verbose(1) << "Looking for starting triangle \n";
   int i=0;
   atom* Walker;
@@ -1227 +1262 @@
       max_coordinate[i] =-1;
     }
-printf("Molecule mol is there and has %d Atoms \n", mol->AtomCount);
+cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
   Walker = mol->start;
   while (Walker->next != mol->end)
@@ -1242 +1277 @@
     }
 
-printf("Found starting atom \n");
+cout << Verbose(1) << "Found starting atom \n";
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
   const int k=2;
@@ -1253 +1288 @@
       Walker = Walker->next;
     }
-printf("%d \n", Walker->nr);
-  double Storage[3];
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=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[2]=999999.;  // This will be an angle looking for the third point.
-printf("%d \n", mol->NumberOfBondsPerAtom[Walker->nr]);
+cout << Verbose(1) << "Number of start atom: " << Walker->nr << endl;
+  double Storage[2];
+  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.
+cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[Walker->nr] << endl;
 
   for (i=1; i< (mol->NumberOfBondsPerAtom[Walker->nr]); i++)
     {
-      Walker2 = mol->start;
-      Walker2 = Walker2->next;
-      while (Walker2->nr != (mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr) )
-        {
-          Walker2 = Walker2->next;
-        }
-
-      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
-    }
-
-  Walker2 = mol->start;
-
-  while (Walker2->nr != int(Storage[0]))
-    {
-      Walker2 = Walker2->next;
-    }
-
+          Walker2 = mol->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
+    }
+
+  Walker2 = Opt_Candidate;
+  Opt_Candidate=NULL;
   helper.CopyVector(&(Walker->x));
   helper.SubtractVector(&(Walker2->x));
   Oben.ProjectOntoPlane(&helper);
   helper.VectorProduct(&Oben);
-  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-2.;   // 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[2]= -10.;  // This will be an angle looking for the third point.
+  Storage[0]=-2.;       // This will indicate the quadrant.
+  Storage[1]= 9999999.; // This will be an angle looking for the third point.
 
   Chord.CopyVector(&(Walker->x)); // bring into calling function
   Chord.SubtractVector(&(Walker2->x));
 
-printf("Looking for third point candidates \n");
-       Find_next_suitable_point(Walker, Walker2, (mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][0]->rightatom : mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom), 0, &Chord, &helper, &Oben, Storage, RADIUS, mol);
-  Walker3 = mol->start;
-  while (Walker3->nr != int(Storage[0]))
-    {
-      Walker3 = Walker3->next;
-    }
-
-  //Starting Triangle is Walker, Walker2, index Storage[0]
+  cout << Verbose(1) << "Looking for third point candidates \n";
+  Find_next_suitable_point(Walker, Walker2, mol->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker), 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
+
+
+  //Starting Triangle is Walker, Walker2, Opt_Candidate
 
   AddPoint(Walker);
   AddPoint(Walker2);
-  AddPoint(Walker3);
+  AddPoint(Opt_Candidate);
+
+  cout << Verbose(1) << "The found starting triangle consists of " << *Walker << ", " << *Walker2 << " and " << *Opt_Candidate << "." << endl;
 
   BPS[0] = new class BoundaryPointSet(Walker);
@@ -1308 +1329 @@
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
-  BPS[1] = new class BoundaryPointSet(Walker3);
+  BPS[1] = new class BoundaryPointSet(Opt_Candidate);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
   BPS[0] = new class BoundaryPointSet(Walker);
@@ -1314 +1335 @@
   BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
 
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  Tesselation::BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
   TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
   TrianglesOnBoundaryCount++;
@@ -1324 +1345 @@
     };
 
-       BTS->GetNormalVector(*BTS->NormalVector);
-
-       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
+       BTS->GetNormalVector(BTS->NormalVector);
+
+       if( BTS->NormalVector.ScalarProduct(&Oben)<0)
          {
-           BTS->NormalVector->Scale(-1);
+           BTS->NormalVector.Scale(-1);
          }
 };
 
 
-void Find_non_convex_border(Tesselation Tess, molecule* mol)
-{
-        printf("Entering finding of non convex hull. \n");
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol)
+{
+        cout << Verbose(1) << "Entering finding of non convex hull. " << endl;
+        cout << flush;
   const double RADIUS =6;
+  LineMap::iterator baseline;
   Tess.Find_starting_triangle(mol, RADIUS);
 
-  for (LineMap::iterator baseline = Tess.LinesOnBoundary.begin(); baseline != Tess.LinesOnBoundary.end(); baseline++)
+  baseline = Tess.LinesOnBoundary.begin();
+  while (baseline != Tess.LinesOnBoundary.end()) {
     if (baseline->second->TrianglesCount == 1)
       {
-                Tess.Find_next_suitable_triangle(mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
+        cout << Verbose(1) << "Begin of Tesselation ... " << endl;
+                Tess.Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
+        cout << Verbose(1) << "End of Tesselation ... " << endl;
       }
     else
       {
-        printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
+        cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
       }
-};
+  baseline++;
+  }
+
+};