source: molecuilder/src/linkedcell.cpp@ 6dd8d3

/** \file linkedcell.cpp
 *
 * Function implementations for the class LinkedCell.
 *
 */


#include "atom.hpp"
#include "helpers.hpp"
#include "linkedcell.hpp"
#include "log.hpp"
#include "molecule.hpp"
#include "tesselation.hpp"
#include "vector.hpp"

// ========================================================= class LinkedCell ===========================================


/** 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 *set LCNodeSet class with all LCNode's
 * \param RADIUS edge length of cells
 */
LinkedCell::LinkedCell(const PointCloud * const set, const double radius)
{
  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 == NULL) || (set->IsEmpty())) {
    DoeLog(1) && (eLog()<< Verbose(1) << "set is NULL or contains no linked cell nodes!" << endl);
    return;
  }
  // 1. find max and min per axis of atoms
  set->GoToFirst();
  Walker = set->GetPoint();
  for (int i=0;i<NDIM;i++) {
    max.x[i] = Walker->node->x[i];
    min.x[i] = Walker->node->x[i];
  }
  set->GoToFirst();
  while (!set->IsEnd()) {
    Walker = set->GetPoint();
    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];
    }
    set->GoToNext();
  }
  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) {
    DoeLog(1) && (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 ... ";
  set->GoToFirst();
  while (!set->IsEnd()) {
    Walker = set->GetPoint();
    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;
    set->GoToNext();
  }
  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()) {
    DoeLog(1) && (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) {
    DoeLog(1) && (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;
};

/** 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() const
{
  bool status = true;
  for(int i=0;i<NDIM;i++)
    status = status && ((n[i] >=0) && (n[i] < N[i]));
  if (!status)
  DoeLog(1) && (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;
};


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

/** 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 LinkedCell::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(const TesselPoint * const Walker) const
{
  bool status = false;
  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];
  if (CheckBounds()) {
    for (LinkedNodes::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
      status = status || ((*Runner) == Walker);
    return status;
  } else {
    DoeLog(1) && (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;
  }
};

/** 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];
  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] << "] ";
  }
  //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;
};


/** 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
{
  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 {
    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);
            }
          }
        }
      }
  return TesselList;
};