Changeset ef016c for molecuilder/src

Timestamp:

Apr 9, 2010, 1:39:03 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Children:

e760c1

Parents:

bd8561

Message:

changes to LinkedCell::GetNeighbourBounds(), LinkedCell::GetallNeighbours() and LinkedCell::GetPointsInsideSphere()

• LinkedCell::GetNeighbourBounds() has a step up to which step-nearest neighbours cells are included, heeding boundaries
• LinkedCell::GetallNeighbours() uses GetNeighbourBounds()
• LinkedCell::GetPointsInsideSphere() uses GetallNeighbours() and then removes all points outside sphere from this list.

Location:

molecuilder/src

Files:

• linkedcell.cpp (modified) (4 diffs)
• linkedcell.hpp (modified) (1 diff)

molecuilder/src/linkedcell.cpp

@@ -244 +244 @@
 };
 
+/** 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();
+};
+
 /** Calculates the index for a given LCNode *Walker.
  * \param *Walker LCNode to set index tos
@@ -269 +281 @@
  * \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] << "] ";
+void LinkedCell::GetNeighbourBounds(int lower[NDIM], int upper[NDIM], int step) const
+{
+  for (int i=0;i<NDIM;i++) {
+    lower[i] = 0;
+    for (int s=step; s>0;--s)
+      if ((n[i]-s) >= 0) {
+        lower[i] = n[i]-s;
+        break;
+      }
+    upper[i] = 0;
+    for (int s=step; s>0;--s)
+      if ((n[i]+s) < N[i]) {
+        upper[i] = n[i]+s;
+        break;
+      }
+    //Log() << Verbose(0) << " [" << lower[i] << "," << upper[i] << "] ";
     // check for this axis whether the point is outside of our grid
     if (n[i] < 0)
@@ -285 +307 @@
 };
 
-/** 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;
-};
-
 /** 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 = 0) const
-{
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+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
-  //Log() << Verbose(1) << "LC Intervals:";
-  const int step = (distance == 0) ? 1 : (int)floor(distance/RADIUS)+1;
-  for (int i=0;i<NDIM;i++) {
-    Nlower[i] = ((N[i]-step) >= 0) ? N[i]-step : 0;
-    Nupper[i] = ((N[i]+step) < N[i]) ? N[i]+step : N[i]-step;
-    //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
-  }
+  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]++)
@@ -344 +345 @@
 LinkedCell::LinkedNodes* LinkedCell::GetPointsInsideSphere(const double radius, const Vector * const center) const
 {
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
   const double radiusSquared = radius*radius;
   TesselPoint *Walker = NULL;
   LinkedNodes *TesselList = new LinkedNodes;
-
-  if (SetIndexToVector(center)) {
-    for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = n[i];
-    //Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << index << "." << endl;
-  } else {
+  LinkedNodes *NeighbourList = NULL;
+
+  // set index of LC to center of sphere
+  if (!SetIndexToVector(center)) {
     DoeLog(1) && (eLog()<< Verbose(1) << "Vector " << *center << " is outside of LinkedCell's bounding box." << endl);
     return TesselList;
   }
-  // then go through the current and all neighbouring cells and check the contained points for possible candidates
-  //Log() << Verbose(1) << "LC Intervals:";
-  for (int i=0;i<NDIM;i++) {
-    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-    Nupper[i] = ((N[i]+1) < N[i]) ? N[i]+1 : N[i]-1;
-    //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
-  }
-  //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;
-            if ((center->DistanceSquared(Walker->node) - radiusSquared) > MYEPSILON) {
-              TesselList->push_back(Walker);
-            }
-          }
-        }
+
+  // gather all neighbours first, then look who fulfills distance criteria
+  NeighbourList = GetallNeighbours(radius);
+  if (NeighbourList != NULL) {
+    for (LinkedNodes::const_iterator Runner = NeighbourList->begin(); Runner != NeighbourList->end(); Runner++) {
+      Walker = *Runner;
+      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;
 };

molecuilder/src/linkedcell.hpp

@@ -63 +63 @@
     bool CheckBounds()const ;
     bool CheckBounds(const int relative[NDIM])const ;
-    void GetNeighbourBounds(int lower[NDIM], int upper[NDIM])const ;
+    void GetNeighbourBounds(int lower[NDIM], int upper[NDIM], int step = 1)const ;
 
-    LinkedCell::LinkedNodes* GetallNeighbours(const double distance) const;
+    LinkedCell::LinkedNodes* GetallNeighbours(const double distance = 0) const;
     LinkedCell::LinkedNodes* GetPointsInsideSphere(const double radius, const Vector * const center) const;