Changeset 5f612ee for src/boundary.cpp


Ignore:
Timestamp:
Apr 27, 2010, 2:25:42 PM (15 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
632bc3
Parents:
13d5a9 (diff), c695c9 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'Analysis_PairCorrelation' into StructureRefactoring

Conflicts:

molecuilder/src/Makefile.am
molecuilder/src/World.cpp
molecuilder/src/World.hpp
molecuilder/src/boundary.cpp
molecuilder/src/builder.cpp
molecuilder/src/log.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/periodentafel.cpp
molecuilder/src/tesselation.cpp
molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
molecuilder/src/unittests/Makefile.am
molecuilder/src/unittests/bondgraphunittest.cpp
molecuilder/src/unittests/gslvectorunittest.cpp
molecuilder/src/unittests/logunittest.cpp
molecuilder/src/unittests/tesselation_boundarytriangleunittest.hpp
molecuilder/src/vector.cpp
molecuilder/tests/Tesselations/defs.in

Conflicts have been many and too numerous to listen here, just the few general cases

  • new molecule() replaced by World::getInstance().createMolecule()
  • new atom() replaced by World::getInstance().createAtom() where appropriate.
  • Some DoLog()s added interfered with changes to the message produced by Log() << Verbose(.) << ...
  • DoLog() has been erroneously added to TestRunner.cpp as well, there cout is appropriate
  • ...

Additionally, there was a bug in atom::clone(), sort was set to atom::nr of the atom to clone not of the clone itself. This caused a failure of the fragmentation.

This merge has been fully checked from a clean build directory with subsequent configure,make all install and make check.
It configures, compiles and runs all test cases and the test suite without errors.

Signed-off-by: Frederik Heber <heber@…>

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    • Property mode changed from 100755 to 100644
    r13d5a9 r5f612ee  
    1818#include "tesselation.hpp"
    1919#include "tesselationhelpers.hpp"
     20#include "World.hpp"
    2021
    2122#include<gsl/gsl_poly.h>
    … …  
    5758  } else {
    5859    BoundaryPoints = BoundaryPtr;
    59     Log() << Verbose(0) << "Using given boundary points set." << endl;
     60    DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
    6061  }
    6162  // determine biggest "diameter" of cluster for each axis
    … …  
    163164    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
    164165
    165     Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
     166    DoLog(1) && (Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl);
    166167
    167168    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    … …  
    184185        angle = 2. * M_PI - angle;
    185186      }
    186       Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
     187      DoLog(1) && (Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
    187188      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
    188189      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
    189         Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
    190         Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
    191         Log() << Verbose(2) << "New vector: " << *Walker << endl;
     190        DoLog(2) && (Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl);
     191        DoLog(2) && (Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl);
     192        DoLog(2) && (Log() << Verbose(2) << "New vector: " << *Walker << endl);
    192193        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
    193194        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
    194195          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
    195196          BoundaryTestPair.first->second.second = Walker;
    196           Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
     197          DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
    197198        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
    198199          helper.CopyVector(&Walker->x);
    … …  
    203204          if (helper.NormSquared() < oldhelperNorm) {
    204205            BoundaryTestPair.first->second.second = Walker;
    205             Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
     206            DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl);
    206207          } else {
    207             Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl;
     208            DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl);
    208209          }
    209210        } else {
    210           Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
     211          DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
    211212        }
    212213      }
    … …  
    227228    // 3c. throw out points whose distance is less than the mean of left and right neighbours
    228229    bool flag = false;
    229     Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
     230    DoLog(1) && (Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl);
    230231    do { // do as long as we still throw one out per round
    231232      flag = false;
    … …  
    282283          const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
    283284          //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           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;
     285          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);
    285286          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
    286287            // throw out point
    287             Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
     288            DoLog(1) && (Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl);
    288289            BoundaryPoints[axis].erase(runner);
    289290            flag = true;
    … …  
    320321      BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
    321322  } else {
    322       Log() << Verbose(0) << "Using given boundary points set." << endl;
     323      DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
    323324  }
    324325
    … …  
    326327  for (int axis=0; axis < NDIM; axis++)
    327328    {
    328       Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
     329      DoLog(1) && (Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl);
    329330      int i=0;
    330331      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    331332        if (runner != BoundaryPoints[axis].begin())
    332           Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
     333          DoLog(0) && (Log() << Verbose(0) << ", " << i << ": " << *runner->second.second);
    333334        else
    334           Log() << Verbose(0) << i << ": " << *runner->second.second;
     335          DoLog(0) && (Log() << Verbose(0) << i << ": " << *runner->second.second);
    335336        i++;
    336337      }
    337       Log() << Verbose(0) << endl;
     338      DoLog(0) && (Log() << Verbose(0) << endl);
    338339    }
    339340
    … …  
    342343    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
    343344        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
    344           eLog() << Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl;
    345 
    346   Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
     345          DoeLog(2) && (eLog()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl);
     346
     347  DoLog(0) && (Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl);
    347348  // now we have the whole set of edge points in the BoundaryList
    348349
    … …  
    362363  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
    363364  if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
    364     eLog() << Verbose(1) << "Insertion of straddling points failed!" << endl;
    365 
    366   Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
     365    DoeLog(1) && (eLog()<< Verbose(1) << "Insertion of straddling points failed!" << endl);
     366
     367  DoLog(0) && (Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl);
    367368
    368369  // 4. Store triangles in tecplot file
    … …  
    395396    for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
    396397      line = LineRunner->second;
    397       Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
     398      DoLog(1) && (Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl);
    398399      if (!line->CheckConvexityCriterion()) {
    399         Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
     400        DoLog(1) && (Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl);
    400401
    401402        // flip the line
    402403        if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
    403           eLog() << Verbose(1) << "Correction of concave baselines failed!" << endl;
     404          DoeLog(1) && (eLog()<< Verbose(1) << "Correction of concave baselines failed!" << endl);
    404405        else {
    405406          TesselStruct->FlipBaseline(line);
    406           Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
     407          DoLog(1) && (Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl);
    407408        }
    408409      }
    … …  
    414415//    Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
    415416
    416   Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
     417  DoLog(0) && (Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl);
    417418
    418419  // 4. Store triangles in tecplot file
    … …  
    456457
    457458  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    458     eLog() << Verbose(1) << "TesselStruct is empty." << endl;
     459    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty." << endl);
    459460    return false;
    460461  }
    … …  
    462463  PointMap::iterator PointRunner;
    463464  while (!TesselStruct->PointsOnBoundary.empty()) {
    464     Log() << Verbose(1) << "Remaining points are: ";
     465    DoLog(1) && (Log() << Verbose(1) << "Remaining points are: ");
    465466    for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
    466       Log() << Verbose(0) << *(PointSprinter->second) << "\t";
    467     Log() << Verbose(0) << endl;
     467      DoLog(0) && (Log() << Verbose(0) << *(PointSprinter->second) << "\t");
     468    DoLog(0) && (Log() << Verbose(0) << endl);
    468469
    469470    PointRunner = TesselStruct->PointsOnBoundary.begin();
    … …  
    521522  // check whether there is something to work on
    522523  if (TesselStruct == NULL) {
    523     eLog() << Verbose(1) << "TesselStruct is empty!" << endl;
     524    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty!" << endl);
    524525    return volume;
    525526  }
    … …  
    537538      PointAdvance++;
    538539      point = PointRunner->second;
    539       Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
     540      DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
    540541      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    541542        line = LineRunner->second;
    542         Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
     543        DoLog(1) && (Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl);
    543544        if (!line->CheckConvexityCriterion()) {
    544545          // remove the point if needed
    545           Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
     546          DoLog(1) && (Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl);
    546547          volume += TesselStruct->RemovePointFromTesselatedSurface(point);
    547548          sprintf(dummy, "-first-%d", ++run);
    … …  
    564565      LineAdvance++;
    565566      line = LineRunner->second;
    566       Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
     567      DoLog(1) && (Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl);
    567568      // take highest of both lines
    568569      if (TesselStruct->IsConvexRectangle(line) == NULL) {
    … …  
    605606
    606607  // end
    607   Log() << Verbose(0) << "Volume is " << volume << "." << endl;
     608  DoLog(0) && (Log() << Verbose(0) << "Volume is " << volume << "." << endl);
    608609  return volume;
    609610};
    … …  
    734735      totalmass += Walker->type->mass;
    735736  }
    736   Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
    737   Log() << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
     737  DoLog(0) && (Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl);
     738  DoLog(0) && (Log() << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
    738739
    739740  // solve cubic polynomial
    740   Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
     741  DoLog(1) && (Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl);
    741742  if (IsAngstroem)
    742743    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
    743744  else
    744745    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
    745   Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
     746  DoLog(1) && (Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
    746747
    747748  double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
    748   Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
     749  DoLog(1) && (Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
    749750  if (minimumvolume > cellvolume) {
    750     eLog() << Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl;
    751     Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
     751    DoeLog(1) && (eLog()<< Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl);
     752    DoLog(0) && (Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl);
    752753    for (int i = 0; i < NDIM; i++)
    753754      BoxLengths.x[i] = GreatestDiameter[i];
    … …  
    761762    double x2 = 0.;
    762763    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       Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
     764      DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl);
    764765    else {
    765       Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
     766      DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl);
    766767      x0 = x2; // sorted in ascending order
    767768    }
    … …  
    774775
    775776    // set new box dimensions
    776     Log() << Verbose(0) << "Translating to box with these boundaries." << endl;
     777    DoLog(0) && (Log() << Verbose(0) << "Translating to box with these boundaries." << endl);
    777778    mol->SetBoxDimension(&BoxLengths);
    778779    mol->CenterInBox();
    … …  
    780781  // update Box of atoms by boundary
    781782  mol->SetBoxDimension(&BoxLengths);
    782   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;
     783  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);
    783784};
    784785
    … …  
    790791 * \param *filler molecule which the box is to be filled with
    791792 * \param configuration contains box dimensions
     793 * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
    792794 * \param distance[NDIM] distance between filling molecules in each direction
    793795 * \param boundary length of boundary zone between molecule and filling mollecules
    … …  
    798800 * \return *mol pointer to new molecule with filled atoms
    799801 */
    800 molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
     802molecule * 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)
    801803{
    802804        Info FunctionInfo(__func__);
    … …  
    805807  int N[NDIM];
    806808  int n[NDIM];
    807   double *M =  ReturnFullMatrixforSymmetric(filler->cell_size);
     809  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    808810  double Rotations[NDIM*NDIM];
     811  double *MInverse = InverseMatrix(M);
    809812  Vector AtomTranslations;
    810813  Vector FillerTranslations;
    811814  Vector FillerDistance;
     815  Vector Inserter;
    812816  double FillIt = false;
    813817  atom *Walker = NULL;
    814818  bond *Binder = NULL;
    815   int i = 0;
    816   LinkedCell *LCList[List->ListOfMolecules.size()];
    817819  double phi[NDIM];
    818   class Tesselation *TesselStruct[List->ListOfMolecules.size()];
    819 
    820   i=0;
    821   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    822     Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
    823     LCList[i] = new LinkedCell((*ListRunner), 10.); // get linked cell list
    824     Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
    825     TesselStruct[i] = NULL;
    826     FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
    827     i++;
    828   }
     820  map<molecule *, Tesselation *> TesselStruct;
     821  map<molecule *, LinkedCell *> LCList;
     822
     823  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
     824    if ((*ListRunner)->AtomCount > 0) {
     825      DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl);
     826      LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list
     827      DoLog(1) && (Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl);
     828      TesselStruct[(*ListRunner)] = NULL;
     829      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell *&)LCList[(*ListRunner)], 5., NULL);
     830    }
    829831
    830832  // Center filler at origin
    831   filler->CenterOrigin();
     833  filler->CenterEdge(&Inserter);
    832834  filler->Center.Zero();
     835  DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
     836  Binder = filler->first;
     837  while(Binder->next != filler->last) {
     838    Binder = Binder->next;
     839    DoLog(2) && (Log() << Verbose(2) << "  " << *Binder << endl);
     840  }
    833841
    834842  filler->CountAtoms();
    … …  
    840848  for(int i=0;i<NDIM;i++)
    841849    N[i] = (int) ceil(1./FillerDistance.x[i]);
    842   Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl;
     850  DoLog(1) && (Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl);
    843851
    844852  // initialize seed of random number generator to current time
    … …  
    852860        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    853861        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(0) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
    862             FillIt = false;
     862        // create molecule random translation vector ...
     863        for (int i=0;i<NDIM;i++)
     864          FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
     865        DoLog(2) && (Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl);
     866
     867        // go through all atoms
     868        for (int i=0;i<filler->AtomCount;i++)
     869          CopyAtoms[i] = NULL;
     870        Walker = filler->start;
     871        while (Walker->next != filler->end) {
     872          Walker = Walker->next;
     873
     874          // create atomic random translation vector ...
     875          for (int i=0;i<NDIM;i++)
     876            AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
     877
     878          // ... and rotation matrix
     879          if (DoRandomRotation) {
     880            for (int i=0;i<NDIM;i++) {
     881              phi[i] = rand()/(RAND_MAX/(2.*M_PI));
     882            }
     883
     884            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
     885            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
     886            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
     887            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
     888            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
     889            Rotations[7] =               sin(phi[1])                                                ;
     890            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
     891            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
     892            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
     893          }
     894
     895          // ... and put at new position
     896          Inserter.CopyVector(&(Walker->x));
     897          if (DoRandomRotation)
     898            Inserter.MatrixMultiplication(Rotations);
     899          Inserter.AddVector(&AtomTranslations);
     900          Inserter.AddVector(&FillerTranslations);
     901          Inserter.AddVector(&CurrentPosition);
     902
     903          // check whether inserter is inside box
     904          Inserter.MatrixMultiplication(MInverse);
     905          FillIt = true;
     906          for (int i=0;i<NDIM;i++)
     907            FillIt = FillIt && (Inserter.x[i] >= -MYEPSILON) && ((Inserter.x[i]-1.) <= MYEPSILON);
     908          Inserter.MatrixMultiplication(M);
     909
     910          // Check whether point is in- or outside
     911          for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
     912            // get linked cell list
     913            if (TesselStruct[(*ListRunner)] != NULL) {
     914              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
     915              FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
     916            }
     917          }
     918          // insert into Filling
     919          if (FillIt) {
     920            DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl);
     921            // copy atom ...
     922            CopyAtoms[Walker->nr] = Walker->clone();
     923            CopyAtoms[Walker->nr]->x.CopyVector(&Inserter);
     924            Filling->AddAtom(CopyAtoms[Walker->nr]);
     925            DoLog(4) && (Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl);
    863926          } else {
    864             const double distance = (TesselStruct[i]->GetDistanceSquaredToSurface(CurrentPosition, LCList[i]));
    865             FillIt = FillIt && (distance > boundary*boundary);
    866             if (FillIt) {
    867               Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is outer point." << endl;
    868             } else {
    869               Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point or within boundary." << endl;
    870               break;
    871             }
    872             i++;
     927            DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl);
     928            CopyAtoms[Walker->nr] = NULL;
     929            continue;
    873930          }
    874931        }
    875 
    876         if (FillIt) {
    877           // fill in Filler
    878           Log() << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
    879 
    880           // create molecule random translation vector ...
    881           for (int i=0;i<NDIM;i++)
    882             FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
    883           Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
    884 
    885           // go through all atoms
    886           Walker = filler->start;
    887           while (Walker->next != filler->end) {
    888             Walker = Walker->next;
    889             // copy atom ...
    890             CopyAtoms[Walker->nr] = Walker->clone();
    891 
    892             // create atomic random translation vector ...
    893             for (int i=0;i<NDIM;i++)
    894               AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
    895 
    896             // ... and rotation matrix
    897             if (DoRandomRotation) {
    898               for (int i=0;i<NDIM;i++) {
    899                 phi[i] = rand()/(RAND_MAX/(2.*M_PI));
    900               }
    901 
    902               Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
    903               Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
    904               Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
    905               Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
    906               Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
    907               Rotations[7] =               sin(phi[1])                                                ;
    908               Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
    909               Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
    910               Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
    911             }
    912 
    913             // ... and put at new position
    914             if (DoRandomRotation)
    915               CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
    916             CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
    917             CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
    918             CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
    919 
    920             // insert into Filling
    921 
    922             // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
    923             Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
    924             Filling->AddAtom(CopyAtoms[Walker->nr]);
    925           }
    926 
    927           // go through all bonds and add as well
    928           Binder = filler->first;
    929           while(Binder->next != filler->last) {
    930             Binder = Binder->next;
    931             Log() << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
     932        // go through all bonds and add as well
     933        Binder = filler->first;
     934        while(Binder->next != filler->last) {
     935          Binder = Binder->next;
     936          if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
     937            Log()  << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
    932938            Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
    933939          }
    934         } else {
    935           // leave empty
    936           Log() << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
    937940        }
    938941      }
    939942  Free(&M);
    940 
    941   // output to file
    942   TesselStruct[0]->LastTriangle = NULL;
    943   StoreTrianglesinFile(Filling, TesselStruct[0], "Tesselated", ".dat");
    944 
    945   for (size_t i=0;i<List->ListOfMolecules.size();i++) {
    946         delete(LCList[i]);
    947         delete(TesselStruct[i]);
    948   }
     943  Free(&MInverse);
     944
    949945  return Filling;
    950946};
    … …  
    965961  bool freeLC = false;
    966962  bool status = false;
    967   CandidateForTesselation *baseline=0;
    968   LineMap::iterator testline;
     963  CandidateForTesselation *baseline = NULL;
    969964  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
    970965  bool TesselationFailFlag = false;
    971   BoundaryTriangleSet *T = NULL;
    972966
    973967  if (TesselStruct == NULL) {
    974     Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
     968    DoLog(1) && (Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl);
    975969    TesselStruct= new Tesselation;
    976970  } else {
    977971    delete(TesselStruct);
    978     Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
     972    DoLog(1) && (Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl);
    979973    TesselStruct = new Tesselation;
    980974  }
    … …  
    987981
    988982  // 1. get starting triangle
    989   TesselStruct->FindStartingTriangle(RADIUS, LCList);
     983  if (!TesselStruct->FindStartingTriangle(RADIUS, LCList)) {
     984    DoeLog(0) && (eLog() << Verbose(0) << "No valid starting triangle found." << endl);
     985    //performCriticalExit();
     986  }
     987  if (filename != NULL) {
     988    if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
     989      TesselStruct->Output(filename, mol);
     990    }
     991  }
    990992
    991993  // 2. expand from there
    992994  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
    993     // 2a. fill all new OpenLines
    994     Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
     995    (cerr << "There are " <<  TesselStruct->TrianglesOnBoundary.size() << " triangles and " << TesselStruct->OpenLines.size() << " open lines to scan for candidates." << endl);
     996    // 2a. print OpenLines without candidates
     997    DoLog(1) && (Log() << Verbose(1) << "There are the following open lines to scan for a candidates:" << endl);
    995998    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
    996       Log() << Verbose(2) << *(Runner->second) << endl;
    997 
    998     for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
    999       baseline = Runner->second;
    1000       if (baseline->pointlist.empty()) {
    1001         T = (((baseline->BaseLine->triangles.begin()))->second);
    1002         Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
    1003         TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
    1004       }
    1005     }
    1006 
    1007     // 2b. search for smallest ShortestAngle among all candidates
     999      if (Runner->second->pointlist.empty())
     1000        DoLog(1) && (Log() << Verbose(1) << " " << *(Runner->second) << endl);
     1001
     1002    // 2b. find best candidate for each OpenLine
     1003    TesselationFailFlag = TesselStruct->FindCandidatesforOpenLines(RADIUS, LCList);
     1004
     1005    // 2c. print OpenLines with candidates again
     1006    DoLog(1) && (Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl);
     1007    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
     1008      DoLog(1) && (Log() << Verbose(1) << " " << *(Runner->second) << endl);
     1009
     1010    // 2d. search for smallest ShortestAngle among all candidates
    10081011    double ShortestAngle = 4.*M_PI;
    1009     Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
    1010     for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
    1011       Log() << Verbose(2) << *(Runner->second) << endl;
    1012 
    10131012    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
    10141013      if (Runner->second->ShortestAngle < ShortestAngle) {
    10151014        baseline = Runner->second;
    10161015        ShortestAngle = baseline->ShortestAngle;
    1017         //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
     1016        DoLog(1) && (Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl);
    10181017      }
    10191018    }
     1019    // 2e. if we found one, add candidate
    10201020    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
    10211021      OneLoopWithoutSuccessFlag = false;
    10221022    else {
    1023       TesselStruct->AddCandidateTriangle(*baseline);
    1024     }
    1025 
    1026     // write temporary envelope
     1023      TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
     1024    }
     1025
     1026    // 2f. write temporary envelope
    10271027    if (filename != NULL) {
    10281028      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
    … …  
    10591059  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
    10601060
    1061   // correct degenerated polygons
    1062   TesselStruct->CorrectAllDegeneratedPolygons();
    1063 
    1064   // check envelope for consistency
    1065   status = CheckListOfBaselines(TesselStruct);
     1061//  // correct degenerated polygons
     1062//  TesselStruct->CorrectAllDegeneratedPolygons();
     1063//
     1064//  // check envelope for consistency
     1065//  status = CheckListOfBaselines(TesselStruct);
    10661066
    10671067  // write final envelope
Note: See TracChangeset for help on using the changeset viewer.