Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    r7f4bee r4fc93f  
    1010#include "element.hpp"
    1111#include "helpers.hpp"
     12#include "info.hpp"
    1213#include "linkedcell.hpp"
    1314#include "log.hpp"
    … …  
    3334double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
    3435{
     36        Info FunctionInfo(__func__);
    3537  // get points on boundary of NULL was given as parameter
    3638  bool BoundaryFreeFlag = false;
    … …  
    5355  } else {
    5456    BoundaryPoints = BoundaryPtr;
    55     Log() << Verbose(1) << "Using given boundary points set." << endl;
     57    Log() << Verbose(0) << "Using given boundary points set." << endl;
    5658  }
    5759  // determine biggest "diameter" of cluster for each axis
    … …  
    6769          //Log() << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
    6870          for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    69               //Log() << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
     71              //Log() << Verbose(1) << "Current runner is " << *(runner->second.second) << "." << endl;
    7072              // seek for the neighbours pair where the Othercomponent sign flips
    7173              Neighbour = runner;
    … …  
    8284                  DistanceVector.CopyVector(&runner->second.second->x);
    8385                  DistanceVector.SubtractVector(&Neighbour->second.second->x);
    84                   //Log() << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
     86                  //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
    8587                } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
    8688                  OldComponent) - DistanceVector.x[Othercomponent] / fabs(
    … …  
    9193                    OtherNeighbour = BoundaryPoints[axis].end();
    9294                  OtherNeighbour--;
    93                   //Log() << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
     95                  //Log() << Verbose(1) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
    9496                  // now we have found the pair: Neighbour and OtherNeighbour
    9597                  OtherVector.CopyVector(&runner->second.second->x);
    9698                  OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
    97                   //Log() << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
    98                   //Log() << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
     99                  //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
     100                  //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
    99101                  // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
    100102                  w1 = fabs(OtherVector.x[Othercomponent]);
    … …  
    103105                      * OtherVector.x[component]) / (w1 + w2));
    104106                  // mark if it has greater diameter
    105                   //Log() << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
     107                  //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
    106108                  GreatestDiameter[component] = (GreatestDiameter[component]
    107109                      > tmp) ? GreatestDiameter[component] : tmp;
    108110                } //else
    109               //Log() << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
     111              //Log() << Verbose(1) << "Saw no sign flip, probably top or bottom node." << endl;
    110112            }
    111113        }
    … …  
    135137Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
    136138{
     139        Info FunctionInfo(__func__);
    137140  atom *Walker = NULL;
    138141  PointMap PointsOnBoundary;
    … …  
    149152  double angle = 0.;
    150153
    151   Log() << Verbose(1) << "Finding all boundary points." << endl;
    152154  // 3a. Go through every axis
    153155  for (int axis = 0; axis < NDIM; axis++) {
    … …  
    176178        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
    177179
    178       //Log() << Verbose(0) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
     180      //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
    179181      if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
    180182        angle = 2. * M_PI - angle;
    181183      }
    182       Log() << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
     184      Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
    183185      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
    184186      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
    … …  
    210212    // printing all inserted for debugging
    211213    //    {
    212     //      Log() << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     214    //      Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
    213215    //      int i=0;
    214216    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    … …  
    249251          SideA.SubtractVector(MolCenter);
    250252          SideA.ProjectOntoPlane(&AxisVector);
    251           //          Log() << Verbose(0) << "SideA: ";
    252           //          SideA.Output(out);
    253           //          Log() << Verbose(0) << endl;
     253          //          Log() << Verbose(1) << "SideA: " << SideA << endl;
    254254
    255255          SideB.CopyVector(&right->second.second->x);
    256256          SideB.SubtractVector(MolCenter);
    257257          SideB.ProjectOntoPlane(&AxisVector);
    258           //          Log() << Verbose(0) << "SideB: ";
    259           //          SideB.Output(out);
    260           //          Log() << Verbose(0) << endl;
     258          //          Log() << Verbose(1) << "SideB: " << SideB << endl;
    261259
    262260          SideC.CopyVector(&left->second.second->x);
    263261          SideC.SubtractVector(&right->second.second->x);
    264262          SideC.ProjectOntoPlane(&AxisVector);
    265           //          Log() << Verbose(0) << "SideC: ";
    266           //          SideC.Output(out);
    267           //          Log() << Verbose(0) << endl;
     263          //          Log() << Verbose(1) << "SideC: " << SideC << endl;
    268264
    269265          SideH.CopyVector(&runner->second.second->x);
    270266          SideH.SubtractVector(MolCenter);
    271267          SideH.ProjectOntoPlane(&AxisVector);
    272           //          Log() << Verbose(0) << "SideH: ";
    273           //          SideH.Output(out);
    274           //          Log() << Verbose(0) << endl;
     268          //          Log() << Verbose(1) << "SideH: " << SideH << endl;
    275269
    276270          // calculate each length
    … …  
    285279          const double delta = SideC.Angle(&SideH);
    286280          const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
    287           //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;
     281          //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;
    288282          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;
    289283          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
    … …  
    311305void FindConvexBorder(const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
    312306{
     307        Info FunctionInfo(__func__);
    313308  bool BoundaryFreeFlag = false;
    314309  Boundaries *BoundaryPoints = NULL;
    315 
    316   Log() << Verbose(1) << "Begin of FindConvexBorder" << endl;
    317310
    318311  if (TesselStruct != NULL) // free if allocated
    … …  
    325318      BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
    326319  } else {
    327       Log() << Verbose(1) << "Using given boundary points set." << endl;
     320      Log() << Verbose(0) << "Using given boundary points set." << endl;
    328321  }
    329322
    … …  
    331324  for (int axis=0; axis < NDIM; axis++)
    332325    {
    333       Log() << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     326      Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
    334327      int i=0;
    335328      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    … …  
    347340    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
    348341        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
    349           Log() << Verbose(3) << "WARNING: Point " << *(runner->second.second) << " is already present!" << endl;
    350 
    351   Log() << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
     342          eLog() << Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl;
     343
     344  Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
    352345  // now we have the whole set of edge points in the BoundaryList
    353346
    … …  
    367360  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
    368361  if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
    369     Log() << Verbose(1) << "Insertion of straddling points failed!" << endl;
    370 
    371   Log() << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
     362    eLog() << Verbose(1) << "Insertion of straddling points failed!" << endl;
     363
     364  Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
    372365
    373366  // 4. Store triangles in tecplot file
    … …  
    406399        // flip the line
    407400        if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
    408           Log() << Verbose(1) << "ERROR: Correction of concave baselines failed!" << endl;
     401          eLog() << Verbose(1) << "Correction of concave baselines failed!" << endl;
    409402        else {
    410403          TesselStruct->FlipBaseline(line);
    … …  
    419412//    Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
    420413
    421   Log() << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
     414  Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
    422415
    423416  // 4. Store triangles in tecplot file
    … …  
    445438  if (BoundaryFreeFlag)
    446439    delete[] (BoundaryPoints);
    447 
    448   Log() << Verbose(1) << "End of FindConvexBorder" << endl;
    449440};
    450441
    … …  
    458449bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
    459450{
     451        Info FunctionInfo(__func__);
    460452  int i=0;
    461453  char number[MAXSTRINGSIZE];
    462454
    463455  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    464     Log() << Verbose(2) << "ERROR: TesselStruct is empty." << endl;
     456    eLog() << Verbose(1) << "TesselStruct is empty." << endl;
    465457    return false;
    466458  }
    … …  
    468460  PointMap::iterator PointRunner;
    469461  while (!TesselStruct->PointsOnBoundary.empty()) {
    470     Log() << Verbose(2) << "Remaining points are: ";
     462    Log() << Verbose(1) << "Remaining points are: ";
    471463    for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
    472464      Log() << Verbose(0) << *(PointSprinter->second) << "\t";
    … …  
    511503double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
    512504{
     505        Info FunctionInfo(__func__);
    513506  double volume = 0;
    514507  class BoundaryPointSet *point = NULL;
    … …  
    524517  int run = 0;
    525518
    526   Log() << Verbose(0) << "Begin of ConvexizeNonconvexEnvelope" << endl;
    527 
    528519  // check whether there is something to work on
    529520  if (TesselStruct == NULL) {
    530     Log() << Verbose(1) << "ERROR: TesselStruct is empty!" << endl;
     521    eLog() << Verbose(1) << "TesselStruct is empty!" << endl;
    531522    return volume;
    532523  }
    … …  
    547538      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    548539        line = LineRunner->second;
    549         Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
     540        Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    550541        if (!line->CheckConvexityCriterion()) {
    551542          // remove the point if needed
    … …  
    612603
    613604  // end
    614   Log() << Verbose(1) << "Volume is " << volume << "." << endl;
    615   Log() << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
     605  Log() << Verbose(0) << "Volume is " << volume << "." << endl;
    616606  return volume;
    617607};
    … …  
    627617double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration)
    628618{
     619        Info FunctionInfo(__func__);
    629620  bool IsAngstroem = configuration->GetIsAngstroem();
    630621  double volume = 0.;
    … …  
    633624
    634625  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
    635   Log() << Verbose(1)
    636       << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
    637       << endl;
    638626  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
    639627    { // go through every triangle, calculate volume of its pyramid with CoG as peak
    … …  
    650638      const double h = x.Norm(); // distance of CoG to triangle
    651639      const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
    652       Log() << Verbose(2) << "Area of triangle is " << setprecision(10) << G << " "
     640      Log() << Verbose(1) << "Area of triangle is " << setprecision(10) << G << " "
    653641          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
    654642          << h << " and the volume is " << PyramidVolume << " "
    … …  
    672660void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
    673661{
     662        Info FunctionInfo(__func__);
    674663  // 4. Store triangles in tecplot file
    675664  if (filename != NULL) {
    … …  
    679668      OutputName.append(TecplotSuffix);
    680669      ofstream *tecplot = new ofstream(OutputName.c_str());
    681       WriteTecplotFile(tecplot, TesselStruct, mol, 0);
     670      WriteTecplotFile(tecplot, TesselStruct, mol, -1);
    682671      tecplot->close();
    683672      delete(tecplot);
    … …  
    706695void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
    707696{
     697        Info FunctionInfo(__func__);
    708698  bool IsAngstroem = true;
    709699  double *GreatestDiameter = NULL;
    … …  
    756746  Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
    757747  if (minimumvolume > cellvolume) {
    758     eLog() << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
     748    eLog() << Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl;
    759749    Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
    760750    for (int i = 0; i < NDIM; i++)
    … …  
    806796molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
    807797{
     798        Info FunctionInfo(__func__);
    808799  molecule *Filling = new molecule(filler->elemente);
    809800  Vector CurrentPosition;
    … …  
    823814  class Tesselation *TesselStruct[List->ListOfMolecules.size()];
    824815
    825   Log() << Verbose(0) << "Begin of FillBoxWithMolecule" << endl;
    826 
    827816  i=0;
    828817  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    … …  
    849838  for(int i=0;i<NDIM;i++) {
    850839    N[i] = (int) ceil(1./FillerDistance.x[i]);
    851     Log() << Verbose(0) << N[i];
     840    Log() << Verbose(1) << N[i];
    852841    if (i != NDIM-1)
    853       Log() << Verbose(0)<< ", ";
     842      Log() << Verbose(1)<< ", ";
    854843    else
    855       Log() << Verbose(0) << "." << endl;
     844      Log() << Verbose(1) << "." << endl;
    856845  }
    857846
    … …  
    870859          // get linked cell list
    871860          if (TesselStruct[i] == NULL) {
    872             Log() << Verbose(1) << "ERROR: TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
     861            eLog() << Verbose(1) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
    873862            FillIt = false;
    874863          } else {
    … …  
    885874          for (int i=0;i<NDIM;i++)
    886875            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
    887           Log() << Verbose(3) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
     876          Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
    888877
    889878          // go through all atoms
    … …  
    946935        delete(TesselStruct[i]);
    947936  }
    948   Log() << Verbose(0) << "End of FillBoxWithMolecule" << endl;
    949 
    950937  return Filling;
    951938};
    … …  
    959946 * \param RADIUS radius of the virtual sphere
    960947 * \param *filename filename prefix for output of vertex data
    961  */
    962 void FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
     948 * \return true - tesselation successful, false - tesselation failed
     949 */
     950bool FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
    963951{
     952        Info FunctionInfo(__func__);
    964953  bool freeLC = false;
    965   LineMap::iterator baseline;
     954  bool status = false;
     955  CandidateForTesselation *baseline;
    966956  LineMap::iterator testline;
    967   bool OneLoopWithoutSuccessFlag = false;  // marks whether we went once through all baselines without finding any without two triangles
     957  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
    968958  bool TesselationFailFlag = false;
    969 
    970   Log() << Verbose(1) << "Entering search for non convex hull. " << endl;
     959  BoundaryTriangleSet *T = NULL;
     960
    971961  if (TesselStruct == NULL) {
    972962    Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
    … …  
    978968  }
    979969
    980   Log() << Verbose(0) << "Begin of FindNonConvexBorder\n";
    981 
    982970  // initialise Linked Cell
    983971  if (LCList == NULL) {
    … …  
    990978
    991979  // 2. expand from there
    992   baseline = TesselStruct->LinesOnBoundary.begin();
    993   baseline++; // skip first line
    994   while ((baseline != TesselStruct->LinesOnBoundary.end()) || (OneLoopWithoutSuccessFlag)) {
    995     if (baseline->second->triangles.size() == 1) {
    996       // 3. find next triangle
    997       TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one.
    998       OneLoopWithoutSuccessFlag = OneLoopWithoutSuccessFlag || TesselationFailFlag;
    999       if (!TesselationFailFlag)
    1000         eLog() << Verbose(0) << "WARNING: FindNextSuitableTriangle failed." << endl;
    1001 
    1002       // write temporary envelope
    1003       if (filename != NULL) {
    1004         if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
    1005           TesselStruct->Output(filename, mol);
    1006         }
     980  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
     981    // 2a. fill all new OpenLines
     982    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
     983    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
     984      Log() << Verbose(2) << *(Runner->second) << endl;
     985
     986    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
     987      baseline = Runner->second;
     988      if (baseline->pointlist.empty()) {
     989        T = (((baseline->BaseLine->triangles.begin()))->second);
     990        Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
     991        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
    1007992      }
    1008       baseline = TesselStruct->LinesOnBoundary.end();
    1009       Log() << Verbose(2) << "Baseline set to end." << endl;
    1010     } else {
    1011       //Log() << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->triangles.size() << " triangles adjacent" << endl;
    1012       if (baseline->second->triangles.size() != 2)
    1013         Log() << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
    1014     }
    1015 
    1016     if ((baseline == TesselStruct->LinesOnBoundary.end()) && (OneLoopWithoutSuccessFlag)) {
    1017       baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
     993    }
     994
     995    // 2b. search for smallest ShortestAngle among all candidates
     996    double ShortestAngle = 4.*M_PI;
     997    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
     998    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
     999      Log() << Verbose(2) << *(Runner->second) << endl;
     1000
     1001    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
     1002      if (Runner->second->ShortestAngle < ShortestAngle) {
     1003        baseline = Runner->second;
     1004        ShortestAngle = baseline->ShortestAngle;
     1005        //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
     1006      }
     1007    }
     1008    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
    10181009      OneLoopWithoutSuccessFlag = false;
    1019     }
    1020     baseline++;
    1021   }
    1022   // check envelope for consistency
    1023   CheckListOfBaselines(TesselStruct);
    1024 
    1025   // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
    1026   //->InsertStraddlingPoints(mol, LCList);
     1010    else {
     1011      TesselStruct->AddCandidateTriangle(*baseline);
     1012    }
     1013
     1014    // write temporary envelope
     1015    if (filename != NULL) {
     1016      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
     1017        TesselStruct->Output(filename, mol);
     1018      }
     1019    }
     1020  }
     1021//  // check envelope for consistency
     1022//  status = CheckListOfBaselines(TesselStruct);
     1023//
     1024//  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
     1025//  //->InsertStraddlingPoints(mol, LCList);
    10271026//  mol->GoToFirst();
    10281027//  class TesselPoint *Runner = NULL;
    … …  
    10391038//  }
    10401039
    1041   // Purges surplus triangles.
    1042   TesselStruct->RemoveDegeneratedTriangles();
     1040//  // Purges surplus triangles.
     1041//  TesselStruct->RemoveDegeneratedTriangles();
    10431042
    10441043  // check envelope for consistency
    1045   CheckListOfBaselines(TesselStruct);
     1044  status = CheckListOfBaselines(TesselStruct);
    10461045
    10471046  // write final envelope
    … …  
    10511050  if (freeLC)
    10521051    delete(LCList);
    1053   Log() << Verbose(0) << "End of FindNonConvexBorder\n";
     1052
     1053  return status;
    10541054};
    10551055
    … …  
    10631063Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
    10641064{
     1065        Info FunctionInfo(__func__);
    10651066  Vector *Center = new Vector;
    10661067  Center->Zero();
Note: See TracChangeset for help on using the changeset viewer.