source: src/linkedcell.hpp@ 1a6bda

/*
 * linkedcell.hpp
 *
 *  If the linked cell should be usable, the class has to inherit LCNodeSet and the nodes (containing the Vectors) have to inherit LCNode. This works well
 *  for molecule and atom classes.
 *
 *  Created on: Aug 3, 2009
 *      Author: heber
 */

#ifndef LINKEDCELL_HPP_
#define LINKEDCELL_HPP_

using namespace std;

/*********************************************** includes ***********************************/

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <list>
#include <vector>
#include <map>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Helpers/defs.hpp"
#include "World.hpp"
#include "LinearAlgebra/Vector.hpp"

/****************************************** forward declarations *****************************/

class PointCloud;
class TesselPoint;

/********************************************** definitions *********************************/

//< Upper bound for number of cell nodes used in sensibility check on allocation.
enum { MAX_LINKEDCELLNODES = 1000000 };

/********************************************** declarations *******************************/

/** Linked Cell class for containing Vectors in real space efficiently.
 */
class LinkedCell {
private:

public:
  class LinkedNodes : public std::list<TesselPoint *> {
  public:
    LinkedNodes();
    ~LinkedNodes();

    TesselPoint * getValue (const_iterator &rhs) const;
    TesselPoint * getValue (iterator &rhs) const;
  };
  //typedef list<TesselPoint *> LinkedNodes;


    Vector max;       // upper boundary
    Vector min;       // lower boundary
    LinkedNodes *LC;  // linked cell list
    double RADIUS;    // cell edge length
    int N[NDIM];      // number of cells per axis
    mutable int n[NDIM];      // temporary variable for current cell per axis
    mutable int index;        // temporary index variable , access by index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];

    LinkedCell();
    LinkedCell(const PointCloud &set, const double radius);
    template <class T> LinkedCell(const T &set, const double radius) :
      LC(NULL),
      RADIUS(radius),
      index(-1)
    {
      class TesselPoint *Walker = NULL;
      for(int i=0;i<NDIM;i++)
        N[i] = 0;
      max.Zero();
      min.Zero();
      DoLog(1) && (Log() << Verbose(1) << "Begin of LinkedCell" << endl);
      if (set.begin() == set.end()) {
        DoeLog(1) && (eLog()<< Verbose(1) << "set contains no linked cell nodes!" << endl);
        return;
      }
      // 1. find max and min per axis of atoms
      typename T::const_iterator Runner = set.begin();
      for (int i=0;i<NDIM;i++) {
        max[i] = set.getValue(Runner)->at(i);
        min[i] = set.getValue(Runner)->at(i);
      }
      for (typename T::const_iterator Runner = set.begin(); Runner != set.end(); Runner++) {
        Walker = set.getValue(Runner);
        for (int i=0;i<NDIM;i++) {
          if (max[i] < Walker->at(i))
            max[i] = Walker->at(i);
          if (min[i] > Walker->at(i))
            min[i] = Walker->at(i);
        }
      }
      DoLog(2) && (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] = static_cast<int>(floor((max[i] - min[i])/RADIUS)+1);
      }
      DoLog(2) && (Log() << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl);

      // 3. allocate the lists
      DoLog(2) && (Log() << Verbose(2) << "Allocating cells ... ");
      if (LC != NULL) {
        DoeLog(1) && (eLog()<< Verbose(1) << "Linked Cell list is already allocated, I do nothing." << endl);
        return;
      }
      ASSERT(N[0]*N[1]*N[2] < MAX_LINKEDCELLNODES, "Number linked of linked cell nodes exceded hard-coded limit, use greater edge length!");
      LC = new LinkedNodes[N[0]*N[1]*N[2]];
      for (index=0;index<N[0]*N[1]*N[2];index++) {
        LC [index].clear();
      }
      DoLog(0) && (Log() << Verbose(0) << "done."  << endl);

      // 4. put each atom into its respective cell
      DoLog(2) && (Log() << Verbose(2) << "Filling cells ... ");
      for (typename T::const_iterator Runner = set.begin(); Runner != set.end(); Runner++) {
        Walker = set.getValue(Runner);
        for (int i=0;i<NDIM;i++) {
          n[i] = static_cast<int>(floor((Walker->at(i) - min[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;
      }
      DoLog(0) && (Log() << Verbose(0) << "done."  << endl);
      DoLog(1) && (Log() << Verbose(1) << "End of LinkedCell" << endl);
    };


    ~LinkedCell();
    const LinkedCell::LinkedNodes* GetCurrentCell()const ;
    const LinkedCell::LinkedNodes* GetRelativeToCurrentCell(const int relative[NDIM])const ;
    bool SetIndexToNode(const TesselPoint * const Walker)const ;
    bool SetIndexToVector(const Vector &x)const ;
    double SetClosestIndexToOutsideVector(const Vector * const x) const;
    bool CheckBounds()const ;
    bool CheckBounds(const int relative[NDIM])const ;
    void GetNeighbourBounds(int lower[NDIM], int upper[NDIM], int step = 1)const ;

    LinkedCell::LinkedNodes* GetallNeighbours(const double distance = 0) const;
    LinkedCell::LinkedNodes* GetPointsInsideSphere(const double radius, const Vector * const center) const;
    // not implemented yet
    bool AddNode(Vector *Walker);
    bool DeleteNode(Vector *Walker);
    bool MoveNode(Vector *Walker);
};

#endif /*LINKEDCELL_HPP_*/