Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    • Property mode changed from 100644 to 100755
    rb5c2d7 r4fc93f  
    1717#include "tesselation.hpp"
    1818#include "tesselationhelpers.hpp"
    19 #include "World.hpp"
    2019
    2120#include<gsl/gsl_poly.h>
    22 #include<time.h>
    2321
    2422// ========================================== F U N C T I O N S =================================
    … …  
    5755  } else {
    5856    BoundaryPoints = BoundaryPtr;
    59     DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
     57    Log() << Verbose(0) << "Using given boundary points set." << endl;
    6058  }
    6159  // determine biggest "diameter" of cluster for each axis
    … …  
    163161    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
    164162
    165     DoLog(1) && (Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl);
     163    Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
    166164
    167165    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    … …  
    184182        angle = 2. * M_PI - angle;
    185183      }
    186       DoLog(1) && (Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
     184      Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
    187185      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
    188186      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
    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);
     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;
    192190        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
    193191        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
    194192          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
    195193          BoundaryTestPair.first->second.second = Walker;
    196           DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
     194          Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
    197195        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
    198196          helper.CopyVector(&Walker->x);
    … …  
    203201          if (helper.NormSquared() < oldhelperNorm) {
    204202            BoundaryTestPair.first->second.second = Walker;
    205             DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl);
     203            Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
    206204          } else {
    207             DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl);
     205            Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl;
    208206          }
    209207        } else {
    210           DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
     208          Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
    211209        }
    212210      }
    … …  
    227225    // 3c. throw out points whose distance is less than the mean of left and right neighbours
    228226    bool flag = false;
    229     DoLog(1) && (Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl);
     227    Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
    230228    do { // do as long as we still throw one out per round
    231229      flag = false;
    … …  
    282280          const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
    283281          //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;
    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);
     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;
    285283          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
    286284            // throw out point
    287             DoLog(1) && (Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl);
     285            Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
    288286            BoundaryPoints[axis].erase(runner);
    289287            flag = true;
    … …  
    320318      BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
    321319  } else {
    322       DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
     320      Log() << Verbose(0) << "Using given boundary points set." << endl;
    323321  }
    324322
    … …  
    326324  for (int axis=0; axis < NDIM; axis++)
    327325    {
    328       DoLog(1) && (Log() << Verbose(1) << "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;
    329327      int i=0;
    330328      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    331329        if (runner != BoundaryPoints[axis].begin())
    332           DoLog(0) && (Log() << Verbose(0) << ", " << i << ": " << *runner->second.second);
     330          Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
    333331        else
    334           DoLog(0) && (Log() << Verbose(0) << i << ": " << *runner->second.second);
     332          Log() << Verbose(0) << i << ": " << *runner->second.second;
    335333        i++;
    336334      }
    337       DoLog(0) && (Log() << Verbose(0) << endl);
     335      Log() << Verbose(0) << endl;
    338336    }
    339337
    … …  
    342340    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
    343341        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
    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);
     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;
    347345  // now we have the whole set of edge points in the BoundaryList
    348346
    … …  
    362360  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
    363361  if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
    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);
     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;
    367365
    368366  // 4. Store triangles in tecplot file
    … …  
    395393    for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
    396394      line = LineRunner->second;
    397       DoLog(1) && (Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl);
     395      Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
    398396      if (!line->CheckConvexityCriterion()) {
    399         DoLog(1) && (Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl);
     397        Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
    400398
    401399        // flip the line
    402400        if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
    403           DoeLog(1) && (eLog()<< Verbose(1) << "Correction of concave baselines failed!" << endl);
     401          eLog() << Verbose(1) << "Correction of concave baselines failed!" << endl;
    404402        else {
    405403          TesselStruct->FlipBaseline(line);
    406           DoLog(1) && (Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl);
     404          Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
    407405        }
    408406      }
    … …  
    414412//    Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
    415413
    416   DoLog(0) && (Log() << Verbose(0) << "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;
    417415
    418416  // 4. Store triangles in tecplot file
    … …  
    456454
    457455  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    458     DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty." << endl);
     456    eLog() << Verbose(1) << "TesselStruct is empty." << endl;
    459457    return false;
    460458  }
    … …  
    462460  PointMap::iterator PointRunner;
    463461  while (!TesselStruct->PointsOnBoundary.empty()) {
    464     DoLog(1) && (Log() << Verbose(1) << "Remaining points are: ");
     462    Log() << Verbose(1) << "Remaining points are: ";
    465463    for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
    466       DoLog(0) && (Log() << Verbose(0) << *(PointSprinter->second) << "\t");
    467     DoLog(0) && (Log() << Verbose(0) << endl);
     464      Log() << Verbose(0) << *(PointSprinter->second) << "\t";
     465    Log() << Verbose(0) << endl;
    468466
    469467    PointRunner = TesselStruct->PointsOnBoundary.begin();
    … …  
    521519  // check whether there is something to work on
    522520  if (TesselStruct == NULL) {
    523     DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty!" << endl);
     521    eLog() << Verbose(1) << "TesselStruct is empty!" << endl;
    524522    return volume;
    525523  }
    … …  
    537535      PointAdvance++;
    538536      point = PointRunner->second;
    539       DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
     537      Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
    540538      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    541539        line = LineRunner->second;
    542         DoLog(1) && (Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl);
     540        Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    543541        if (!line->CheckConvexityCriterion()) {
    544542          // remove the point if needed
    545           DoLog(1) && (Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl);
     543          Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
    546544          volume += TesselStruct->RemovePointFromTesselatedSurface(point);
    547545          sprintf(dummy, "-first-%d", ++run);
    … …  
    564562      LineAdvance++;
    565563      line = LineRunner->second;
    566       DoLog(1) && (Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl);
     564      Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
    567565      // take highest of both lines
    568566      if (TesselStruct->IsConvexRectangle(line) == NULL) {
    … …  
    605603
    606604  // end
    607   DoLog(0) && (Log() << Verbose(0) << "Volume is " << volume << "." << endl);
     605  Log() << Verbose(0) << "Volume is " << volume << "." << endl;
    608606  return volume;
    609607};
    … …  
    656654 * \param *out output stream for debugging
    657655 * \param *mol molecule with atoms and bonds
    658  * \param *TesselStruct Tesselation with boundary triangles
     656 * \param *&TesselStruct Tesselation with boundary triangles
    659657 * \param *filename prefix of filename
    660658 * \param *extraSuffix intermediate suffix
    661659 */
    662 void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
     660void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
    663661{
    664662        Info FunctionInfo(__func__);
    … …  
    734732      totalmass += Walker->type->mass;
    735733  }
    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);
     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;
    738736
    739737  // solve cubic polynomial
    740   DoLog(1) && (Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl);
     738  Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
    741739  if (IsAngstroem)
    742740    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
    743741  else
    744742    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
    745   DoLog(1) && (Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
     743  Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
    746744
    747745  double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
    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);
     746  Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
    749747  if (minimumvolume > cellvolume) {
    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);
     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;
    752750    for (int i = 0; i < NDIM; i++)
    753751      BoxLengths.x[i] = GreatestDiameter[i];
    … …  
    761759    double x2 = 0.;
    762760    if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1], BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
    763       DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl);
     761      Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
    764762    else {
    765       DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl);
     763      Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
    766764      x0 = x2; // sorted in ascending order
    767765    }
    … …  
    774772
    775773    // set new box dimensions
    776     DoLog(0) && (Log() << Verbose(0) << "Translating to box with these boundaries." << endl);
     774    Log() << Verbose(0) << "Translating to box with these boundaries." << endl;
    777775    mol->SetBoxDimension(&BoxLengths);
    778776    mol->CenterInBox();
    … …  
    780778  // update Box of atoms by boundary
    781779  mol->SetBoxDimension(&BoxLengths);
    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);
     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;
    783781};
    784782
    … …  
    790788 * \param *filler molecule which the box is to be filled with
    791789 * \param configuration contains box dimensions
    792  * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
    793790 * \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
    796791 * \param RandAtomDisplacement maximum distance for random displacement per atom
    797792 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
    … …  
    799794 * \return *mol pointer to new molecule with filled atoms
    800795 */
    801 molecule * 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)
     796molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
    802797{
    803798        Info FunctionInfo(__func__);
    … …  
    806801  int N[NDIM];
    807802  int n[NDIM];
    808   double *M =  ReturnFullMatrixforSymmetric(World::get()->cell_size);
     803  double *M =  ReturnFullMatrixforSymmetric(filler->cell_size);
    809804  double Rotations[NDIM*NDIM];
    810   double *MInverse = InverseMatrix(M);
    811805  Vector AtomTranslations;
    812806  Vector FillerTranslations;
    813807  Vector FillerDistance;
    814   Vector Inserter;
    815808  double FillIt = false;
    816809  atom *Walker = NULL;
    817810  bond *Binder = NULL;
     811  int i = 0;
     812  LinkedCell *LCList[List->ListOfMolecules.size()];
    818813  double phi[NDIM];
    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     }
     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  }
    830826
    831827  // Center filler at origin
    832   filler->CenterEdge(&Inserter);
     828  filler->CenterOrigin();
    833829  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   }
    840830
    841831  filler->CountAtoms();
    … …  
    845835  FillerDistance.Init(distance[0], distance[1], distance[2]);
    846836  FillerDistance.InverseMatrixMultiplication(M);
    847   for(int i=0;i<NDIM;i++)
     837  Log() << Verbose(1) << "INFO: Grid steps are ";
     838  for(int i=0;i<NDIM;i++) {
    848839    N[i] = (int) ceil(1./FillerDistance.x[i]);
    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) );
     840    Log() << Verbose(1) << N[i];
     841    if (i != NDIM-1)
     842      Log() << Verbose(1)<< ", ";
     843    else
     844      Log() << Verbose(1) << "." << endl;
     845  }
    853846
    854847  // go over [0,1]^3 filler grid
    … …  
    859852        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    860853        CurrentPosition.MatrixMultiplication(M);
    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 ...
     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;
     863          } else {
     864            FillIt = FillIt && (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
     865            i++;
     866          }
     867        }
     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 ...
    874874          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);
     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;
    920882            // copy atom ...
    921883            CopyAtoms[Walker->nr] = new atom(Walker);
    922             CopyAtoms[Walker->nr]->x.CopyVector(&Inserter);
     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;
    923917            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);
    925           } else {
    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;
    929918          }
    930         }
    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;
     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;
    937925            Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
    938926          }
     927        } else {
     928          // leave empty
     929          Log() << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
    939930        }
    940931      }
    941932  Free(&M);
    942   Free(&MInverse);
    943 
     933  for (size_t i=0;i<List->ListOfMolecules.size();i++) {
     934        delete(LCList[i]);
     935        delete(TesselStruct[i]);
     936  }
    944937  return Filling;
    945938};
    … …  
    960953  bool freeLC = false;
    961954  bool status = false;
    962   CandidateForTesselation *baseline = NULL;
     955  CandidateForTesselation *baseline;
     956  LineMap::iterator testline;
    963957  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
    964958  bool TesselationFailFlag = false;
     959  BoundaryTriangleSet *T = NULL;
    965960
    966961  if (TesselStruct == NULL) {
    967     DoLog(1) && (Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl);
     962    Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
    968963    TesselStruct= new Tesselation;
    969964  } else {
    970965    delete(TesselStruct);
    971     DoLog(1) && (Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl);
     966    Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
    972967    TesselStruct = new Tesselation;
    973968  }
    … …  
    980975
    981976  // 1. get starting triangle
    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   }
     977  TesselStruct->FindStartingTriangle(RADIUS, LCList);
    991978
    992979  // 2. expand from there
    993980  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
    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);
     981    // 2a. fill all new OpenLines
     982    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
    997983    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
    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);
     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
     996    double ShortestAngle = 4.*M_PI;
     997    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
    1006998    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
    1010     double ShortestAngle = 4.*M_PI;
     999      Log() << Verbose(2) << *(Runner->second) << endl;
     1000
    10111001    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
    10121002      if (Runner->second->ShortestAngle < ShortestAngle) {
    10131003        baseline = Runner->second;
    10141004        ShortestAngle = baseline->ShortestAngle;
    1015         DoLog(1) && (Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl);
     1005        //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
    10161006      }
    10171007    }
    1018     // 2e. if we found one, add candidate
    10191008    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
    10201009      OneLoopWithoutSuccessFlag = false;
    10211010    else {
    1022       TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
    1023     }
    1024 
    1025     // 2f. write temporary envelope
     1011      TesselStruct->AddCandidateTriangle(*baseline);
     1012    }
     1013
     1014    // write temporary envelope
    10261015    if (filename != NULL) {
    10271016      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
    … …  
    10551044  status = CheckListOfBaselines(TesselStruct);
    10561045
    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 
    10661046  // write final envelope
    10671047  CalculateConcavityPerBoundaryPoint(TesselStruct);
Note: See TracChangeset for help on using the changeset viewer.