Index: src/LinkedCell/LinkedCell_Model.cpp
===================================================================
--- src/LinkedCell/LinkedCell_Model.cpp	(revision d8296152e374a74c483120eb411121989ba27a98)
+++ src/LinkedCell/LinkedCell_Model.cpp	(revision d8296152e374a74c483120eb411121989ba27a98)
@@ -0,0 +1,388 @@
+/*
+ * 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;
+}
Index: src/LinkedCell/LinkedCell_Model.hpp
===================================================================
--- src/LinkedCell/LinkedCell_Model.hpp	(revision d8296152e374a74c483120eb411121989ba27a98)
+++ src/LinkedCell/LinkedCell_Model.hpp	(revision d8296152e374a74c483120eb411121989ba27a98)
@@ -0,0 +1,114 @@
+/*
+ * LinkedCell_Model.hpp
+ *
+ *  Created on: Nov 15, 2011
+ *      Author: heber
+ */
+
+#ifndef LINKEDCELL_MODEL_HPP_
+#define LINKEDCELL_MODEL_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <boost/multi_array.hpp>
+#include <map>
+
+#include "LinearAlgebra/defs.hpp"
+#include "LinearAlgebra/RealSpaceMatrix.hpp"
+#include "LinkedCell/types.hpp"
+
+class Box;
+class IPointCloud;
+class LinkedCell_ModelTest;
+class TesselPoint;
+class Vector;
+
+namespace LinkedCell {
+
+  /** This is the model in the MVC ansatz for the LinkedCell structure.
+   *
+   * \sa linkedcell
+   *
+   * The model represents a certain linked cell structure with a specific mesh
+   * size (i.e. edge length). Note that all coordinates are internally always
+   * transformed to [0,1]^3, hence we require the domain (\ref Box) for the
+   * transformation matrices.
+   *
+   * This class stores internally list of particles in three-dimensional arrays.
+   *
+   */
+  class LinkedCell_Model
+  {
+    //!> grant unit test access to internal structure
+    friend class ::LinkedCell_ModelTest;
+  public:
+    LinkedCell_Model(const double radius, const Box &_domain);
+    LinkedCell_Model(IPointCloud &set, const double radius, const Box &_domain);
+    ~LinkedCell_Model();
+
+    typedef std::pair<tripleIndex,tripleIndex> LinkedCellNeighborhoodBounds;
+
+    const tripleIndex getIndexToVector(const Vector &position) const;
+    const LinkedCellNeighborhoodBounds getNeighborhoodBounds(const tripleIndex &index, const tripleIndex &step = NearestNeighbors) const;
+    const tripleIndex getStep(const double distance) const;
+    const LinkedCell& getCell(const tripleIndex &index) const;
+    bool checkArrayBounds(const tripleIndex &index) const;
+    void applyBoundaryConditions(tripleIndex &index) const;
+
+    void addNode(TesselPoint *&Walker);
+    void deleteNode(const TesselPoint *Walker);
+    void moveNode(const TesselPoint *Walker);
+
+    void setPartition(double distance);
+
+    //!> static indices to get nearest neighbors as default argument
+    static tripleIndex NearestNeighbors;
+
+  private:
+    void AllocateCells();
+    void Reset();
+    void insertPointCloud(IPointCloud &set);
+
+    LinkedCellArray::index getSize(const size_t dim) const;
+
+    typedef LinkedCellArray::size_type size_type;
+    typedef LinkedCellArray::iterator iterator3;
+    typedef boost::subarray_gen<LinkedCellArray, NDIM-1>::type::iterator iterator2;
+    typedef boost::subarray_gen<LinkedCellArray, NDIM-2>::type::iterator iterator1;
+
+    typedef std::map<TesselPoint *, LinkedCell *> MapPointToCell;
+
+    //!> internal shape of multi_array
+    size_type *internal_Sizes;
+
+    //!> internal index of current cell, this makes going through linked cell faster
+    tripleIndex internal_index;
+
+    //!> Lookup map to get the cell for a given TesselPoint
+    MapPointToCell CellLookup;
+
+    //!> edge length per axis
+    boost::array<double, 3> EdgeLength;
+
+    //!> Linked cell array
+    LinkedCellArray N;
+
+    //!> matrix to divide Box with
+    RealSpaceMatrix Dimensions;
+
+    //!> matrix to transform normal position to obtain partitioned position
+    RealSpaceMatrix Partition;
+
+    //!> Box with matrix transformations and boundary conditions
+    const Box &domain;
+  };
+
+}
+
+// inlined functions
+#include "LinkedCell_Model_inline.hpp"
+
+#endif /* LINKEDCELL_MODEL_HPP_ */
Index: src/LinkedCell/LinkedCell_Model_inline.hpp
===================================================================
--- src/LinkedCell/LinkedCell_Model_inline.hpp	(revision d8296152e374a74c483120eb411121989ba27a98)
+++ src/LinkedCell/LinkedCell_Model_inline.hpp	(revision d8296152e374a74c483120eb411121989ba27a98)
@@ -0,0 +1,51 @@
+/*
+ * LinkedCell_Model_inline.hpp
+ *
+ *  Created on: Nov 18, 2011
+ *      Author: heber
+ */
+
+#ifndef LINKEDCELL_MODEL_INLINE_HPP_
+#define LINKEDCELL_MODEL_INLINE_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// have all includes only needed here extra
+#ifdef HAVE_INLINE
+#include "CodePatterns/Assert.hpp"
+#endif
+
+/** Returns a reference to the cell indicated by LinkedCell_Model::internal_index.
+ *
+ * \return LinkedCell ref to current cell
+ */
+#ifdef HAVE_INLINE
+inline
+#endif
+const LinkedCell::LinkedCell& LinkedCell::LinkedCell_Model::getCell(const tripleIndex &index) const
+{
+  return *N(index);
+}
+
+
+/** Returns size of array for given \a dim.
+ *
+ * @param dim desired dimension
+ * @return size of array along dimension
+ */
+#ifdef HAVE_INLINE
+inline
+#endif
+LinkedCell::LinkedCellArray::index LinkedCell::LinkedCell_Model::getSize(const size_t dim) const
+{
+  ASSERT((dim >= 0) && (dim < NDIM),
+      "LinkedCell_Model::getSize() - dimension "
+      +toString(dim)+" is out of bounds.");
+  return N.shape()[dim];
+}
+
+
+#endif /* LINKEDCELL_MODEL_INLINE_HPP_ */
Index: src/LinkedCell/unittests/LinkedCell_ModelUnitTest.cpp
===================================================================
--- src/LinkedCell/unittests/LinkedCell_ModelUnitTest.cpp	(revision d8296152e374a74c483120eb411121989ba27a98)
+++ src/LinkedCell/unittests/LinkedCell_ModelUnitTest.cpp	(revision d8296152e374a74c483120eb411121989ba27a98)
@@ -0,0 +1,255 @@
+/*
+ * 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_ModelUnitTest.cpp
+ *
+ *  Created on: Nov 17, 2011
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <list>
+
+#include "Box.hpp"
+#include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
+#include "LinearAlgebra/RealSpaceMatrix.hpp"
+#include "LinkedCell/LinkedCell.hpp"
+#include "LinkedCell/LinkedCell_Model.hpp"
+#include "LinkedCell/PointCloudAdaptor.hpp"
+
+#include "LinkedCell_ModelUnitTest.hpp"
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+/********************************************** Test classes **************************************/
+
+#define DOMAINLENGTH 20.
+#define EDGELENGTH 2.
+#define NUMBERCELLS DOMAINLENGTH/EDGELENGTH
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( LinkedCell_ModelTest );
+
+
+void LinkedCell_ModelTest::setUp()
+{
+  // failing asserts should be thrown
+  ASSERT_DO(Assert::Throw);
+
+  setVerbosity(3);
+
+  // create diag(20.) matrix
+  RealSpaceMatrix BoxM;
+  BoxM.setIdentity();
+  BoxM *= DOMAINLENGTH;
+
+  // create Box with this matrix
+  domain = new Box(BoxM);
+
+  // create LinkedCell structure with this Box
+  LC = new LinkedCell::LinkedCell_Model(EDGELENGTH, *domain);
+
+  // create a list of nodes and add to LCImpl
+  std::vector< Vector > VectorList;
+  for (size_t i=0;i<((size_t)floor(NUMBERCELLS));++i)
+    VectorList.push_back(Vector((double)i*EDGELENGTH,(double)i*EDGELENGTH,(double)i*EDGELENGTH));
+  for (size_t i=0;i<VectorList.size();++i) {
+    TesselPoint * Walker = new TesselPoint();
+    Walker->setName(std::string("Walker")+toString(i));
+    Walker->setPosition(VectorList[i]);
+    NodeList.insert(Walker);
+  }
+}
+
+
+void LinkedCell_ModelTest::tearDown()
+{
+  delete LC;
+  delete domain;
+
+  // remove all nodes again
+  for (PointSet::iterator iter = NodeList.begin();
+      !NodeList.empty();
+      iter = NodeList.begin()) {
+    delete *iter;
+    NodeList.erase(iter);
+  }
+}
+
+/** UnitTest for correct construction
+ */
+void LinkedCell_ModelTest::AllocationTest()
+{
+  // check that first cell is allocated
+  LinkedCell::tripleIndex index;
+  index[0] = index[1] = index[2] = 0;
+  CPPUNIT_ASSERT(LC->N(index) != NULL);
+
+  // check that very last cell is allocated
+  index[0] = index[1] = index[2] = (size_t)floor(NUMBERCELLS)-1;
+  CPPUNIT_ASSERT(LC->N(index) != NULL);
+
+}
+
+/** UnitTest for getSize()
+ */
+void LinkedCell_ModelTest::getSizeTest()
+{
+  // check getSize()
+  for(size_t i=0; i<NDIM; ++i)
+    CPPUNIT_ASSERT_EQUAL((LinkedCell::LinkedCellArray::index)floor(NUMBERCELLS), LC->getSize(i));
+#ifndef NDEBUG
+  std::cout << "The following assertion is intended and is not a failure of the code." << std::endl;
+  CPPUNIT_ASSERT_THROW( LC->getSize(4), Assert::AssertionFailure);
+#endif
+}
+
+/** UnitTest for Reset()
+ */
+void LinkedCell_ModelTest::ResetTest()
+{
+  LC->Reset();
+
+  for(size_t i=0; i<NDIM; ++i)
+    CPPUNIT_ASSERT_EQUAL((LinkedCell::LinkedCellArray::index)0, LC->getSize(i));
+}
+
+
+/** UnitTest for insertPointCloud()
+ */
+void LinkedCell_ModelTest::insertPointCloudTest()
+{
+
+  // create the linked cell structure
+  PointCloudAdaptor< PointSet > cloud(&NodeList, std::string("NodeList"));
+  LC->insertPointCloud(cloud);
+
+  // test structure
+  CPPUNIT_ASSERT_EQUAL(((size_t)floor(NUMBERCELLS)), LC->CellLookup.size());
+  LinkedCell::tripleIndex index;
+  for (size_t i=0;i<((size_t)floor(NUMBERCELLS));++i) {
+    index[0] = index[1] = index[2] = i;
+    // assert that in the destined cell we have one Walker
+    CPPUNIT_ASSERT_EQUAL((size_t)1, LC->N(index)->size());
+  }
+  index[0] = 9;
+  index[1] = index[2] = 0;
+  // assert that in the destined cell we have one Walker
+  CPPUNIT_ASSERT_EQUAL((size_t)0, LC->N(index)->size());
+
+}
+
+/** UnitTest for insertPointCloud()
+ */
+void LinkedCell_ModelTest::setPartitionTest()
+{
+  RealSpaceMatrix Pmatrix = LC->Partition;
+  RealSpaceMatrix Dmatrix = LC->Dimensions;
+
+  LC->Reset();
+
+  LC->setPartition(2.*EDGELENGTH);
+
+  Pmatrix *= 0.5;
+  Dmatrix *= 0.5;
+
+  CPPUNIT_ASSERT_EQUAL(Pmatrix, LC->Partition);
+  CPPUNIT_ASSERT_EQUAL(Dmatrix, LC->Dimensions);
+}
+
+/** UnitTest for insertPointCloud()
+ */
+void LinkedCell_ModelTest::getStepTest()
+{
+  // zero distance
+  LinkedCell::tripleIndex index = LC->getStep(0.);
+  for (size_t i = 0; i<NDIM; ++i)
+    CPPUNIT_ASSERT( (size_t)0 == index[i]);
+  // check all possible shells on boundary
+  for (double length = EDGELENGTH; length < DOMAINLENGTH; length+=EDGELENGTH) {
+    LinkedCell::tripleIndex index = LC->getStep(length);
+    for (size_t i = 0; i<NDIM; ++i) {
+      std::cout << (size_t)(length/EDGELENGTH) << " ==" << index[i] << std::endl;
+      CPPUNIT_ASSERT( (LinkedCell::LinkedCellArray::index)(length/EDGELENGTH) == index[i]);
+    }
+  }
+  // check all possible shells at half interval
+  for (double length = 0.5 * EDGELENGTH; length < DOMAINLENGTH; length+=EDGELENGTH) {
+    LinkedCell::tripleIndex index = LC->getStep(length);
+    for (size_t i = 0; i<NDIM; ++i)
+      CPPUNIT_ASSERT( (LinkedCell::LinkedCellArray::index)ceil(length/EDGELENGTH) == index[i]);
+  }
+}
+
+/** UnitTest for insertPointCloud()
+ */
+void LinkedCell_ModelTest::getIndexToVectorTest()
+{
+  {
+    const Vector test(0.,0.,0.);
+    const LinkedCell::tripleIndex index = LC->getIndexToVector(test);
+    for (size_t i = 0; i<NDIM; ++i)
+      CPPUNIT_ASSERT( (size_t)0 == index[i]);
+  }
+  {
+    const Vector test(DOMAINLENGTH/2.,DOMAINLENGTH/2.,DOMAINLENGTH/2.);
+    const LinkedCell::tripleIndex index = LC->getIndexToVector(test);
+    for (size_t i = 0; i<NDIM; ++i)
+      CPPUNIT_ASSERT( (size_t)floor(DOMAINLENGTH/EDGELENGTH/2.) == index[i]);
+  }
+  {
+    const Vector test(DOMAINLENGTH/2.,DOMAINLENGTH/3.,DOMAINLENGTH/4.);
+    const LinkedCell::tripleIndex index = LC->getIndexToVector(test);
+    for (size_t i = 0; i<NDIM; ++i)
+      CPPUNIT_ASSERT( (LinkedCell::LinkedCellArray::index)floor(DOMAINLENGTH/EDGELENGTH/(double)(i+2.)) == index[i]);
+  }
+}
+
+/** UnitTest for insertPointCloud()
+ */
+void LinkedCell_ModelTest::nodeTest()
+{
+  // create point
+  TesselPoint * Walker = new TesselPoint();
+  Walker->setName(std::string("Walker9"));
+  Walker->setPosition(Vector(9.8,7.6,5.4));
+  PointCloudAdaptor< PointSet > cloud(&NodeList, std::string("NodeList"));
+  LC->insertPointCloud(cloud);
+
+  // check addNode
+  LC->addNode(Walker);
+  CPPUNIT_ASSERT_EQUAL( NodeList.size()+1, LC->CellLookup.size());
+  LinkedCell::tripleIndex index1 = LC->getIndexToVector(Walker->getPosition());
+  const LinkedCell::tripleIndex &index2 = LC->CellLookup[Walker]->getIndices();
+  CPPUNIT_ASSERT(index1 == index2);
+
+  // check moveNode
+#ifndef NDEBUG
+  std::cout << "The following assertion is intended and is not a failure of the code." << std::endl;
+  CPPUNIT_ASSERT_THROW( LC->moveNode(Walker), Assert::AssertionFailure );
+#endif
+
+  // check deleteNode
+  LC->deleteNode(Walker);
+  CPPUNIT_ASSERT_EQUAL( NodeList.size(), LC->CellLookup.size());
+
+  delete Walker;
+}
Index: src/LinkedCell/unittests/LinkedCell_ModelUnitTest.hpp
===================================================================
--- src/LinkedCell/unittests/LinkedCell_ModelUnitTest.hpp	(revision d8296152e374a74c483120eb411121989ba27a98)
+++ src/LinkedCell/unittests/LinkedCell_ModelUnitTest.hpp	(revision d8296152e374a74c483120eb411121989ba27a98)
@@ -0,0 +1,61 @@
+/*
+ * LinkedCell_ModelUnitTest.hpp
+ *
+ *  Created on: Nov 17, 2011
+ *      Author: heber
+ */
+
+#ifndef LINKEDCELL_MODELUNITTEST_HPP_
+#define LINKEDCELL_MODELUNITTEST_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+class Box;
+
+namespace LinkedCell {
+  class LinkedCell_Model;
+}
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/********************************************** Test classes **************************************/
+
+class LinkedCell_ModelTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( LinkedCell_ModelTest) ;
+    CPPUNIT_TEST ( AllocationTest );
+    CPPUNIT_TEST ( getSizeTest );
+    CPPUNIT_TEST ( ResetTest );
+    CPPUNIT_TEST ( insertPointCloudTest );
+    CPPUNIT_TEST ( setPartitionTest );
+    CPPUNIT_TEST ( getStepTest );
+    CPPUNIT_TEST ( getIndexToVectorTest );
+    CPPUNIT_TEST ( nodeTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+
+      void AllocationTest();
+      void getSizeTest();
+      void ResetTest();
+      void insertPointCloudTest();
+      void setPartitionTest();
+      void getStepTest();
+      void getIndexToVectorTest();
+      void nodeTest();
+
+private:
+      typedef std::set<TesselPoint *> PointSet;
+      PointSet NodeList;
+
+      Box *domain;
+      LinkedCell::LinkedCell_Model *LC;
+};
+
+
+#endif /* LINKEDCELL_MODELUNITTEST_HPP_ */
Index: src/LinkedCell/unittests/Makefile.am
===================================================================
--- src/LinkedCell/unittests/Makefile.am	(revision 02ce363b06a8f5ac471e055cfd494a84e197010f)
+++ src/LinkedCell/unittests/Makefile.am	(revision d8296152e374a74c483120eb411121989ba27a98)
@@ -4,14 +4,20 @@
 LINKEDCELLTESTSSOURCES = \
 	../LinkedCell/LinkedCell.cpp \
+	../LinkedCell/LinkedCell_Model.cpp \
 	../LinkedCell/unittests/LinkedCellUnitTest.cpp \
+	../LinkedCell/unittests/LinkedCell_ModelUnitTest.cpp \
 	../LinkedCell/unittests/linkedcellUnitTest.cpp
 
 LINKEDCELLTESTSHEADERS = \
 	../LinkedCell/LinkedCell.hpp \
+	../LinkedCell/LinkedCell_Model.hpp \
+	../LinkedCell/PointCloudAdaptor.hpp \
 	../LinkedCell/unittests/LinkedCellUnitTest.hpp \
+	../LinkedCell/unittests/LinkedCell_ModelUnitTest.hpp \
 	../LinkedCell/unittests/linkedcellUnitTest.hpp
 
 LINKEDCELLTESTS = \
   LinkedCellUnitTest \
+  LinkedCell_ModelUnitTest \
   linkedcellUnitTest
 
@@ -43,4 +49,13 @@
 	$(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la
 
+LinkedCell_ModelUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
+	../LinkedCell/unittests/LinkedCell_ModelUnitTest.cpp \
+	../LinkedCell/unittests/LinkedCell_ModelUnitTest.hpp \
+	../LinkedCell/LinkedCell.cpp \
+	../LinkedCell/LinkedCell_Model.cpp \
+	../LinkedCell/LinkedCell.hpp \
+	../LinkedCell/LinkedCell_Model.hpp \
+	../LinkedCell/PointCloudAdaptor.hpp
+LinkedCell_ModelUnitTest_LDADD = $(LINKEDCELLLIBS)
 
 
