source: src/linkedcell.cpp@ d103d3

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2011 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include "atom.hpp"
#include "linkedcell.hpp"
#include "CodePatterns/Verbose.hpp"
#include "CodePatterns/Range.hpp"
#include "CodePatterns/Log.hpp"
#include "molecule.hpp"
#include "IPointCloud.hpp"
#include "Tesselation/tesselation.hpp"
#include "LinearAlgebra/Vector.hpp"

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

/** Constructor for class LinkedCell.
 */
LinkedCell::LinkedCell() :
  LC(NULL),
  RADIUS(0.),
  index(-1)
{
  for(int i=0;i<NDIM;i++)
    N[i] = 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(IPointCloud & set, const double radius) :
  LC(NULL),
  RADIUS(radius),
  index(-1)
{
  TesselPoint *Walker = NULL;

  for(int i=0;i<NDIM;i++)
    N[i] = 0;
  max.Zero();
  min.Zero();
  LOG(1, "Begin of LinkedCell");
  if (set.IsEmpty()) {
    ELOG(1, "set is NULL or contains no linked cell nodes!");
    return;
  }
  // 1. find max and min per axis of atoms
  set.GoToFirst();
  Walker = set.GetPoint();
  for (int i=0;i<NDIM;i++) {
    max[i] = Walker->at(i);
    min[i] = Walker->at(i);
  }
  set.GoToFirst();
  while (!set.IsEnd()) {
    Walker = set.GetPoint();
    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);
    }
    set.GoToNext();
  }
  LOG(2, "Bounding box is " << min << " and " << max << ".");

  // 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);
  }
  LOG(2, "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << ".");

  // 3. allocate the lists
  LOG(2, "INFO: Allocating cells ... ");
  if (LC != NULL) {
    ELOG(1, "Linked Cell list is already allocated, I do nothing.");
    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 TesselPointSTLList[N[0]*N[1]*N[2]];
  for (index=0;index<N[0]*N[1]*N[2];index++) {
    LC [index].clear();
  }
  LOG(0, "done.");

  // 4. put each atom into its respective cell
  LOG(2, "INFO: Filling cells ... ");
  set.GoToFirst();
  while (!set.IsEnd()) {
    Walker = set.GetPoint();
    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(2, *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << ".");
    set.GoToNext();
  }
  LOG(0, "done.");
  LOG(1, "End of LinkedCell");
};


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

/** 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)
//    ELOG(1, "indices are out of bounds!");
  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 TesselPointSTLList* 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 TesselPointSTLList* 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;
  }
};

/** 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 & x) const
{
  for (int i=0;i<NDIM;i++)
    n[i] = (int)floor((x.at(i) - min[i])/RADIUS);

  return CheckBounds();
};

/** 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] = static_cast<int>(floor((Walker->at(i) - min[i])/RADIUS));
  }
  index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
  if (CheckBounds()) {
    for (TesselPointSTLList::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
      status = status || ((*Runner) == Walker);
    return status;
  } else {
    ELOG(1, "Node at " << *Walker << " is out of bounds.");
    return false;
  }
};

/** Calculates the interval bounds of the linked cell grid.
 * \param lower lower bounds
 * \param upper upper bounds
 * \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]-step;
    if (lower[i] < 0)
      lower[i] = 0;
    if (lower[i] >= N[i])
      lower[i] = N[i]-1;
    upper[i] = n[i]+step;
    if (upper[i] >= N[i])
      upper[i] = N[i]-1;
    if (upper[i] < 0)
      upper[i] = 0;
    //LOG(0, "axis " << i << " has bounds [" << lower[i] << "," << upper[i] << "]");
  }
};

/** 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
 */
TesselPointSTLList* LinkedCell::GetallNeighbours(const double distance) const
{
  int Nlower[NDIM], Nupper[NDIM];
  TesselPoint *Walker = NULL;
  TesselPointSTLList *TesselList = new TesselPointSTLList;

  // 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);

  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 TesselPointSTLList *List = GetCurrentCell();
        //LOG(1, "Current cell is " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << ".");
        if (List != NULL) {
          for (TesselPointSTLList::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->at(i) - min[i])/RADIUS);
    if (n[i] < 0)
      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[i] = min[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[i] = min[i]+RADIUS*((double)n[i]+c[i]);
        // set up distance vector
        Distance = (*x) - 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 = center -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
 */
TesselPointSTLList* LinkedCell::GetPointsInsideSphere(const double radius, const Vector * const center) const
{
  const double radiusSquared = radius*radius;
  TesselPoint *Walker = NULL;
  TesselPointSTLList *TesselList = new TesselPointSTLList;
  TesselPointSTLList *NeighbourList = NULL;

  // set index of LC to center of sphere
  const double dist = SetClosestIndexToOutsideVector(center);
  if (dist > 2.*radius) {
    ELOG(1, "Vector " << *center << " is too far away from any atom in LinkedCell's bounding box.");
    return TesselList;
  } else
    LOG(1, "Distance of closest cell to center of sphere with radius " << radius << " is " << dist << ".");

  // gather all neighbours first, then look who fulfills distance criteria
  NeighbourList = GetallNeighbours(2.*radius-dist);
  //LOG(1, "I found " << NeighbourList->size() << " neighbours to check.");
  if (NeighbourList != NULL) {
    for (TesselPointSTLList::const_iterator Runner = NeighbourList->begin(); Runner != NeighbourList->end(); Runner++) {
      Walker = *Runner;
      //LOG(1, "Current neighbour is at " << *Walker->node << ".");
      if ((Walker->DistanceSquared(*center) - radiusSquared) < MYEPSILON) {
        TesselList->push_back(Walker);
      }
    }
    delete(NeighbourList);
  } else
    ELOG(2, "Around vector " << *center << " there are no atoms.");
  return TesselList;
};