Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    • Property mode changed from 100755 to 100644
    r4fc93f rb5c2d7  
    1717#include "tesselation.hpp"
    1818#include "tesselationhelpers.hpp"
     19#include "World.hpp"
    1920
    2021#include<gsl/gsl_poly.h>
     22#include<time.h>
    2123
    2224// ========================================== F U N C T I O N S =================================
    … …  
    5557  } else {
    5658    BoundaryPoints = BoundaryPtr;
    57     Log() << Verbose(0) << "Using given boundary points set." << endl;
     59    DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
    5860  }
    5961  // determine biggest "diameter" of cluster for each axis
    … …  
    161163    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
    162164
    163     Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
     165    DoLog(1) && (Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl);
    164166
    165167    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    … …  
    182184        angle = 2. * M_PI - angle;
    183185      }
    184       Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
     186      DoLog(1) && (Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
    185187      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
    186188      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
    187         Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
    188         Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
    189         Log() << Verbose(2) << "New vector: " << *Walker << endl;
     189        DoLog(2) && (Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl);
     190        DoLog(2) && (Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl);
     191        DoLog(2) && (Log() << Verbose(2) << "New vector: " << *Walker << endl);
    190192        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
    191193        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
    192194          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
    193195          BoundaryTestPair.first->second.second = Walker;
    194           Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
     196          DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
    195197        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
    196198          helper.CopyVector(&Walker->x);
    … …  
    201203          if (helper.NormSquared() < oldhelperNorm) {
    202204            BoundaryTestPair.first->second.second = Walker;
    203             Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
     205            DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl);
    204206          } else {
    205             Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl;
     207            DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl);
    206208          }
    207209        } else {
    208           Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
     210          DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
    209211        }
    210212      }
    … …  
    225227    // 3c. throw out points whose distance is less than the mean of left and right neighbours
    226228    bool flag = false;
    227     Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
     229    DoLog(1) && (Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl);
    228230    do { // do as long as we still throw one out per round
    229231      flag = false;
    … …  
    280282          const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
    281283          //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;
    282           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;
     284          DoLog(1) && (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);
    283285          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
    284286            // throw out point
    285             Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
     287            DoLog(1) && (Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl);
    286288            BoundaryPoints[axis].erase(runner);
    287289            flag = true;
    … …  
    318320      BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
    319321  } else {
    320       Log() << Verbose(0) << "Using given boundary points set." << endl;
     322      DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
    321323  }
    322324
    … …  
    324326  for (int axis=0; axis < NDIM; axis++)
    325327    {
    326       Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     328      DoLog(1) && (Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl);
    327329      int i=0;
    328330      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    329331        if (runner != BoundaryPoints[axis].begin())
    330           Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
     332          DoLog(0) && (Log() << Verbose(0) << ", " << i << ": " << *runner->second.second);
    331333        else
    332           Log() << Verbose(0) << i << ": " << *runner->second.second;
     334          DoLog(0) && (Log() << Verbose(0) << i << ": " << *runner->second.second);
    333335        i++;
    334336      }
    335       Log() << Verbose(0) << endl;
     337      DoLog(0) && (Log() << Verbose(0) << endl);
    336338    }
    337339
    … …  
    340342    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
    341343        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
    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;
     344          DoeLog(2) && (eLog()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl);
     345
     346  DoLog(0) && (Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl);
    345347  // now we have the whole set of edge points in the BoundaryList
    346348
    … …  
    360362  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
    361363  if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
    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;
     364    DoeLog(1) && (eLog()<< Verbose(1) << "Insertion of straddling points failed!" << endl);
     365
     366  DoLog(0) && (Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl);
    365367
    366368  // 4. Store triangles in tecplot file
    … …  
    393395    for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
    394396      line = LineRunner->second;
    395       Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
     397      DoLog(1) && (Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl);
    396398      if (!line->CheckConvexityCriterion()) {
    397         Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
     399        DoLog(1) && (Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl);
    398400
    399401        // flip the line
    400402        if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
    401           eLog() << Verbose(1) << "Correction of concave baselines failed!" << endl;
     403          DoeLog(1) && (eLog()<< Verbose(1) << "Correction of concave baselines failed!" << endl);
    402404        else {
    403405          TesselStruct->FlipBaseline(line);
    404           Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
     406          DoLog(1) && (Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl);
    405407        }
    406408      }
    … …  
    412414//    Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
    413415
    414   Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
     416  DoLog(0) && (Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl);
    415417
    416418  // 4. Store triangles in tecplot file
    … …  
    454456
    455457  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    456     eLog() << Verbose(1) << "TesselStruct is empty." << endl;
     458    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty." << endl);
    457459    return false;
    458460  }
    … …  
    460462  PointMap::iterator PointRunner;
    461463  while (!TesselStruct->PointsOnBoundary.empty()) {
    462     Log() << Verbose(1) << "Remaining points are: ";
     464    DoLog(1) && (Log() << Verbose(1) << "Remaining points are: ");
    463465    for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
    464       Log() << Verbose(0) << *(PointSprinter->second) << "\t";
    465     Log() << Verbose(0) << endl;
     466      DoLog(0) && (Log() << Verbose(0) << *(PointSprinter->second) << "\t");
     467    DoLog(0) && (Log() << Verbose(0) << endl);
    466468
    467469    PointRunner = TesselStruct->PointsOnBoundary.begin();
    … …  
    519521  // check whether there is something to work on
    520522  if (TesselStruct == NULL) {
    521     eLog() << Verbose(1) << "TesselStruct is empty!" << endl;
     523    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty!" << endl);
    522524    return volume;
    523525  }
    … …  
    535537      PointAdvance++;
    536538      point = PointRunner->second;
    537       Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
     539      DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
    538540      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    539541        line = LineRunner->second;
    540         Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
     542        DoLog(1) && (Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl);
    541543        if (!line->CheckConvexityCriterion()) {
    542544          // remove the point if needed
    543           Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
     545          DoLog(1) && (Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl);
    544546          volume += TesselStruct->RemovePointFromTesselatedSurface(point);
    545547          sprintf(dummy, "-first-%d", ++run);
    … …  
    562564      LineAdvance++;
    563565      line = LineRunner->second;
    564       Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
     566      DoLog(1) && (Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl);
    565567      // take highest of both lines
    566568      if (TesselStruct->IsConvexRectangle(line) == NULL) {
    … …  
    603605
    604606  // end
    605   Log() << Verbose(0) << "Volume is " << volume << "." << endl;
     607  DoLog(0) && (Log() << Verbose(0) << "Volume is " << volume << "." << endl);
    606608  return volume;
    607609};
    … …  
    654656 * \param *out output stream for debugging
    655657 * \param *mol molecule with atoms and bonds
    656  * \param *&TesselStruct Tesselation with boundary triangles
     658 * \param *TesselStruct Tesselation with boundary triangles
    657659 * \param *filename prefix of filename
    658660 * \param *extraSuffix intermediate suffix
    659661 */
    660 void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
     662void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
    661663{
    662664        Info FunctionInfo(__func__);
    … …  
    732734      totalmass += Walker->type->mass;
    733735  }
    734   Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
    735   Log() << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
     736  DoLog(0) && (Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl);
     737  DoLog(0) && (Log() << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
    736738
    737739  // solve cubic polynomial
    738   Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
     740  DoLog(1) && (Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl);
    739741  if (IsAngstroem)
    740742    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
    741743  else
    742744    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
    743   Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
     745  DoLog(1) && (Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
    744746
    745747  double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
    746   Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
     748  DoLog(1) && (Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
    747749  if (minimumvolume > cellvolume) {
    748     eLog() << Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl;
    749     Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
     750    DoeLog(1) && (eLog()<< Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl);
     751    DoLog(0) && (Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl);
    750752    for (int i = 0; i < NDIM; i++)
    751753      BoxLengths.x[i] = GreatestDiameter[i];
    … …  
    759761    double x2 = 0.;
    760762    if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1], BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
    761       Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
     763      DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl);
    762764    else {
    763       Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
     765      DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl);
    764766      x0 = x2; // sorted in ascending order
    765767    }
    … …  
    772774
    773775    // set new box dimensions
    774     Log() << Verbose(0) << "Translating to box with these boundaries." << endl;
     776    DoLog(0) && (Log() << Verbose(0) << "Translating to box with these boundaries." << endl);
    775777    mol->SetBoxDimension(&BoxLengths);
    776778    mol->CenterInBox();
    … …  
    778780  // update Box of atoms by boundary
    779781  mol->SetBoxDimension(&BoxLengths);
    780   Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
     782  DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
    781783};
    782784
    … …  
    788790 * \param *filler molecule which the box is to be filled with
    789791 * \param configuration contains box dimensions
     792 * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
    790793 * \param distance[NDIM] distance between filling molecules in each direction
     794 * \param boundary length of boundary zone between molecule and filling mollecules
     795 * \param epsilon distance to surface which is not filled
    791796 * \param RandAtomDisplacement maximum distance for random displacement per atom
    792797 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
    … …  
    794799 * \return *mol pointer to new molecule with filled atoms
    795800 */
    796 molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
     801molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
    797802{
    798803        Info FunctionInfo(__func__);
    … …  
    801806  int N[NDIM];
    802807  int n[NDIM];
    803   double *M =  ReturnFullMatrixforSymmetric(filler->cell_size);
     808  double *M =  ReturnFullMatrixforSymmetric(World::get()->cell_size);
    804809  double Rotations[NDIM*NDIM];
     810  double *MInverse = InverseMatrix(M);
    805811  Vector AtomTranslations;
    806812  Vector FillerTranslations;
    807813  Vector FillerDistance;
     814  Vector Inserter;
    808815  double FillIt = false;
    809816  atom *Walker = NULL;
    810817  bond *Binder = NULL;
    811   int i = 0;
    812   LinkedCell *LCList[List->ListOfMolecules.size()];
    813818  double phi[NDIM];
    814   class Tesselation *TesselStruct[List->ListOfMolecules.size()];
    815 
    816   i=0;
    817   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    818     Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
    819     LCList[i] = new LinkedCell((*ListRunner), 5.); // get linked cell list
    820     if (TesselStruct[i] == NULL) {
    821       Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
    822       FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
    823     }
    824     i++;
    825   }
     819  map<molecule *, Tesselation *> TesselStruct;
     820  map<molecule *, LinkedCell *> LCList;
     821
     822  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
     823    if ((*ListRunner)->AtomCount > 0) {
     824      DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl);
     825      LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list
     826      DoLog(1) && (Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl);
     827      TesselStruct[(*ListRunner)] = NULL;
     828      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell *&)LCList[(*ListRunner)], 5., NULL);
     829    }
    826830
    827831  // Center filler at origin
    828   filler->CenterOrigin();
     832  filler->CenterEdge(&Inserter);
    829833  filler->Center.Zero();
     834  DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
     835  Binder = filler->first;
     836  while(Binder->next != filler->last) {
     837    Binder = Binder->next;
     838    DoLog(2) && (Log() << Verbose(2) << "  " << *Binder << endl);
     839  }
    830840
    831841  filler->CountAtoms();
    … …  
    835845  FillerDistance.Init(distance[0], distance[1], distance[2]);
    836846  FillerDistance.InverseMatrixMultiplication(M);
    837   Log() << Verbose(1) << "INFO: Grid steps are ";
    838   for(int i=0;i<NDIM;i++) {
     847  for(int i=0;i<NDIM;i++)
    839848    N[i] = (int) ceil(1./FillerDistance.x[i]);
    840     Log() << Verbose(1) << N[i];
    841     if (i != NDIM-1)
    842       Log() << Verbose(1)<< ", ";
    843     else
    844       Log() << Verbose(1) << "." << endl;
    845   }
     849  DoLog(1) && (Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl);
     850
     851  // initialize seed of random number generator to current time
     852  srand ( time(NULL) );
    846853
    847854  // go over [0,1]^3 filler grid
    … …  
    852859        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    853860        CurrentPosition.MatrixMultiplication(M);
    854         Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;
    855         // Check whether point is in- or outside
    856         FillIt = true;
    857         i=0;
    858         for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    859           // get linked cell list
    860           if (TesselStruct[i] == NULL) {
    861             eLog() << Verbose(1) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
    862             FillIt = false;
     861        // create molecule random translation vector ...
     862        for (int i=0;i<NDIM;i++)
     863          FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
     864        DoLog(2) && (Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl);
     865
     866        // go through all atoms
     867        for (int i=0;i<filler->AtomCount;i++)
     868          CopyAtoms[i] = NULL;
     869        Walker = filler->start;
     870        while (Walker->next != filler->end) {
     871          Walker = Walker->next;
     872
     873          // create atomic random translation vector ...
     874          for (int i=0;i<NDIM;i++)
     875            AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
     876
     877          // ... and rotation matrix
     878          if (DoRandomRotation) {
     879            for (int i=0;i<NDIM;i++) {
     880              phi[i] = rand()/(RAND_MAX/(2.*M_PI));
     881            }
     882
     883            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
     884            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
     885            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
     886            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
     887            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
     888            Rotations[7] =               sin(phi[1])                                                ;
     889            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
     890            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
     891            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
     892          }
     893
     894          // ... and put at new position
     895          Inserter.CopyVector(&(Walker->x));
     896          if (DoRandomRotation)
     897            Inserter.MatrixMultiplication(Rotations);
     898          Inserter.AddVector(&AtomTranslations);
     899          Inserter.AddVector(&FillerTranslations);
     900          Inserter.AddVector(&CurrentPosition);
     901
     902          // check whether inserter is inside box
     903          Inserter.MatrixMultiplication(MInverse);
     904          FillIt = true;
     905          for (int i=0;i<NDIM;i++)
     906            FillIt = FillIt && (Inserter.x[i] >= -MYEPSILON) && ((Inserter.x[i]-1.) <= MYEPSILON);
     907          Inserter.MatrixMultiplication(M);
     908
     909          // Check whether point is in- or outside
     910          for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
     911            // get linked cell list
     912            if (TesselStruct[(*ListRunner)] != NULL) {
     913              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
     914              FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
     915            }
     916          }
     917          // insert into Filling
     918          if (FillIt) {
     919            DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl);
     920            // copy atom ...
     921            CopyAtoms[Walker->nr] = new atom(Walker);
     922            CopyAtoms[Walker->nr]->x.CopyVector(&Inserter);
     923            Filling->AddAtom(CopyAtoms[Walker->nr]);
     924            DoLog(4) && (Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl);
    863925          } else {
    864             FillIt = FillIt && (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
    865             i++;
     926            DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl);
     927            CopyAtoms[Walker->nr] = NULL;
     928            continue;
    866929          }
    867930        }
    868 
    869         if (FillIt) {
    870           // fill in Filler
    871           Log() << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
    872 
    873           // create molecule random translation vector ...
    874           for (int i=0;i<NDIM;i++)
    875             FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
    876           Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
    877 
    878           // go through all atoms
    879           Walker = filler->start;
    880           while (Walker->next != filler->end) {
    881             Walker = Walker->next;
    882             // copy atom ...
    883             CopyAtoms[Walker->nr] = new atom(Walker);
    884 
    885             // create atomic random translation vector ...
    886             for (int i=0;i<NDIM;i++)
    887               AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
    888 
    889             // ... and rotation matrix
    890             if (DoRandomRotation) {
    891               for (int i=0;i<NDIM;i++) {
    892                 phi[i] = rand()/(RAND_MAX/(2.*M_PI));
    893               }
    894 
    895               Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
    896               Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
    897               Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
    898               Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
    899               Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
    900               Rotations[7] =               sin(phi[1])                                                ;
    901               Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
    902               Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
    903               Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
    904             }
    905 
    906             // ... and put at new position
    907             if (DoRandomRotation)
    908               CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
    909             CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
    910             CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
    911             CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
    912 
    913             // insert into Filling
    914 
    915             // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
    916             Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
    917             Filling->AddAtom(CopyAtoms[Walker->nr]);
    918           }
    919 
    920           // go through all bonds and add as well
    921           Binder = filler->first;
    922           while(Binder->next != filler->last) {
    923             Binder = Binder->next;
    924             Log() << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
     931        // go through all bonds and add as well
     932        Binder = filler->first;
     933        while(Binder->next != filler->last) {
     934          Binder = Binder->next;
     935          if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
     936            Log()  << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
    925937            Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
    926938          }
    927         } else {
    928           // leave empty
    929           Log() << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
    930939        }
    931940      }
    932941  Free(&M);
    933   for (size_t i=0;i<List->ListOfMolecules.size();i++) {
    934         delete(LCList[i]);
    935         delete(TesselStruct[i]);
    936   }
     942  Free(&MInverse);
     943
    937944  return Filling;
    938945};
    … …  
    953960  bool freeLC = false;
    954961  bool status = false;
    955   CandidateForTesselation *baseline;
    956   LineMap::iterator testline;
     962  CandidateForTesselation *baseline = NULL;
    957963  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
    958964  bool TesselationFailFlag = false;
    959   BoundaryTriangleSet *T = NULL;
    960965
    961966  if (TesselStruct == NULL) {
    962     Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
     967    DoLog(1) && (Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl);
    963968    TesselStruct= new Tesselation;
    964969  } else {
    965970    delete(TesselStruct);
    966     Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
     971    DoLog(1) && (Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl);
    967972    TesselStruct = new Tesselation;
    968973  }
    … …  
    975980
    976981  // 1. get starting triangle
    977   TesselStruct->FindStartingTriangle(RADIUS, LCList);
     982  if (!TesselStruct->FindStartingTriangle(RADIUS, LCList)) {
     983    DoeLog(0) && (eLog() << Verbose(0) << "No valid starting triangle found." << endl);
     984    //performCriticalExit();
     985  }
     986  if (filename != NULL) {
     987    if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
     988      TesselStruct->Output(filename, mol);
     989    }
     990  }
    978991
    979992  // 2. expand from there
    980993  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;
     994    (cerr << "There are " <<  TesselStruct->TrianglesOnBoundary.size() << " triangles and " << TesselStruct->OpenLines.size() << " open lines to scan for candidates." << endl);
     995    // 2a. print OpenLines without candidates
     996    DoLog(1) && (Log() << Verbose(1) << "There are the following open lines to scan for a candidates:" << endl);
    983997    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.
    992       }
    993     }
    994 
    995     // 2b. search for smallest ShortestAngle among all candidates
     998      if (Runner->second->pointlist.empty())
     999        DoLog(1) && (Log() << Verbose(1) << " " << *(Runner->second) << endl);
     1000
     1001    // 2b. find best candidate for each OpenLine
     1002    TesselationFailFlag = TesselStruct->FindCandidatesforOpenLines(RADIUS, LCList);
     1003
     1004    // 2c. print OpenLines with candidates again
     1005    DoLog(1) && (Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl);
     1006    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
     1007      DoLog(1) && (Log() << Verbose(1) << " " << *(Runner->second) << endl);
     1008
     1009    // 2d. search for smallest ShortestAngle among all candidates
    9961010    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 
    10011011    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
    10021012      if (Runner->second->ShortestAngle < ShortestAngle) {
    10031013        baseline = Runner->second;
    10041014        ShortestAngle = baseline->ShortestAngle;
    1005         //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
     1015        DoLog(1) && (Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl);
    10061016      }
    10071017    }
     1018    // 2e. if we found one, add candidate
    10081019    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
    10091020      OneLoopWithoutSuccessFlag = false;
    10101021    else {
    1011       TesselStruct->AddCandidateTriangle(*baseline);
    1012     }
    1013 
    1014     // write temporary envelope
     1022      TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
     1023    }
     1024
     1025    // 2f. write temporary envelope
    10151026    if (filename != NULL) {
    10161027      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
    … …  
    10441055  status = CheckListOfBaselines(TesselStruct);
    10451056
     1057  // store before correction
     1058  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
     1059
     1060//  // correct degenerated polygons
     1061//  TesselStruct->CorrectAllDegeneratedPolygons();
     1062//
     1063//  // check envelope for consistency
     1064//  status = CheckListOfBaselines(TesselStruct);
     1065
    10461066  // write final envelope
    10471067  CalculateConcavityPerBoundaryPoint(TesselStruct);
Note: See TracChangeset for help on using the changeset viewer.