Changes in src/linkedcell.cpp [caf4ba:717e0c]

File:

• src/linkedcell.cpp (modified) (14 diffs)

src/linkedcell.cpp

@@ +1 @@
+/** \file linkedcell.cpp
+ *
+ * Function implementations for the class LinkedCell.
+ *
+ */
+
+
+#include "atom.hpp"
+#include "helpers.hpp"
 #include "linkedcell.hpp"
-#include "molecules.hpp"
+#include "log.hpp"
+#include "molecule.hpp"
 #include "tesselation.hpp"
+#include "vector.hpp"
 
 // ========================================================= class LinkedCell ===========================================
@@ -23 +34 @@
  * \param RADIUS edge length of cells
  */
-LinkedCell::LinkedCell(PointCloud *set, double radius)
+LinkedCell::LinkedCell(const PointCloud * const set, const double radius)
 {
   TesselPoint *Walker = NULL;
@@ -34 +45 @@
   max.Zero();
   min.Zero();
-  cout << Verbose(1) << "Begin of LinkedCell" << endl;
-  if ((set == NULL) || (set->IsEmpty())) {
-    cerr << "ERROR: set is NULL or contains no linked cell nodes!" << endl;
+  Log() << Verbose(1) << "Begin of LinkedCell" << endl;
+  if (set->IsEmpty()) {
+    eLog() << Verbose(1) << "set contains no linked cell nodes!" << endl;
     return;
   }
@@ -47 +58 @@
   }
   set->GoToFirst();
-  while (!set->IsLast()) {
+  while (!set->IsEnd()) {
     Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
@@ -57 +68 @@
     set->GoToNext();
   }
-  cout << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
+  Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
 
   // 2. find then number of cells per axis
@@ -63 +74 @@
     N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
   }
-  cout << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl;
+  Log() << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl;
 
   // 3. allocate the lists
-  cout << Verbose(2) << "Allocating cells ... ";
+  Log() << Verbose(2) << "Allocating cells ... ";
   if (LC != NULL) {
-    cout << Verbose(1) << "ERROR: Linked Cell list is already allocated, I do nothing." << endl;
+    eLog() << Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl;
     return;
   }
@@ -75 +86 @@
     LC [index].clear();
   }
-  cout << "done."  << endl;
+  Log() << Verbose(0) << "done."  << endl;
 
   // 4. put each atom into its respective cell
-  cout << Verbose(2) << "Filling cells ... ";
+  Log() << Verbose(2) << "Filling cells ... ";
   set->GoToFirst();
-  while (!set->IsLast()) {
+  while (!set->IsEnd()) {
     Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
@@ -87 +98 @@
     index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
     LC[index].push_back(Walker);
-    //cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
+    //Log() << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
     set->GoToNext();
   }
-  cout << "done."  << endl;
-  cout << Verbose(1) << "End of LinkedCell" << endl;
+  Log() << Verbose(0) << "done."  << endl;
+  Log() << Verbose(1) << "End of LinkedCell" << endl;
+};
+
+
+/** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
+ * \param *set LCNodeSet class with all LCNode's
+ * \param RADIUS edge length of cells
+ */
+LinkedCell::LinkedCell(LinkedNodes *set, const double radius)
+{
+  class TesselPoint *Walker = NULL;
+  RADIUS = radius;
+  LC = NULL;
+  for(int i=0;i<NDIM;i++)
+    N[i] = 0;
+  index = -1;
+  max.Zero();
+  min.Zero();
+  Log() << Verbose(1) << "Begin of LinkedCell" << endl;
+  if (set->empty()) {
+    eLog() << Verbose(1) << "set contains no linked cell nodes!" << endl;
+    return;
+  }
+  // 1. find max and min per axis of atoms
+  LinkedNodes::iterator Runner = set->begin();
+  for (int i=0;i<NDIM;i++) {
+    max.x[i] = (*Runner)->node->x[i];
+    min.x[i] = (*Runner)->node->x[i];
+  }
+  for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
+    Walker = *Runner;
+    for (int i=0;i<NDIM;i++) {
+      if (max.x[i] < Walker->node->x[i])
+        max.x[i] = Walker->node->x[i];
+      if (min.x[i] > Walker->node->x[i])
+        min.x[i] = Walker->node->x[i];
+    }
+  }
+  Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
+
+  // 2. find then number of cells per axis
+  for (int i=0;i<NDIM;i++) {
+    N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
+  }
+  Log() << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl;
+
+  // 3. allocate the lists
+  Log() << Verbose(2) << "Allocating cells ... ";
+  if (LC != NULL) {
+    eLog() << Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl;
+    return;
+  }
+  LC = new LinkedNodes[N[0]*N[1]*N[2]];
+  for (index=0;index<N[0]*N[1]*N[2];index++) {
+    LC [index].clear();
+  }
+  Log() << Verbose(0) << "done."  << endl;
+
+  // 4. put each atom into its respective cell
+  Log() << Verbose(2) << "Filling cells ... ";
+  for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
+    Walker = *Runner;
+    for (int i=0;i<NDIM;i++) {
+      n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
+    }
+    index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
+    LC[index].push_back(Walker);
+    //Log() << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
+  }
+  Log() << Verbose(0) << "done."  << endl;
+  Log() << Verbose(1) << "End of LinkedCell" << endl;
 };
 
@@ -112 +193 @@
  * \return if all in intervals - true, else -false
  */
-bool LinkedCell::CheckBounds()
+bool LinkedCell::CheckBounds() const
 {
   bool status = true;
@@ -118 +199 @@
     status = status && ((n[i] >=0) && (n[i] < N[i]));
   if (!status)
-  cerr << "ERROR: indices are out of bounds!" << endl;
+  eLog() << Verbose(1) << "indices are out of bounds!" << endl;
+  return status;
+};
+
+/** Checks whether LinkedCell::n[] plus relative offset is each within [0,N[]].
+ * Note that for this check we don't admonish if out of bounds.
+ * \param relative[NDIM] relative offset to current cell
+ * \return if all in intervals - true, else -false
+ */
+bool LinkedCell::CheckBounds(const int relative[NDIM]) const
+{
+  bool status = true;
+  for(int i=0;i<NDIM;i++)
+    status = status && ((n[i]+relative[i] >=0) && (n[i]+relative[i] < N[i]));
   return status;
 };
@@ -126 +220 @@
  * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
  */
-LinkedNodes* LinkedCell::GetCurrentCell()
+const LinkedNodes* LinkedCell::GetCurrentCell() const
 {
   if (CheckBounds()) {
@@ -136 +230 @@
 };
 
+/** Returns a pointer to the current cell.
+ * \param relative[NDIM] offset for each axis with respect to the current cell LinkedCell::n[NDIM]
+ * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[]+relative[] are out of bounds.
+ */
+const LinkedNodes* LinkedCell::GetRelativeToCurrentCell(const int relative[NDIM]) const
+{
+  if (CheckBounds(relative)) {
+    index = (n[0]+relative[0]) * N[1] * N[2] + (n[1]+relative[1]) * N[2] + (n[2]+relative[2]);
+    return (&(LC[index]));
+  } else {
+    return NULL;
+  }
+};
+
 /** Calculates the index for a given LCNode *Walker.
  * \param *Walker LCNode to set index tos
  * \return if the atom is also found in this cell - true, else - false
  */
-bool LinkedCell::SetIndexToNode(TesselPoint *Walker)
+bool LinkedCell::SetIndexToNode(const TesselPoint * const Walker) const
 {
   bool status = false;
@@ -152 +260 @@
     return status;
   } else {
-    cerr << Verbose(1) << "ERROR: Node at " << *Walker << " is out of bounds." << endl;
+    eLog() << Verbose(1) << "Node at " << *Walker << " is out of bounds." << endl;
     return false;
+  }
+};
+
+/** Calculates the interval bounds of the linked cell grid.
+ * \param *lower lower bounds
+ * \param *upper upper bounds
+ */
+void LinkedCell::GetNeighbourBounds(int lower[NDIM], int upper[NDIM]) const
+{
+  for (int i=0;i<NDIM;i++) {
+    lower[i] = ((n[i]-1) >= 0) ? n[i]-1 : 0;
+    upper[i] = ((n[i]+1) < N[i]) ? n[i]+1 : N[i]-1;
+    //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+    // check for this axis whether the point is outside of our grid
+    if (n[i] < 0)
+      upper[i] = lower[i];
+    if (n[i] > N[i])
+      lower[i] = upper[i];
+
+    //Log() << Verbose(0) << "axis " << i << " has bounds [" << lower[i] << "," << upper[i] << "]" << endl;
   }
 };
@@ -161 +289 @@
  * \return Vector is inside bounding box - true, else - false
  */
-bool LinkedCell::SetIndexToVector(Vector *x)
+bool LinkedCell::SetIndexToVector(const Vector * const x) const
 {
   bool status = true;