source: src/LinkedCell/LinkedCell_Model.cpp@ 5344e4

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

/*
 * LinkedCell_Model.cpp
 *
 *  Created on: Nov 15, 2011
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "LinkedCell_Model.hpp"

#include <algorithm>
#include <boost/multi_array.hpp>
#include <limits>

#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinkedCell/IPointCloud.hpp"
#include "LinkedCell/LinkedCell.hpp"
#include "Atom/TesselPoint.hpp"

#include "LinkedCell_Model_inline.hpp"

// initialize static entities
LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::NearestNeighbors;

/** Constructor of LinkedCell_Model.
 *
 * @param radius desired maximum neighborhood distance
 * @param _domain Box instance with domain size and boundary conditions
 */
LinkedCell::LinkedCell_Model::LinkedCell_Model(const double radius, const Box &_domain) :
    internal_Sizes(NULL),
    domain(_domain)
{
  // set default argument
  NearestNeighbors[0] = NearestNeighbors[1] = NearestNeighbors[2] = 1;

  setPartition(radius);

  // allocate linked cell structure
  AllocateCells();
}

/** Constructor of LinkedCell_Model.
 *
 * @oaram set set of points to place into the linked cell structure
 * @param radius desired maximum neighborhood distance
 * @param _domain Box instance with domain size and boundary conditions
 */
LinkedCell::LinkedCell_Model::LinkedCell_Model(IPointCloud &set, const double radius, const Box &_domain) :
    internal_Sizes(NULL),
    domain(_domain)
{
  Info info(__func__);

  // get the partition of the domain
  setPartition(radius);

  // allocate linked cell structure
  AllocateCells();

  insertPointCloud(set);
}

/** Destructor of class LinkedCell_Model.
 *
 */
LinkedCell::LinkedCell_Model::~LinkedCell_Model()
{
  Reset();
}


/** Allocates as much cells per axis as required by
 * LinkedCell_Model::BoxPartition.
 *
 */
void LinkedCell::LinkedCell_Model::AllocateCells()
{
  // resize array
  tripleIndex index;
  for (int i=0;i<NDIM;i++)
    index[i] = static_cast<LinkedCellArray::index>(Dimensions.at(i,i));
  N.resize(index);
  ASSERT(getSize(0)*getSize(1)*getSize(2) < MAX_LINKEDCELLNODES,
      "LinkedCell_Model::AllocateCells() - Number linked of linked cell nodes exceded hard-coded limit, use greater edge length!");

  // allocate LinkedCell instances
  for(index[0] = 0; index[0] != static_cast<LinkedCellArray::index>(Dimensions.at(0,0)); ++index[0]) {
    for(index[1] = 0; index[1] != static_cast<LinkedCellArray::index>(Dimensions.at(1,1)); ++index[1]) {
      for(index[2] = 0; index[2] != static_cast<LinkedCellArray::index>(Dimensions.at(2,2)); ++index[2]) {
        LOG(5, "INFO: Creating cell at " << index[0] << " " << index[1] << " " << index[2] << ".");
        N(index) = new LinkedCell(index);
      }
    }
  }
}

/** Frees all Linked Cell instances and sets array dimensions to (0,0,0).
 *
 */
void LinkedCell::LinkedCell_Model::Reset()
{
  // free all LinkedCell instances
  for(iterator3 iter3 = N.begin(); iter3 != N.end(); ++iter3) {
    for(iterator2 iter2 = (*iter3).begin(); iter2 != (*iter3).end(); ++iter2) {
      for(iterator1 iter1 = (*iter2).begin(); iter1 != (*iter2).end(); ++iter1) {
        delete *iter1;
      }
    }
  }
  // set dimensions to zero
  N.resize(boost::extents[0][0][0]);
}

/** Inserts all points contained in \a set.
 *
 * @param set set with points to insert into linked cell structure
 */
void LinkedCell::LinkedCell_Model::insertPointCloud(IPointCloud &set)
{
  if (set.IsEmpty()) {
    ELOG(1, "set is NULL or contains no linked cell nodes!");
    return;
  }

  // put each atom into its respective cell
  set.GoToFirst();
  while (!set.IsEnd()) {
    TesselPoint *Walker = set.GetPoint();
    addNode(Walker);
    set.GoToNext();
  }
}

/** Calculates the required edge length for the given desired distance.
 *
 * We need to make some matrix transformations in order to obtain the required
 * edge lengths per axis. Goal is guarantee that whatever the shape of the
 * domain that always all points at least up to \a distance away are contained
 * in the nearest neighboring cells.
 *
 * @param distance distance of this linked cell array
 */
void LinkedCell::LinkedCell_Model::setPartition(double distance)
{
  // generate box matrix of desired edge length
  RealSpaceMatrix neighborhood;
  neighborhood.setIdentity();
  neighborhood *= distance;

  // obtain refs to both domain matrix transformations
  //const RealSpaceMatrix &M = domain.getM();
  const RealSpaceMatrix &Minv = domain.getMinv();

  RealSpaceMatrix output = Minv * neighborhood;

  std::cout << Minv << " * " << neighborhood << " = " << output << std::endl;

  Dimensions = output.invert();
  Partition = Minv*Dimensions; //

  std::cout << "Dimensions are then " << Dimensions << std::endl;
  std::cout << "Partition matrix is then " << Partition << std::endl;
}

/** Returns the number of required neighbor-shells to get all neighboring points
 * in the given \a distance.
 *
 * @param distance radius of neighborhood sphere
 * @return number of LinkedCell's per dimension to get all neighbors
 */
const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getStep(const double distance) const
{
  tripleIndex index;
  index[0] = index[1] = index[2] = 0;

  if (fabs(distance) < std::numeric_limits<double>::min())
    return index;
  // generate box matrix of desired edge length
  RealSpaceMatrix neighborhood;
  neighborhood.setIdentity();
  neighborhood *= distance;

  const RealSpaceMatrix output = Partition * neighborhood;

  //std::cout << "GetSteps: " << Partition << " * " << neighborhood << " = " << output << std::endl;

  const RealSpaceMatrix steps = output;
  for (size_t i =0; i<NDIM; ++i)
    index[i] = ceil(steps.at(i,i));
  LOG(2, "INFO: number of shells are ("+toString(index[0])+","+toString(index[1])+","+toString(index[2])+").");

  return index;
}

/** Calculates the index of the cell \a position would belong to.
 *
 * @param position position whose associated cell to calculate
 * @return index of the cell
 */
const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getIndexToVector(const Vector &position) const
{
  tripleIndex index;
  Vector x(Partition*position);
  LOG(2, "INFO: Transformed position is " << x << ".");
  for (int i=0;i<NDIM;i++) {
    index[i] = static_cast<LinkedCellArray::index>(floor(x[i]));
  }
  return index;
}

/** Adds a node to the linked cell structure
 *
 * @param Walker node to add
 */
void LinkedCell::LinkedCell_Model::addNode(TesselPoint *&Walker)
{
  tripleIndex index = getIndexToVector(Walker->getPosition());
  LOG(2, "INFO: " << *Walker << " goes into cell " << index[0] << ", " << index[1] << ", " << index[2] << ".");
  LOG(2, "INFO: Cell's indices are "
      << N(index)->getIndex(0) << " "
      << N(index)->getIndex(1) << " "
      << N(index)->getIndex(2) << ".");
  N(index)->addPoint(Walker);
  std::pair<MapPointToCell::iterator, bool> inserter = CellLookup.insert( std::make_pair(Walker, N(index)) );
  ASSERT( inserter.second,
      "LinkedCell_Model::addNode() - Walker "
      +toString(*Walker)+" is already present in cell "
      +toString((inserter.first)->second->getIndex(0))+" "
      +toString((inserter.first)->second->getIndex(1))+" "
      +toString((inserter.first)->second->getIndex(2))+".");
}

/** Removes a node to the linked cell structure
 *
 * We do nothing of Walker is not found
 *
 * @param Walker node to remove
 */
void LinkedCell::LinkedCell_Model::deleteNode(const TesselPoint *Walker)
{
  MapPointToCell::iterator iter = CellLookup.begin();
  for (; iter != CellLookup.end(); ++iter)
    if (iter->first == Walker)
      break;
  ASSERT(iter != CellLookup.end(),
      "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
  if (iter != CellLookup.end()) {
    CellLookup.erase(iter);
    iter->second->deletePoint(Walker);
  }
}

/** Move Walker from current cell to another on position change.
 *
 * @param Walker node who has moved.
 */
void LinkedCell::LinkedCell_Model::moveNode(const TesselPoint *Walker)
{
  ASSERT(0, "LinkedCell_Model::moveNode() - not implemented yet.");
}

/** Checks whether cell indicated by \a relative relative to LinkedCell_Model::internal_index
 * is out of bounds.
 *
 * \note We do not check for boundary conditions of LinkedCell_Model::domain,
 * we only look at the array sizes
 *
 * @param relative index relative to LinkedCell_Model::internal_index.
 * @return true - relative index is still inside bounds, false - outside
 */
bool LinkedCell::LinkedCell_Model::checkArrayBounds(const tripleIndex &index) const
{
  bool status = true;
  for (size_t i=0;i<NDIM;++i) {
    status = status && (
        (index[i] >= 0) &&
        (index[i] < getSize(i))
        );
  }
  return status;
}

/** Corrects \a index according to boundary conditions of LinkedCell_Model::domain .
 *
 * @param index index to correct according to boundary conditions
 */
void LinkedCell::LinkedCell_Model::applyBoundaryConditions(tripleIndex &index) const
{
  for (size_t i=0;i<NDIM;++i) {
    switch (domain.getConditions()[i]) {
    case Box::Wrap:
      if ((index[i] < 0) || (index[i] >= getSize(i)))
        index[i] = (index[i] % getSize(i));
      break;
    case Box::Bounce:
      break;
    case Box::Ignore:
      break;
    default:
      ASSERT(0, "LinkedCell_Model::checkBounds() - unknown boundary conditions.");
      break;
    }
  }
}

/** Calculates the interval bounds of the linked cell grid.
 * \param index index to give relative bounds to
 * \param step how deep to check the neighbouring cells (i.e. number of layers to check)
 * \return pair of tripleIndex indicating lower and upper bounds
 */
const LinkedCell::LinkedCell_Model::LinkedCellNeighborhoodBounds LinkedCell::LinkedCell_Model::getNeighborhoodBounds(
    const tripleIndex &index,
    const tripleIndex &step
    ) const
{
  LinkedCellNeighborhoodBounds neighbors;

  // calculate bounds
  for (size_t i=0;i<NDIM;++i) {
    switch (domain.getConditions()[i]) {
      case Box::Wrap:
        if (index[i] - step[i] >= 0) {
          neighbors.first[i] = index[i] - step[i];
        } else {
          neighbors.first[i] = getSize(i) - (index[i] - step[i]);
        }
        if (index[i] + step[i] < getSize(i)) {
          neighbors.second[i] = index[i] + step[i];
        } else {
          neighbors.second[i] = index[i] + step[i] - getSize(i);
        }
        break;
      case Box::Bounce:
        if (index[i] - step[i] >= 0) {
          neighbors.first[i] = index[i] - step[i];
        } else {
          neighbors.first[i] = - (index[i] - step[i]);
        }
        if (index[i] + step[i] < getSize(i)) {
          neighbors.second[i] = index[i] + step[i];
        } else {
          const size_t shift = index[i] + step[i];
          const size_t mod = shift / getSize(i);
          if ((mod / 2)*2 == mod) // even -> come in from lower bound
            neighbors.second[i] = (shift % getSize(i));
          else // odd -> come in from upper bound
            neighbors.second[i] = getSize(i) - (shift % getSize(i));
        }
        break;
      case Box::Ignore:
        ASSERT(index[i] - step[i] >= 0,
            "LinkedCell_Model::getNeighborhoodBounds() - lower out of bounds on 'Ignore' boundary conditions");
        neighbors.first[i] = index[i] - step[i];
        ASSERT (index[i] + step[i] < getSize(i),
            "LinkedCell_Model::getNeighborhoodBounds() - upper out of bounds on 'Ignore' boundary conditions");
        neighbors.second[i] = index[i] + step[i];
        ASSERT(neighbors.second[i] < neighbors.first[i],
            "LinkedCell_Model::getNeighborhoodBounds() - ["
            +toString(neighbors.first[i])+","+toString(neighbors.second[i])
            +"] does not make sense as array bounds.");
        break;
      default:
        ASSERT(0, "LinkedCell_Model::getNeighborhoodBounds() - unknown boundary conditions.");
        break;
    }
  }

  return neighbors;
}