Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/analysis_correlation.cpp

    rbd61b41 ra67d19  
    1515#include "tesselation.hpp"
    1616#include "tesselationhelpers.hpp"
     17#include "triangleintersectionlist.hpp"
    1718#include "vector.hpp"
    1819#include "verbose.hpp"
     20#include "World.hpp"
    1921
    2022
    … …  
    3436
    3537  if (molecules->ListOfMolecules.empty()) {
    36     eLog() << Verbose(1) <<"No molecule given." << endl;
     38    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    3739    return outmap;
    3840  }
    … …  
    4042  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    4143    if ((*MolWalker)->ActiveFlag) {
    42       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    43       atom *Walker = (*MolWalker)->start;
    44       while (Walker->next != (*MolWalker)->end) {
    45         Walker = Walker->next;
    46         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
     44      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     45      atom *Walker = (*MolWalker)->start;
     46      while (Walker->next != (*MolWalker)->end) {
     47        Walker = Walker->next;
     48        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
    4749        if ((type1 == NULL) || (Walker->type == type1)) {
    4850          for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
    4951            if ((*MolOtherWalker)->ActiveFlag) {
    50               Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
     52              DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    5153              atom *OtherWalker = (*MolOtherWalker)->start;
    5254              while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
    5355                OtherWalker = OtherWalker->next;
    54                 Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
     56                DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
    5557                if (Walker->nr < OtherWalker->nr)
    5658                  if ((type2 == NULL) || (OtherWalker->type == type2)) {
    57                     distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
     59                    distance = Walker->node->PeriodicDistance(OtherWalker->node, World::get()->cell_size);
    5860                    //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
    5961                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
    … …  
    9092
    9193  if (molecules->ListOfMolecules.empty()) {
    92     eLog() << Verbose(1) <<"No molecule given." << endl;
     94    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    9395    return outmap;
    9496  }
    … …  
    9698  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    9799    if ((*MolWalker)->ActiveFlag) {
    98       double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
     100      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
    99101      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    100       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    101       atom *Walker = (*MolWalker)->start;
    102       while (Walker->next != (*MolWalker)->end) {
    103         Walker = Walker->next;
    104         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
     102      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     103      atom *Walker = (*MolWalker)->start;
     104      while (Walker->next != (*MolWalker)->end) {
     105        Walker = Walker->next;
     106        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
    105107        if ((type1 == NULL) || (Walker->type == type1)) {
    106108          periodicX.CopyVector(Walker->node);
    … …  
    115117                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
    116118                  if ((*MolOtherWalker)->ActiveFlag) {
    117                     Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
     119                    DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    118120                    atom *OtherWalker = (*MolOtherWalker)->start;
    119121                    while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
    120122                      OtherWalker = OtherWalker->next;
    121                       Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
     123                      DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
    122124                      if (Walker->nr < OtherWalker->nr)
    123125                        if ((type2 == NULL) || (OtherWalker->type == type2)) {
    … …  
    162164
    163165  if (molecules->ListOfMolecules.empty()) {
    164     Log() << Verbose(1) <<"No molecule given." << endl;
     166    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    165167    return outmap;
    166168  }
    … …  
    168170  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    169171    if ((*MolWalker)->ActiveFlag) {
    170       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    171       atom *Walker = (*MolWalker)->start;
    172       while (Walker->next != (*MolWalker)->end) {
    173         Walker = Walker->next;
    174         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
     172      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     173      atom *Walker = (*MolWalker)->start;
     174      while (Walker->next != (*MolWalker)->end) {
     175        Walker = Walker->next;
     176        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
    175177        if ((type == NULL) || (Walker->type == type)) {
    176           distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
    177           Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
     178          distance = Walker->node->PeriodicDistance(point, World::get()->cell_size);
     179          DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    178180          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
    179181        }
    … …  
    202204
    203205  if (molecules->ListOfMolecules.empty()) {
    204     Log() << Verbose(1) <<"No molecule given." << endl;
     206    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    205207    return outmap;
    206208  }
    … …  
    208210  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    209211    if ((*MolWalker)->ActiveFlag) {
    210       double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
     212      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
    211213      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    212       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    213       atom *Walker = (*MolWalker)->start;
    214       while (Walker->next != (*MolWalker)->end) {
    215         Walker = Walker->next;
    216         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
     214      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     215      atom *Walker = (*MolWalker)->start;
     216      while (Walker->next != (*MolWalker)->end) {
     217        Walker = Walker->next;
     218        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
    217219        if ((type == NULL) || (Walker->type == type)) {
    218220          periodicX.CopyVector(Walker->node);
    … …  
    226228                checkX.MatrixMultiplication(FullMatrix);
    227229                distance = checkX.Distance(point);
    228                 Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
     230                DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    229231                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
    230232              }
    … …  
    255257
    256258  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    257     Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
     259    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    258260    return outmap;
    259261  }
    … …  
    261263  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    262264    if ((*MolWalker)->ActiveFlag) {
    263       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    264       atom *Walker = (*MolWalker)->start;
    265       while (Walker->next != (*MolWalker)->end) {
    266         Walker = Walker->next;
    267         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
     265      DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
     266      atom *Walker = (*MolWalker)->start;
     267      while (Walker->next != (*MolWalker)->end) {
     268        Walker = Walker->next;
     269        //Log() << Verbose(1) << "Current atom is " << *Walker << "." << endl;
    268270        if ((type == NULL) || (Walker->type == type)) {
    269           triangle = Surface->FindClosestTriangleToPoint(Walker->node, LC );
    270           if (triangle != NULL) {
    271             distance = DistanceToTrianglePlane(Walker->node, triangle);
    272             outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
    273           }
    274         }
    275       }
    276     }
     271          TriangleIntersectionList Intersections(Walker->node,Surface,LC);
     272          distance = Intersections.GetSmallestDistance();
     273          triangle = Intersections.GetClosestTriangle();
     274          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
     275        }
     276      }
     277    } else
     278      DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
     279
    277280
    278281  return outmap;
    … …  
    304307
    305308  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    306     Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
     309    DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    307310    return outmap;
    308311  }
    309312  outmap = new CorrelationToSurfaceMap;
    310   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    311     if ((*MolWalker)->ActiveFlag) {
    312       double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
     313  double ShortestDistance = 0.;
     314  BoundaryTriangleSet *ShortestTriangle = NULL;
     315  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     316    if ((*MolWalker)->ActiveFlag) {
     317      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
    313318      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    314       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    315       atom *Walker = (*MolWalker)->start;
    316       while (Walker->next != (*MolWalker)->end) {
    317         Walker = Walker->next;
    318         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
     319      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
     320      atom *Walker = (*MolWalker)->start;
     321      while (Walker->next != (*MolWalker)->end) {
     322        Walker = Walker->next;
     323        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
    319324        if ((type == NULL) || (Walker->type == type)) {
    320325          periodicX.CopyVector(Walker->node);
    321326          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    322327          // go through every range in xyz and get distance
     328          ShortestDistance = -1.;
    323329          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    324330            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    … …  
    327333                checkX.AddVector(&periodicX);
    328334                checkX.MatrixMultiplication(FullMatrix);
    329                 triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
    330                 if (triangle != NULL) {
    331                   distance = DistanceToTrianglePlane(&checkX, triangle);
    332                   outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
     335                TriangleIntersectionList Intersections(&checkX,Surface,LC);
     336                distance = Intersections.GetSmallestDistance();
     337                triangle = Intersections.GetClosestTriangle();
     338                if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
     339                  ShortestDistance = distance;
     340                  ShortestTriangle = triangle;
    333341                }
    334           }
     342              }
     343          // insert
     344          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
     345          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
    335346        }
    336347      }
    … …  
    362373{
    363374  Info FunctionInfo(__func__);
    364   *file << "# BinStart\tCount" << endl;
     375  *file << "BinStart\tCount" << endl;
    365376  for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    366     *file << runner->first << "\t" << runner->second << endl;
     377    *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
    367378  }
    368379};
    … …  
    375386{
    376387  Info FunctionInfo(__func__);
    377   *file << "# BinStart\tAtom1\tAtom2" << endl;
     388  *file << "BinStart\tAtom1\tAtom2" << endl;
    378389  for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    379     *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
     390    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
    380391  }
    381392};
    … …  
    388399{
    389400  Info FunctionInfo(__func__);
    390   *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
     401  *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
    391402  for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    392403    *file << runner->first;
    393404    for (int i=0;i<NDIM;i++)
    394       *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
     405      *file << "\t" << setprecision(8) << (runner->second.first->node->x[i] - runner->second.second->x[i]);
    395406    *file << endl;
    396407  }
    … …  
    404415{
    405416  Info FunctionInfo(__func__);
    406   *file << "# BinStart\tTriangle" << endl;
    407   for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    408     *file << runner->first << "\t" << *(runner->second.second) << endl;
    409   }
    410 };
    411 
     417  *file << "BinStart\tTriangle" << endl;
     418  if (!map->empty())
     419    for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     420      *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
     421    }
     422};
     423
Note: See TracChangeset for help on using the changeset viewer.