Changes in src/analysis_correlation.cpp [a67d19:1614174]

File:

• src/analysis_correlation.cpp (modified) (24 diffs)

src/analysis_correlation.cpp

@@ -10 +10 @@
 #include "analysis_correlation.hpp"
 #include "element.hpp"
-#include "info.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
-#include "triangleintersectionlist.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
-#include "World.hpp"
 
 
@@ -31 +28 @@
 PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
 {
-  Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
 
   if (molecules->ListOfMolecules.empty()) {
-    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
+    eLog() << Verbose(1) <<"No molecule given." << endl;
     return outmap;
   }
@@ -42 +38 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
+      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type1 == NULL) || (Walker->type == type1)) {
           for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
             if ((*MolOtherWalker)->ActiveFlag) {
-              DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
+              Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
               atom *OtherWalker = (*MolOtherWalker)->start;
               while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
                 OtherWalker = OtherWalker->next;
-                DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
+                Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
                 if (Walker->nr < OtherWalker->nr)
                   if ((type2 == NULL) || (OtherWalker->type == type2)) {
-                    distance = Walker->node->PeriodicDistance(OtherWalker->node, World::get()->cell_size);
+                    distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
                     //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
@@ -81 +77 @@
 PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
 {
-  Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
@@ -92 +87 @@
 
   if (molecules->ListOfMolecules.empty()) {
-    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
+    eLog() << Verbose(1) <<"No molecule given." << endl;
     return outmap;
   }
@@ -98 +93 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
+      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
-      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
+      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type1 == NULL) || (Walker->type == type1)) {
           periodicX.CopyVector(Walker->node);
@@ -117 +112 @@
                 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
                   if ((*MolOtherWalker)->ActiveFlag) {
-                    DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
+                    Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
                     atom *OtherWalker = (*MolOtherWalker)->start;
                     while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
                       OtherWalker = OtherWalker->next;
-                      DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
+                      Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
                       if (Walker->nr < OtherWalker->nr)
                         if ((type2 == NULL) || (OtherWalker->type == type2)) {
@@ -159 +154 @@
 CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
 {
-  Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
 
   if (molecules->ListOfMolecules.empty()) {
-    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
+    Log() << Verbose(1) <<"No molecule given." << endl;
     return outmap;
   }
@@ -170 +164 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
+      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          distance = Walker->node->PeriodicDistance(point, World::get()->cell_size);
-          DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
+          distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
+          Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
         }
@@ -196 +190 @@
 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
 {
-  Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
@@ -204 +197 @@
 
   if (molecules->ListOfMolecules.empty()) {
-    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
+    Log() << Verbose(1) <<"No molecule given." << endl;
     return outmap;
   }
@@ -210 +203 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
+      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
-      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
+      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
           periodicX.CopyVector(Walker->node);
@@ -228 +221 @@
                 checkX.MatrixMultiplication(FullMatrix);
                 distance = checkX.Distance(point);
-                DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
+                Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
                 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
               }
@@ -250 +243 @@
 CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
 {
-  Info FunctionInfo(__func__);
   CorrelationToSurfaceMap *outmap = NULL;
   double distance = 0;
@@ -257 +249 @@
 
   if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
-    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
+    Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
     return outmap;
   }
@@ -263 +255 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        //Log() << Verbose(1) << "Current atom is " << *Walker << "." << endl;
+      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          TriangleIntersectionList Intersections(Walker->node,Surface,LC);
-          distance = Intersections.GetSmallestDistance();
-          triangle = Intersections.GetClosestTriangle();
-          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
-        }
-      }
-    } else
-      DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
-
+          triangle = Surface->FindClosestTriangleToPoint(Walker->node, LC );
+          if (triangle != NULL) {
+            distance = DistanceToTrianglePlane(Walker->node, triangle);
+            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
+          }
+        }
+      }
+    }
 
   return outmap;
@@ -297 +288 @@
 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
 {
-  Info FunctionInfo(__func__);
   CorrelationToSurfaceMap *outmap = NULL;
   double distance = 0;
@@ -307 +297 @@
 
   if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
-    DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
+    Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
     return outmap;
   }
   outmap = new CorrelationToSurfaceMap;
-  double ShortestDistance = 0.;
-  BoundaryTriangleSet *ShortestTriangle = NULL;
-  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
-    if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
+    if ((*MolWalker)->ActiveFlag) {
+      double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
-      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
-      atom *Walker = (*MolWalker)->start;
-      while (Walker->next != (*MolWalker)->end) {
-        Walker = Walker->next;
-        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
+      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
           periodicX.CopyVector(Walker->node);
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
-          ShortestDistance = -1.;
           for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
@@ -333 +320 @@
                 checkX.AddVector(&periodicX);
                 checkX.MatrixMultiplication(FullMatrix);
-                TriangleIntersectionList Intersections(&checkX,Surface,LC);
-                distance = Intersections.GetSmallestDistance();
-                triangle = Intersections.GetClosestTriangle();
-                if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
-                  ShortestDistance = distance;
-                  ShortestTriangle = triangle;
+                triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
+                if (triangle != NULL) {
+                  distance = DistanceToTrianglePlane(&checkX, triangle);
+                  outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
                 }
-              }
-          // insert
-          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
-          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
+          }
         }
       }
@@ -353 +335 @@
 };
 
-/** Returns the index of the bin for a given value.
+/** Returns the start of the bin for a given value.
  * \param value value whose bin to look for
  * \param BinWidth width of bin
  * \param BinStart first bin
  */
-int GetBin ( const double value, const double BinWidth, const double BinStart )
-{
-  Info FunctionInfo(__func__);
-  int bin =(int) (floor((value - BinStart)/BinWidth));
-  return (bin);
+double GetBin ( const double value, const double BinWidth, const double BinStart )
+{
+  double bin =(double) (floor((value - BinStart)/BinWidth));
+  return (bin*BinWidth+BinStart);
 };
 
@@ -372 +353 @@
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
 {
-  Info FunctionInfo(__func__);
-  *file << "BinStart\tCount" << endl;
+  *file << "# BinStart\tCount" << endl;
   for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
+    *file << runner->first << "\t" << runner->second << endl;
   }
 };
@@ -385 +365 @@
 void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
 {
-  Info FunctionInfo(__func__);
-  *file << "BinStart\tAtom1\tAtom2" << endl;
+  *file << "# BinStart\tAtom1\tAtom2" << endl;
   for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
+    *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
   }
 };
@@ -398 +377 @@
 void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
 {
-  Info FunctionInfo(__func__);
-  *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
+  *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
   for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     *file << runner->first;
     for (int i=0;i<NDIM;i++)
-      *file << "\t" << setprecision(8) << (runner->second.first->node->x[i] - runner->second.second->x[i]);
+      *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
     *file << endl;
   }
@@ -414 +392 @@
 void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
 {
-  Info FunctionInfo(__func__);
-  *file << "BinStart\tTriangle" << endl;
-  if (!map->empty())
-    for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
-      *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
-    }
-};
-
+  *file << "# BinStart\tTriangle" << endl;
+  for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
+    *file << runner->first << "\t" << *(runner->second.second) << endl;
+  }
+};
+