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

File:

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

src/linkedcell.cpp

@@ -45 +45 @@
   max.Zero();
   min.Zero();
-  Log() << Verbose(1) << "Begin of LinkedCell" << endl;
-  if (set->IsEmpty()) {
-    eLog() << Verbose(1) << "set contains no linked cell nodes!" << endl;
+  DoLog(1) && (Log() << Verbose(1) << "Begin of LinkedCell" << endl);
+  if ((set == NULL) || (set->IsEmpty())) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "set is NULL or contains no linked cell nodes!" << endl);
     return;
   }
@@ -68 +68 @@
     set->GoToNext();
   }
-  Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
+  DoLog(2) && (Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl);
 
   // 2. find then number of cells per axis
@@ -74 +74 @@
     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;
+  DoLog(2) && (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 ... ";
+  DoLog(2) && (Log() << Verbose(2) << "Allocating cells ... ");
   if (LC != NULL) {
-    eLog() << Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl);
     return;
   }
@@ -86 +86 @@
     LC [index].clear();
   }
-  Log() << Verbose(0) << "done."  << endl;
+  DoLog(0) && (Log() << Verbose(0) << "done."  << endl);
 
   // 4. put each atom into its respective cell
-  Log() << Verbose(2) << "Filling cells ... ";
+  DoLog(2) && (Log() << Verbose(2) << "Filling cells ... ");
   set->GoToFirst();
   while (!set->IsEnd()) {
@@ -101 +101 @@
     set->GoToNext();
   }
-  Log() << Verbose(0) << "done."  << endl;
-  Log() << Verbose(1) << "End of LinkedCell" << endl;
+  DoLog(0) && (Log() << Verbose(0) << "done."  << endl);
+  DoLog(1) && (Log() << Verbose(1) << "End of LinkedCell" << endl);
 };
 
@@ -120 +120 @@
   max.Zero();
   min.Zero();
-  Log() << Verbose(1) << "Begin of LinkedCell" << endl;
+  DoLog(1) && (Log() << Verbose(1) << "Begin of LinkedCell" << endl);
   if (set->empty()) {
-    eLog() << Verbose(1) << "set contains no linked cell nodes!" << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "set contains no linked cell nodes!" << endl);
     return;
   }
@@ -140 +140 @@
     }
   }
-  Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
+  DoLog(2) && (Log() << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl);
 
   // 2. find then number of cells per axis
@@ -146 +146 @@
     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;
+  DoLog(2) && (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 ... ";
+  DoLog(2) && (Log() << Verbose(2) << "Allocating cells ... ");
   if (LC != NULL) {
-    eLog() << Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl);
     return;
   }
@@ -158 +158 @@
     LC [index].clear();
   }
-  Log() << Verbose(0) << "done."  << endl;
+  DoLog(0) && (Log() << Verbose(0) << "done."  << endl);
 
   // 4. put each atom into its respective cell
-  Log() << Verbose(2) << "Filling cells ... ";
+  DoLog(2) && (Log() << Verbose(2) << "Filling cells ... ");
   for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
     Walker = *Runner;
@@ -171 +171 @@
     //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;
+  DoLog(0) && (Log() << Verbose(0) << "done."  << endl);
+  DoLog(1) && (Log() << Verbose(1) << "End of LinkedCell" << endl);
 };
 
@@ -199 +199 @@
     status = status && ((n[i] >=0) && (n[i] < N[i]));
   if (!status)
-  eLog() << Verbose(1) << "indices are out of bounds!" << endl;
+  DoeLog(1) && (eLog()<< Verbose(1) << "indices are out of bounds!" << endl);
   return status;
 };
@@ -220 +220 @@
  * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
  */
-const LinkedNodes* LinkedCell::GetCurrentCell() const
+const LinkedCell::LinkedNodes* LinkedCell::GetCurrentCell() const
 {
   if (CheckBounds()) {
@@ -234 +234 @@
  * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[]+relative[] are out of bounds.
  */
-const LinkedNodes* LinkedCell::GetRelativeToCurrentCell(const int relative[NDIM]) const
+const LinkedCell::LinkedNodes* LinkedCell::GetRelativeToCurrentCell(const int relative[NDIM]) const
 {
   if (CheckBounds(relative)) {
@@ -242 +242 @@
     return NULL;
   }
+};
+
+/** Set the index to the cell containing a given Vector *x.
+ * \param *x Vector with coordinates
+ * \return Vector is inside bounding box - true, else - false
+ */
+bool LinkedCell::SetIndexToVector(const Vector * const x) const
+{
+  for (int i=0;i<NDIM;i++)
+    n[i] = (int)floor((x->x[i] - min.x[i])/RADIUS);
+
+  return CheckBounds();
 };
 
@@ -260 +272 @@
     return status;
   } else {
-    eLog() << Verbose(1) << "Node at " << *Walker << " is out of bounds." << endl;
+    DoeLog(1) && (eLog()<< Verbose(1) << "Node at " << *Walker << " is out of bounds." << endl);
     return false;
   }
@@ -268 +280 @@
  * \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
+ * \param step how deep to check the neighbouring cells (i.e. number of layers to check)
+ */
+void LinkedCell::GetNeighbourBounds(int lower[NDIM], int upper[NDIM], int step) const
+{
+  for (int i=0;i<NDIM;i++) {
+    lower[i] = n[i];
+    for (int s=step; s>0;--s)
+      if ((n[i]-s) >= 0) {
+        lower[i] = n[i]-s;
+        break;
+      }
+    upper[i] = n[i];
+    for (int s=step; s>0;--s)
+      if ((n[i]+s) < N[i]) {
+        upper[i] = n[i]+s;
+        break;
+      }
+    //Log() << Verbose(0) << "axis " << i << " has bounds [" << lower[i] << "," << upper[i] << "]" << endl;
+  }
+};
+
+/** Returns a list with all neighbours from the current LinkedCell::index.
+ * \param distance (if no distance, then adjacent cells are taken)
+ * \return list of tesselpoints
+ */
+LinkedCell::LinkedNodes* LinkedCell::GetallNeighbours(const double distance) const
+{
+  int Nlower[NDIM], Nupper[NDIM];
+  TesselPoint *Walker = NULL;
+  LinkedNodes *TesselList = new LinkedNodes;
+
+  // then go through the current and all neighbouring cells and check the contained points for possible candidates
+  const int step = (distance == 0) ? 1 : (int)floor(distance/RADIUS + 1.);
+  GetNeighbourBounds(Nlower, Nupper, step);
+
+  //Log() << Verbose(0) << endl;
+  for (n[0] = Nlower[0]; n[0] <= Nupper[0]; n[0]++)
+    for (n[1] = Nlower[1]; n[1] <= Nupper[1]; n[1]++)
+      for (n[2] = Nlower[2]; n[2] <= Nupper[2]; n[2]++) {
+        const LinkedNodes *List = GetCurrentCell();
+        //Log() << Verbose(1) << "Current cell is " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
+        if (List != NULL) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            Walker = *Runner;
+            TesselList->push_back(Walker);
+          }
+        }
+      }
+  return TesselList;
+};
+
+/** Set the index to the cell containing a given Vector *x, which is not inside the LinkedCell's domain
+ * Note that as we have to check distance from every corner of the closest cell, this function is faw more
+ * expensive and if Vector is known to be inside LinkedCell's domain, then SetIndexToVector() should be used.
+ * \param *x Vector with coordinates
+ * \return minimum squared distance of cell to given vector (if inside of domain, distance is 0)
+ */
+double LinkedCell::SetClosestIndexToOutsideVector(const Vector * const x) const
+{
+  for (int i=0;i<NDIM;i++) {
+    n[i] = (int)floor((x->x[i] - min.x[i])/RADIUS);
     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;
-  }
-};
-
-/** Calculates the index for a given Vector *x.
- * \param *x Vector with coordinates
- * \return Vector is inside bounding box - true, else - false
- */
-bool LinkedCell::SetIndexToVector(const Vector * const x) const
-{
-  bool status = true;
-  for (int i=0;i<NDIM;i++) {
-    n[i] = (int)floor((x->x[i] - min.x[i])/RADIUS);
-    if (max.x[i] < x->x[i])
-      status = false;
-    if (min.x[i] > x->x[i])
-      status = false;
-  }
-  return status;
-};
-
+      n[i] = 0;
+    if (n[i] >= N[i])
+      n[i] = N[i]-1;
+  }
+
+  // calculate distance of cell to vector
+  double distanceSquared = 0.;
+  bool outside = true;  // flag whether x is found in- or outside of LinkedCell's domain/closest cell
+  Vector corner; // current corner of closest cell
+  Vector tester; // Vector pointing from corner to center of closest cell
+  Vector Distance;  // Vector from corner of closest cell to x
+
+  Vector center;  // center of the closest cell
+  for (int i=0;i<NDIM;i++)
+    center.x[i] = min.x[i]+((double)n[i]+.5)*RADIUS;
+
+  int c[NDIM];
+  for (c[0]=0;c[0]<=1;c[0]++)
+    for (c[1]=0; c[1]<=1;c[1]++)
+      for (c[2]=0; c[2]<=1;c[2]++) {
+        // set up corner
+        for (int i=0;i<NDIM;i++)
+          corner.x[i] = min.x[i]+RADIUS*((double)n[i]+c[i]);
+        // set up distance vector
+        Distance.CopyVector(x);
+        Distance.SubtractVector(&corner);
+        const double dist = Distance.NormSquared();
+        // check whether distance is smaller
+        if (dist< distanceSquared)
+          distanceSquared = dist;
+        // check whether distance vector goes inside or outside
+        tester.CopyVector(&center);
+        tester.SubtractVector(&corner);
+        if (tester.ScalarProduct(&Distance) < 0)
+          outside = false;
+      }
+  return (outside ? distanceSquared : 0.);
+};
+
+/** Returns a list of all TesselPoint with distance less than \a radius to \a *Center.
+ * \param radius radius of sphere
+ * \param *center center of sphere
+ * \return list of all points inside sphere
+ */
+LinkedCell::LinkedNodes* LinkedCell::GetPointsInsideSphere(const double radius, const Vector * const center) const
+{
+  const double radiusSquared = radius*radius;
+  TesselPoint *Walker = NULL;
+  LinkedNodes *TesselList = new LinkedNodes;
+  LinkedNodes *NeighbourList = NULL;
+
+  // set index of LC to center of sphere
+  const double dist = SetClosestIndexToOutsideVector(center);
+  if (dist > 2.*radius) {
+    DoeLog(1) && (eLog()<< Verbose(1) << "Vector " << *center << " is too far away from any atom in LinkedCell's bounding box." << endl);
+    return TesselList;
+  } else
+    DoLog(1) && (Log() << Verbose(1) << "Distance of closest cell to center of sphere with radius " << radius << " is " << dist << "." << endl);
+
+  // gather all neighbours first, then look who fulfills distance criteria
+  NeighbourList = GetallNeighbours(2.*radius-dist);
+  //Log() << Verbose(1) << "I found " << NeighbourList->size() << " neighbours to check." << endl;
+  if (NeighbourList != NULL) {
+    for (LinkedNodes::const_iterator Runner = NeighbourList->begin(); Runner != NeighbourList->end(); Runner++) {
+      Walker = *Runner;
+      //Log() << Verbose(1) << "Current neighbour is at " << *Walker->node << "." << endl;
+      if ((center->DistanceSquared(Walker->node) - radiusSquared) < MYEPSILON) {
+        TesselList->push_back(Walker);
+      }
+    }
+    delete(NeighbourList);
+  } else
+    DoeLog(2) && (eLog()<< Verbose(2) << "Around vector " << *center << " there are no atoms." << endl);
+  return TesselList;
+};