source: molecuilder/src/linkedcell.cpp@ a87d5e6

#include "linkedcell.hpp"
#include "molecules.hpp"

/** Constructor for class LinkedCell.
 */
LinkedCell::LinkedCell()
{
        LC = NULL;
        for(int i=0;i<NDIM;i++)
                N[i] = 0;
        index = -1;
        RADIUS = 0.;
        max.Zero();
        min.Zero();
};

/** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
 * \param *mol molecule structure with all Atom's
 * \param RADIUS edge length of cells
 */
LinkedCell::LinkedCell(molecule *mol, double radius)
{
        atom *Walker = NULL;

        RADIUS = radius;
        LC = NULL;
        for(int i=0;i<NDIM;i++)
                N[i] = 0;
        index = -1;
        max.Zero();
        min.Zero();
        cout << Verbose(1) << "Begin of LinkedCell" << endl;
        if (mol->start->next == mol->end) {
                cerr << "ERROR: molecule contains no atoms!" << endl;
                return;
        }
        // 1. find max and min per axis of atoms
        Walker = mol->start->next;
        for (int i=0;i<NDIM;i++) {
                max.x[i] = Walker->x.x[i];
                min.x[i] = Walker->x.x[i];
        }
        while (Walker != mol->end) {
                for (int i=0;i<NDIM;i++) {
                        if (max.x[i] < Walker->x.x[i])
                                max.x[i] = Walker->x.x[i];
                        if (min.x[i] > Walker->x.x[i])
                                min.x[i] = Walker->x.x[i];
                }
                Walker = Walker->next;
        }
        cout << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;

        // 2. find then umber of cells per axis
        for (int i=0;i<NDIM;i++) {
                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;

        // 3. allocate the lists
        cout << Verbose(2) << "Allocating cells ... ";
        if (LC != NULL) {
                cout << Verbose(1) << "ERROR: Linked Cell list is already allocated, I do nothing." << endl;
                return;
        }
        LC = new LinkedAtoms[N[0]*N[1]*N[2]];
        for (index=0;index<N[0]*N[1]*N[2];index++) {
                LC [index].clear();
        }
        cout << "done."  << endl;

        // 4. put each atom into its respective cell
  cout << Verbose(2) << "Filling cells ... ";
        Walker = mol->start;
        while (Walker->next != mol->end) {
                Walker = Walker->next;
                for (int i=0;i<NDIM;i++) {
                        n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS);
                }
                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;
        }
  cout << "done."  << endl;
        cout << Verbose(1) << "End of LinkedCell" << endl;
};

/** Destructor for class LinkedCell.
 */
LinkedCell::~LinkedCell()
{
        if (LC != NULL)
        for (index=0;index<N[0]*N[1]*N[2];index++)
                LC[index].clear();
        delete[](LC);
        for(int i=0;i<NDIM;i++)
                N[i] = 0;
        index = -1;
        max.Zero();
        min.Zero();
};

/** Checks whether LinkedCell::n[] is each within [0,N[]].
 * \return if all in intervals - true, else -false
 */
bool LinkedCell::CheckBounds()
{
        bool status = true;
        for(int i=0;i<NDIM;i++)
                status = status && ((n[i] >=0) && (n[i] < N[i]));
        if (!status)
        cerr << "ERROR: indices are out of bounds!" << endl;
        return status;
};


/** Returns a pointer to the current cell.
 * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
 */
LinkedAtoms* LinkedCell::GetCurrentCell()
{
        if (CheckBounds()) {
                index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
                return (&(LC[index]));
        } else {
                return NULL;
        }
};

/** Calculates the index for a given atom *Walker.
 * \param Walker atom to set index to
 * \return if the atom is also found in this cell - true, else - false
 */
bool LinkedCell::SetIndexToAtom(const atom &Walker)
{
        bool status = false;
        for (int i=0;i<NDIM;i++) {
                n[i] = (int)floor((Walker.x.x[i] - min.x[i])/RADIUS);
        }

        index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
        if (CheckBounds()) {
                for (LinkedAtoms::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
                        status = status || ((*Runner) == &Walker);
                return status;
        } else {
                cerr << Verbose(1) << "WARN: Atom "<< Walker << " at " << Walker.x << " 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])
{
  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;
    //cout << " [" << 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];

    //cout << "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 *x)
{
        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;
};